can u pls explain the code in detail?

2 views (last 30 days)
ASHA
ASHA on 7 Apr 2014
Moved: DGM on 12 Feb 2023
function [seg phi] = region_seg(I,init_mask,max_its,alpha,display)
%-- default value for parameter alpha is .1
if(~exist('alpha','var'))
alpha = .2;
end
%-- default behavior is to display intermediate outputs
if(~exist('display','var'))
display = true;
end
%-- ensures image is 2D double matrix
I = im2graydouble(I);
%-- Create a signed distance map (SDF) from mask
phi = mask2phi(init_mask);
%--main loop
for its = 1:max_its % Note: no automatic convergence test
idx = find(phi <= 1.2 & phi >= -1.2); %get the curve's narrow band
%-- find interior and exterior mean
upts = find(phi<=0); % interior points
vpts = find(phi>0); % exterior points
u = sum(I(upts))/(length(upts)+eps); % interior mean
v = sum(I(vpts))/(length(vpts)+eps); % exterior mean
F = (I(idx)-u).^2-(I(idx)-v).^2; % force from image information
curvature = get_curvature(phi,idx); % force from curvature penalty
dphidt = F./max(abs(F)) + alpha*curvature; % gradient descent to minimize energy
%-- maintain the CFL condition
dt = .45/(max(dphidt)+eps);
%-- evolve the curve
phi(idx) = phi(idx) + dt.*dphidt;
%-- Keep SDF smooth
phi = sussman(phi, .5);
%-- intermediate output
if((display>0)&&(mod(its,20) == 0))
showCurveAndPhi(I,phi,its);
end
end
%-- final output
if(display)
showCurveAndPhi(I,phi,its);
end
%-- make mask from SDF
seg = phi<=0; %-- Get mask from levelset
%---------------------------------------------------------------------
%---------------------------------------------------------------------
%-- AUXILIARY FUNCTIONS ----------------------------------------------
%---------------------------------------------------------------------
%---------------------------------------------------------------------
%-- Displays the image with curve superimposed
function showCurveAndPhi(I, phi, i)
imshow(I); hold on;
contour(phi, [0 0], 'g','LineWidth',4);
contour(phi, [0 0], 'k','LineWidth',2);
hold off; title([num2str(i) ' Iterations']); drawnow;
%-- converts a mask to a SDF
function phi = mask2phi(init_a)
phi=bwdist(init_a)-bwdist(1-init_a)+im2double(init_a)-.5;
%-- compute curvature along SDF
function curvature = get_curvature(phi,idx)
[dimy, dimx] = size(phi);
[y x] = ind2sub([dimy,dimx],idx); % get subscripts
%-- get subscripts of neighbors
ym1 = y-1; xm1 = x-1; yp1 = y+1; xp1 = x+1;
%-- bounds checking
ym1(ym1<1) = 1; xm1(xm1<1) = 1;
yp1(yp1>dimy)=dimy; xp1(xp1>dimx) = dimx;
%-- get indexes for 8 neighbors
idup = sub2ind(size(phi),yp1,x);
iddn = sub2ind(size(phi),ym1,x);
idlt = sub2ind(size(phi),y,xm1);
idrt = sub2ind(size(phi),y,xp1);
idul = sub2ind(size(phi),yp1,xm1);
idur = sub2ind(size(phi),yp1,xp1);
iddl = sub2ind(size(phi),ym1,xm1);
iddr = sub2ind(size(phi),ym1,xp1);
%-- get central derivatives of SDF at x,y
phi_x = -phi(idlt)+phi(idrt);
phi_y = -phi(iddn)+phi(idup);
phi_xx = phi(idlt)-2*phi(idx)+phi(idrt);
phi_yy = phi(iddn)-2*phi(idx)+phi(idup);
phi_xy = -0.25*phi(iddl)-0.25*phi(idur)...
+0.25*phi(iddr)+0.25*phi(idul);
phi_x2 = phi_x.^2;
phi_y2 = phi_y.^2;
%-- compute curvature (Kappa)
curvature = ((phi_x2.*phi_yy + phi_y2.*phi_xx - 2*phi_x.*phi_y.*phi_xy)./...
(phi_x2 + phi_y2 +eps).^(3/2)).*(phi_x2 + phi_y2).^(1/2);
%-- Converts image to one channel (grayscale) double
function img = im2graydouble(img)
[dimy, dimx, c] = size(img);
if(isfloat(img)) % image is a double
if(c==3)
img = rgb2gray(uint8(img));
end
else % image is a int
if(c==3)
img = rgb2gray(img);
end
img = double(img);
end
%-- level set re-initialization by the sussman method
function D = sussman(D, dt)
% forward/backward differences
a = D - shiftR(D); % backward
b = shiftL(D) - D; % forward
c = D - shiftD(D); % backward
d = shiftU(D) - D; % forward
a_p = a; a_n = a; % a+ and a-
b_p = b; b_n = b;
c_p = c; c_n = c;
d_p = d; d_n = d;
a_p(a < 0) = 0;
a_n(a > 0) = 0;
b_p(b < 0) = 0;
b_n(b > 0) = 0;
c_p(c < 0) = 0;
c_n(c > 0) = 0;
d_p(d < 0) = 0;
d_n(d > 0) = 0;
dD = zeros(size(D));
D_neg_ind = find(D < 0);
D_pos_ind = find(D > 0);
dD(D_pos_ind) = sqrt(max(a_p(D_pos_ind).^2, b_n(D_pos_ind).^2) ...
+ max(c_p(D_pos_ind).^2, d_n(D_pos_ind).^2)) - 1;
dD(D_neg_ind) = sqrt(max(a_n(D_neg_ind).^2, b_p(D_neg_ind).^2) ...
+ max(c_n(D_neg_ind).^2, d_p(D_neg_ind).^2)) - 1;
D = D - dt .* sussman_sign(D) .* dD;
%-- whole matrix derivatives
function shift = shiftD(M)
shift = shiftR(M')';
function shift = shiftL(M)
shift = [ M(:,2:size(M,2)) M(:,size(M,2)) ];
function shift = shiftR(M)
shift = [ M(:,1) M(:,1:size(M,2)-1) ];
function shift = shiftU(M)
shift = shiftL(M')';
function S = sussman_sign(D)
S = D ./ sqrt(D.^2 + 1);
  1 Comment
Jan
Jan on 7 Apr 2014
@ASHA: Come on, you can imagine that a "detailed explanation" of the complete code is a rather general question. It is impossible to be answered efficiently. Please take the time to narrow down your problem and ask specific questions.

Sign in to comment.

Answers (2)

Walter Roberson
Walter Roberson on 7 Apr 2014
My charge to explain code "in detail" is a minimum of $US450 per hour (for code I agree to work on at all). I estimate there is at least 3 1/2 hours of work in explaining "in detail" what you have posted. More if you want to know why the sussman() function is the way it is instead of just being told what the mechanical steps are that are used within it.
  3 Comments
mohd abdul wahed faisal faisal
Moved: DGM on 12 Feb 2023
for that code that much charge is not worth.
Walter Roberson
Walter Roberson on 15 Jul 2019
Moved: DGM on 12 Feb 2023
Perhaps you have a different understanding of "in detail" than we do. When someone asks to explain code "in detail", I have to assume that I would have to start back with the history of computing, and the theory of electronic and mechanical computing, and work up to specifics about the code in question, with stops along the way on Group Theory, and the Peano Postulates to explain multiplication and addition on computers -- material enough for several textbooks.

Sign in to comment.


lj zhang
lj zhang on 14 Feb 2018
to see the paper Chan & Vese, Active Contours Without Edges, TIP 2001, and you will know anything. It is free of charge.

Products

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!