How can I make this active snake contour more robust??

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

[Stelios Fanourakis]

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

[Stelios Fanourakis]

Anyone can answer me this??

[Jan]

[MOVED from section for answers] Stelios Fanourakis on 12 May 2019 at 0:22

What can I do if I am not getting answers to my questions?

[Jan]

@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.

[Stelios Fanourakis]

@Jan. Thanks for any of your contributions to my questions so far and I hope you continue to do so.

I disagree at what you say for not being the ideal place for me to get helped. I am indeed tremendously helped by many contributors here to almost all of my questions.

To the best of my knowledge, I provide enough and satisfacory information for my needs. For instance, in this case, more robust means to be more robust and ignore different shades of grey. The contour should be less sensitive to gradients and segments only the high gradient descends such as pure white and rest of grey tones. At least the pure white should be purely segmented and involve no other shade of grey. That's what robustness means to me.

Nevertheless, in case sometimes , I stop responding to some of my questions, it means that I have found a solution, either an alternative one or from any other source. Unless, I find a solution, I hope I keep posting here to this helpful and blissful forum.