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

2 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

Jan
on 17 May 2019

@Stelios: As fas as I can see, none of your many questions contained enough information to be answered. You have ignored many questions for clarifications and it seems like you ignore the explanations given as answers. I assume, many contributors of this forum do not care about your questions anymore, because it was not possible to help you in the past. This would mean, that the forum is simply not a suitable place for you to learn Matlab.

In this case that vague expression "more robust" is not clear enought and you should explain statements like "I know it's a matter of tuning some optimization parameters" with any details.

I has been suggested to you to read How to ask a good question, but as far as I remember you did not react to this useful advice yet.

It is my opinion, that if the participation in a forum does not help you to learn Matlab, it is the best strategy to try it somewhere else.

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

Start Hunting!