I guess it must be this line but am not sure how to modify it to be more robust

dphidt = F./max(abs(F)) + alpha*curvature; % gradient descent to minimize energy

9 views (last 30 days)

Show older comments

Hello

I am using a ChanVese active contour from the MathWorks community. I need to make it more robust. I need it to segment only the white sharp lines from the background and not recognise any other shade of gray in an ultrasound image. This way, will help me solve the artifacts issues.

I know it's a matter of tuning some optimization parameters but I don't know which of those. Can you please help?

% Region Based Active Contour Segmentation

%

% seg = region_seg(I,init_mask,max_its,alpha,display)

%

% Inputs: I 2D image

% init_mask Initialization (1 = foreground, 0 = bg)

% max_its Number of iterations to run segmentation for

% alpha (optional) Weight of smoothing term

% higer = smoother. default = 0.2

% display (optional) displays intermediate outputs

% default = true

%

% Outputs: seg Final segmentation mask (1=fg, 0=bg)

%

% Description: This code implements the paper: "Active Contours Without

% Edges" By Chan Vese. This is a nice way to segment images whose

% foregrounds and backgrounds are statistically different and homogeneous.

%

% Example:

% img = imread('tire.tif');

% m = zeros(size(img));

% m(33:33+117,44:44+128) = 1;

% seg = region_seg(img,m,500);

%

% Coded by: Shawn Lankton (www.shawnlankton.com)

%------------------------------------------------------------------------

function seg = 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,'initialmagnification',200,'displayrange',[0 255]); hold on;

contour(phi, [0 0], 'g','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);

active

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

Start Hunting!