C =updateC(Img, u, KB1, KB2, epsilon);

what is the meaning of updateC?

In the last of the program they have given

% Make a function satisfy Neumann boundary condition

if possible anyone please help me by explaining it.

function [phi, b, C]= lse_bfe(phi0,Img, b, Ksigma,KONE, v,timestep,mu,epsilon, iter_lse)
phi=phi0;
KB1 = conv2(b,K,'same');
KB2 = conv2(b.^2,K,'same');
C =updateC(Img, phi, KB1, KB2);

KONE_Img = Img.^2.*KONE;
phi = updateLSF(Img,phi, C, KONE_Img, KB1, KB2, mu, nu, timestep, epsilon, iter_lse);

Hphi=Heaviside(phi,epsilon);
M(:,:,1)=Hphi;
M(:,:,2)=1-Hphi;
b =updateB(Img, C, M,  Ksigma);

% update level set function
function phi = updateLSF(Img, phi0, C, KONE_Img, KB1, KB2, mu, v, timestep, epsilon, iter_lse)
phi=phi0;
Hphi=Heaviside(phi,epsilon);
M(:,:,1)=Hphi;
M(:,:,2)=1-Hphi;
N_class=size(M,3);
e=zeros(size(M));

for k=1:N_class
    e(:,:,k) = KONE_Img - 2*Img.*C(k).*KB1 + C(k)^2*KB2;
end

fork=1:iter_lse
    phi=NeumannBoundCond(phi);
    K=curvature_central(phi);    % div()
    DiracPHI=Dirac(phi,epsilon);
    ImageTerm=-DiracPHI.*(e(:,:,1)-e(:,:,2));
    penalizeTerm=mu*(4*del2(phi)-K);
    lengthTerm=v.*DiracPHI.*K;
    phi=phi+timestep*(lengthTerm+penalizeTerm+ImageTerm);
end

% update b
function b =updateB(Img, C, M,  Ksigma)

PC1=zeros(size(Img));
PC2=PC1;
N_class=size(M,3);
for kk=1:N_class
    PC1=PC1+C(kk)*M(:,:,kk);
    PC2=PC2+C(kk)^2*M(:,:,kk);
end
KNm1 = conv2(PC1.*Img,Ksigma,'same');
KDn1 = conv2(PC2,Ksigma,'same');

b = KNm1./KDn1;

% Update C
function C_new =updateC(Img, phi, Kb1, Kb2, epsilon)
Hu=Heaviside(phi,epsilon);
M(:,:,1)=Hu;
M(:,:,2)=1-Hu;
N_class=size(M,3);
for kk=1:N_class
    Nm2 = Kb1.*Img.*M(:,:,kk);
    Dn2 = Kb2.*M(:,:,kk);
    C_new(kk) = sum(Nm2(:))/sum(Dn2(:));
end

% Make a function satisfy Neumann boundary condition
function g = NeumannBoundCond(f)
[nrow,ncol] = size(f);
g = f;
g([1 nrow],[1 ncol]) = g([3 nrow-2],[3 ncol-2]);  
g([1 nrow],2:end-1) = g([3 nrow-2],2:end-1);          
g(2:end-1,[1 ncol]) = g(2:end-1,[3 ncol-2]);  

function k = curvature_central(phi)
% compute curvature for phi with central difference scheme
[phix,phiy] = gradient(phi);
normDphi = sqrt(phix.^2+phiy.^2+1e-10);
Nx = phix./normDphi;
Ny = phiy./normDphi;
[nxx,junk] = gradient(Nx);
[junk,nyy] = gradient(Ny);
k = nxx+nyy;

function h = Heaviside(x,epsilon)    
h=0.5*(1+(2/pi)*atan(x./epsilon));

function f = Dirac(x, epsilon)    
f=(epsilon/pi)./(epsilon^2.+x.^2);