Segment unfocused background pixels from foreground - Active contours

Hi - I have attached an image of some particles taken from a microscope. It is quite clear what is foreground and what is below the surface and makes up the blurry background. I would like to remove the background, or just isolate the foreground particles for further analysis. I have attempted active contour method (Chan-Vese), however this seems too sensitive and included the background areas as well as my desire particles. Is there a way to implement a less sensitive version of this?

Any advice would be greatly appreciated.

Thanks, Joe

Usage of varibles:
%   input: 
%       I           = any gray/double/RGB input image
%       mask        = initial mask, either customerlized or built-in
%       num_iter    = total number of iterations
%       mu          = weight of length term
%       method      = submethods pick from ('chen','vector','multiphase')
%
%   Types of built-in mask functions
%       'small'     = create a small circular mask
%       'medium'    = create a medium circular mask
%       'large'     = create a large circular mask
%       'whole'     = create a mask with holes around
%       'whole+small' = create a two layer mask with one layer small
%                       circular mask and the other layer with holes around
%                       (only work for method 'multiphase')
%   Types of methods
%       'chen'      = general CV method
%       'vector'    = CV method for vector image
%       'multiphase'= CV method for multiphase (2 phases applied here)
%
%   output: 
%       phi0        = updated level set function

          function seg = chenvese(I,mask,num_iter,mu,method)

              if true
                % code
              end
          %-- Default settings
          %   length term mu = 0.2 and default method = 'chan'
            if(~exist('mu','var')) 
              mu=0.2; 
            end

            if(~exist('method','var')) 
              method = 'chan'; 
            end

          %-- End default settings
          %-- Initializations on input image I and mask
          %  resize original image
             s = 200./min(size(I,1),size(I,2)); % resize scale
             if s<1
                 I = imresize(I,s);
             end

          %   auto mask settings
            if ischar(mask)
                switch lower (mask)
                    case 'small'
                        mask = maskcircle2(I,'small');
                    case 'medium'
                        mask = maskcircle2(I,'medium');
                    case 'large'
                        mask = maskcircle2(I,'large');              
                    case 'whole'
                        mask = maskcircle2(I,'whole'); 
                        %mask = init_mask(I,30);      
                    case 'whole+small'
                        m1 = maskcircle2(I,'whole');
                        m2 = maskcircle2(I,'small');
                        mask = zeros(size(I,1),size(I,2),2);
                        mask(:,:,1) = m1(:,:,1);
                        mask(:,:,2) = m2(:,:,2);
                    otherwise
                        error('unrecognized mask shape name (MASK).');
                end
            else
                if s<1
                    mask = imresize(mask,s);
                end
                if size(mask,1)>size(I,1) || size(mask,2)>size(I,2)
                    error('dimensions of mask unmathch those of the image.')
                end
                switch lower(method)
                    case 'multiphase'
                        if  (size(mask,3) == 1)  
                            error('multiphase requires two masks but only gets one.')
                        end
                end

            end       

          switch lower(method)
              case 'chan'
                  if size(I,3)== 3
                      P = rgb2gray(uint8(I));
                      P = double(P);
                  elseif size(I,3) == 2
                      P = 0.5.*(double(I(:,:,1))+double(I(:,:,2)));
                  else
                      P = double(I);
                  end
                  layer = 1;

              case 'vector'
                  s = 200./min(size(I,1),size(I,2)); % resize scale
                  I = imresize(I,s);
                  mask = imresize(mask,s);
                  layer = size(I,3);
                  if layer == 1
                      display('only one image component for vector image')
                  end
                  P = double(I);

              case 'multiphase'
                  layer = size(I,3);
                  if size(I,1)*size(I,2)>200^2
                      s = 200./min(size(I,1),size(I,2)); % resize scale
                      I = imresize(I,s);
                      mask = imresize(mask,s);
                  end

                  P = double(I);  %P store the original image
              otherwise
                  error('!invalid method')
          end
          %-- End Initializations on input image I and mask
          %--   Core function
          switch lower(method)
              case {'chan','vector'}
                  %-- SDF
                  %   Get the distance map of the initial mask

                  mask = mask(:,:,1);
                  phi0 = bwdist(mask)-bwdist(1-mask)+im2double(mask)-.5; 
                  %   initial force, set to eps to avoid division by zeros
                  force = eps; 
                  %-- End Initialization

                  %-- Display settings
                    figure();
                    subplot(2,2,1); imshow(I); title('Input Image');
                    subplot(2,2,2); contour(flipud(phi0), [0 0], 'r','LineWidth',1); title('initial contour');
                    subplot(2,2,3); title('Segmentation');
                  %-- End Display original image and mask

                  %-- Main loop
                    for n=1:num_iter
                        inidx = find(phi0>=0); % frontground index
                        outidx = find(phi0<0); % background index
                        force_image = 0; % initial image force for each layer 
                        for i=1:layer
                            L = im2double(P(:,:,i)); % get one image component
                            c1 = sum(sum(L.*Heaviside(phi0)))/(length(inidx)+eps); % average inside of Phi0
                            c2 = sum(sum(L.*(1-Heaviside(phi0))))/(length(outidx)+eps); % verage outside of Phi0
                            force_image=-(L-c1).^2+(L-c2).^2+force_image; 
                            % sum Image Force on all components (used for vector image)
                            % if 'chan' is applied, this loop become one sigle code as a
                            % result of layer = 1
                        end

                        % calculate the external force of the image 
                        force = mu*kappa(phi0)./max(max(abs(kappa(phi0))))+1/layer.*force_image;

                        % normalized the force
                        force = force./(max(max(abs(force))));

                        % get stepsize dt
                        dt=0.5;

                        % get parameters for checking whether to stop
                        old = phi0;
                        phi0 = phi0+dt.*force;
                        new = phi0;
                        indicator = checkstop(old,new,dt);

                        % intermediate output
                        if(mod(n,20) == 0) 
                           showphi(I,phi0,n);  
                        end;
                        if indicator % decide to stop or continue 
                            showphi(I,phi0,n);

                            %make mask from SDF
                            seg = phi0<=0; %-- Get mask from levelset

                            subplot(2,2,4); imshow(seg); title('Global Region-Based Segmentation');

                            return;
                        end
                    end;
                    showphi(I,phi0,n);

                    %make mask from SDF
                    seg = phi0<=0; %-- Get mask from levelset

                    subplot(2,2,4); imshow(seg); title('Global Region-Based Segmentation');
              case 'multiphase'
                  %-- Initializations
                  %   Get the distance map of the initial masks
                  mask1 = mask(:,:,1);
                  mask2 = mask(:,:,2);
                  phi1=bwdist(mask1)-bwdist(1-mask1)+im2double(mask1)-.5;%Get phi1 from the initial mask 1
                  phi2=bwdist(mask2)-bwdist(1-mask2)+im2double(mask2)-.5;%Get phi1 from the initial mask 2

                  %-- Display settings
                  figure();
                  subplot(2,2,1); 
                  if layer ~= 1
                      imshow(I); title('Input Image');
                  else
                      imagesc(P); axis image; colormap(gray);title('Input Image');
                  end
                  subplot(2,2,2);
                  hold on
                  contour(flipud(mask1),[0,0],'r','LineWidth',2.5); 
                  contour(flipud(mask1),[0,0],'x','LineWidth',1);
                  contour(flipud(mask2),[0,0],'g','LineWidth',2.5);
                  contour(flipud(mask2),[0,0],'x','LineWidth',1);
                  title('initial contour');
                  hold off
                  subplot(2,2,3); title('Segmentation');
                  %-- End display settings

                  %Main loop
                  for n=1:num_iter
                        %-- Narrow band for each phase
                        nb1 = find(phi1<1.2 & phi1>=-1.2); %narrow band of phi1
                        inidx1 = find(phi1>=0); %phi1 frontground index
                        outidx1 = find(phi1<0); %phi1 background index

                        nb2 = find(phi2<1.2 & phi2>=-1.2); %narrow band of phi2
                        inidx2 = find(phi2>=0); %phi2 frontground index
                        outidx2 = find(phi2<0); %phi2 background index
                        %-- End initiliazaions on narrow band

                        %-- Mean calculations for different partitions
                        %c11 = mean (phi1>0 & phi2>0)
                        %c12 = mean (phi1>0 & phi2<0)
                        %c21 = mean (phi1<0 & phi2>0)
                        %c22 = mean (phi1<0 & phi2<0)

                        cc11 = intersect(inidx1,inidx2); %index belong to (phi1>0 & phi2>0)
                        cc12 = intersect(inidx1,outidx2); %index belong to (phi1>0 & phi2<0)
                        cc21 = intersect(outidx1,inidx2); %index belong to (phi1<0 & phi2>0)
                        cc22 = intersect(outidx1,outidx2); %index belong to (phi1<0 & phi2<0)

                        f_image11 = 0;
                        f_image12 = 0;
                        f_image21 = 0;
                        f_image22 = 0; % initial image force for each layer 

                        for i=1:layer
                            L = im2double(P(:,:,i)); % get one image component

                        if isempty(cc11)
                            c11 = eps;
                        else
                            c11 = mean(L(cc11));
                        end

                        if isempty(cc12)
                            c12 = eps;
                        else
                            c12 = mean(L(cc12)); 
                        end

                        if isempty(cc21)
                            c21 = eps;
                        else
                            c21 = mean(L(cc21));
                        end

                        if isempty(cc22)
                            c22 = eps;
                        else
                            c22 = mean(L(cc22));
                        end

                        %-- End mean calculation

                        %-- Force calculation and normalization
                        % force on each partition

                        f_image11=(L-c11).^2.*Heaviside(phi1).*Heaviside(phi2)+f_image11;
                        f_image12=(L-c12).^2.*Heaviside(phi1).*(1-Heaviside(phi2))+f_image12;
                        f_image21=(L-c21).^2.*(1-Heaviside(phi1)).*Heaviside(phi2)+f_image21;
                        f_image22=(L-c22).^2.*(1-Heaviside(phi1)).*(1-Heaviside(phi2))+f_image22;
                        end

                            % sum Image Force on all components (used for vector image)
                            % if 'chan' is applied, this loop become one sigle code as a
                            % result of layer = 1

                        % calculate the external force of the image 

                        % curvature on phi1
                        curvature1 = mu*kappa(phi1);
                        curvature1 = curvature1(nb1);
                        % image force on phi1
                        fim1 = 1/layer.*(-f_image11(nb1)+f_image21(nb1)-f_image12(nb1)+f_image22(nb1));
                        fim1 = fim1./max(abs(fim1)+eps);

                        % curvature on phi2
                        curvature2 = mu*kappa(phi2);
                        curvature2 = curvature2(nb2);
                        % image force on phi2
                        fim2 = 1/layer.*(-f_image11(nb2)+f_image12(nb2)-f_image21(nb2)+f_image22(nb2));
                        fim2 = fim2./max(abs(fim2)+eps);

                        % force on phi1 and phi2
                        force1 = curvature1+fim1;
                        force2 = curvature2+fim2;
                        %-- End force calculation

                        % detal t
                        dt = 1.5;

                        old(:,:,1) = phi1;
                        old(:,:,2) = phi2;

                        %update of phi1 and phi2
                        phi1(nb1) = phi1(nb1)+dt.*force1;
                        phi2(nb2) = phi2(nb2)+dt.*force2;

                        new(:,:,1) = phi1;
                        new(:,:,2) = phi2;

                        indicator = checkstop(old,new,dt);

                        if indicator 
                           showphi(I, new, n);
                           %make mask from SDF
                           seg11 = (phi1>0 & phi2>0); %-- Get mask from levelset
                           seg12 = (phi1>0 & phi2<0);
                           seg21 = (phi1<0 & phi2>0);
                           seg22 = (phi1<0 & phi2<0);

                           se = strel('disk',1);
                           aa1 = imerode(seg11,se);
                           aa2 = imerode(seg12,se);
                           aa3 = imerode(seg21,se);
                           aa4 = imerode(seg22,se);
                           seg = aa1+2*aa2+3*aa3+4*aa4;

                           subplot(2,2,4); imagesc(seg);axis image;title('Global Region-Based Segmentation');

                            return
                        end
                        % re-initializations
                        phi1 = reinitialization(phi1, 0.6);%sussman(phi1, 0.6);%
                        phi2 = reinitialization(phi2, 0.6);%sussman(phi2,0.6);

                        %intermediate output
                        if(mod(n,20) == 0) 
                           phi(:,:,1) = phi1;
                           phi(:,:,2) = phi2;
                           showphi(I, phi, n);
                        end;
                    end;
                    phi(:,:,1) = phi1;
                    phi(:,:,2) = phi2;
                    showphi(I, phi, n);
                    %make mask from SDF
                  seg11 = (phi1>0 & phi2>0); %-- Get mask from levelset
                  seg12 = (phi1>0 & phi2<0);
                  seg21 = (phi1<0 & phi2>0);
                  seg22 = (phi1<0 & phi2<0);

                  se = strel('disk',1);
                  aa1 = imerode(seg11,se);
                  aa2 = imerode(seg12,se);
                  aa3 = imerode(seg21,se);
                  aa4 = imerode(seg22,se);
                  seg = aa1+2*aa2+3*aa3+4*aa4;
                  %seg = bwlabel(seg);
                  subplot(2,2,4); imagesc(seg);axis image;title('Global Region-Based Segmentation');
          end