--- a
+++ b/functions/ac_seg.m
@@ -0,0 +1,198 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% I : input image
+% init_mask: initial mask
+% max_its : maximum iterations
+% lengthEweight : weight of the length energy 
+% shapeEweight  : weight of the shape energy
+% display : display 
+
+% This code was edited for LV segmentation in MRI from the original code by Shawn Lankton
+
+function [seg,phi] = ac_seg(I,init_mask,max_its,lengthEweight,shapeEweight,display)
+
+  %-- 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);
+  
+  phi_prior=phi;
+  
+  %--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; % region-based force from image information
+    curvature = get_curvature(phi,idx);  % force from curvature penalty
+    
+    % derivite of energy function for prior shape
+    dEdl= -(phi(idx)-phi_prior(idx));
+    
+    % note that d phi/dt= - d E/d phi
+    dphidt = F./max(abs(F)) + lengthEweight*curvature+shapeEweight*dEdl;  % gradient descent to minimize energy
+    
+    %-- maintain the CFL condition
+    dt = .45/(max(dphidt)+eps);
+    
+    %-- evolve the curve
+    phi(idx) = phi(idx) + dt*dphidt;
+
+    indicator=0;
+     
+    %-- Keep SDF smooth
+    phi = sussman(phi, .5);
+
+    %-- intermediate output
+    if((display>0)&&(mod(its,20) == 0)) 
+      showCurveAndPhi(I,phi,its);  
+    end
+    if indicator==1
+        disp(['stopped at:',num2str(its)])
+        break
+    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)
+   max_range=min(200,max(I(:)));
+  imshow(I,'initialmagnification',200,'displayrange',[0 max_range]); 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);    
+
+  
+
+
+
+