Index: trunk/src/transform_field/ima_find_particles.m
===================================================================
--- trunk/src/transform_field/ima_find_particles.m	(revision 671)
+++ trunk/src/transform_field/ima_find_particles.m	(revision 671)
@@ -0,0 +1,117 @@
+% 'ima_remove_particles': removes particles from an image (keeping the local minimum)
+% requires the Matlab image processing toolbox
+
+%------------------------------------------------------------------------
+%%%%  Use the general syntax for transform fields with a single input %%%%
+% OUTPUT: 
+% DataOut:   output field structure 
+
+%INPUT:
+% DataIn:  first input field structure
+%------------------------------------------------------------------------
+function DataOut=ima_find_particles(DataIn)
+%------------------------------------------------------------------------
+DataOut=DataIn;  %default  output field
+if strcmp(DataIn,'*')
+    return
+end
+
+%parameters
+AbsThreshold=150;
+SizePart=3;
+%--------------------------------------------------------- 
+   %A=double(DataIn.A(:,:,3));% take the blue component
+% if ndims(DataIn.A)==3;%color images
+    A=sum(double(DataIn.A),3);% take the sum of color components
+% end
+%%  mask to reduce the  working area (optional)
+Mask=ones(size(A));
+Mask(1:SizePart,:)=0;
+Mask(end-SizePart:end,:)=0;
+Mask(:,1:SizePart)=0;
+Mask(:,end-SizePart:end)=0;
+[Js,Is]=find(A<AbsThreshold &abs(double(DataIn.A(:,:,1))-double(DataIn.A(:,:,3)))<20 & Mask==1);%indices (I,J) of dark pixels
+X=zeros(size(Is));
+Y=zeros(size(Js));
+F=zeros(size(Js));
+for ipart=1:numel(Is)
+    if Mask(Js(ipart),Is(ipart))==1
+        subimage=A(Js(ipart)-SizePart:Js(ipart)+SizePart,Is(ipart)-SizePart:Is(ipart)+SizePart);
+        subimage=max(max(subimage))-subimage;%take negative of the image
+        [vector,F(ipart)] = SUBPIX2DGAUSS (subimage,SizePart+1,SizePart+1);
+        %             X0(ipart)=Is(ipart);%TEST
+        %             Y0(ipart)=Js(ipart);%TEST
+        X(ipart)=Is(ipart)+vector(1);%corrected position
+        Y(ipart)=Js(ipart)+vector(2);
+        Xround=round(X(ipart));
+        Xlow=max(1,Xround-SizePart);
+        Xhigh=min(size(A,2),Xround+SizePart);
+        Yround=round(Y(ipart));
+        Ylow=max(1,Yround-SizePart);
+        Yhigh=min(size(A,1),Yround+SizePart);
+        Mask(Ylow:Yhigh,Xlow:Xhigh)=0;% mask the subregion already treated to
+        % avoid double counting
+    end
+end
+X=X(X>0);
+Y=Y(Y>0);
+huvmat=findobj(allchild(0),'Tag','uvmat');
+if ~isempty(huvmat)
+    haxes=findobj(huvmat,'Tag','PlotAxes');
+    set(haxes,'NextPlot','add')
+    % hold on
+    axes(haxes)
+    plot(X-0.5,size(A,1)-Y+0.5,'+')
+    set(haxes,'NextPlot','replace')
+end
+% hold off
+hmovie=findobj(allchild(0),'Tag','movieaxes');
+if ~isempty(hmovie)
+    set(hmovie,'NextPlot','add')
+    axes(hmovie)
+    plot(X-0.5,size(A,1)-Y+0.5,'+')
+    set(hmovies,'NextPlot','replace')
+end
+
+        
+        %------------------------------------------------------------------------
+% --- Find the maximum of the correlation function after interpolation
+function [vector,F] = SUBPIX2DGAUSS (result_conv,x,y)
+%------------------------------------------------------------------------
+vector=[0 0]; %default
+F=-2;
+peaky=y;
+peakx=x;
+[npy,npx]=size(result_conv);
+if (x <= npx-1) && (y <= npy-1) && (x >= 1) && (y >= 1)
+    F=0;
+    for i=-1:1
+        for j=-1:1
+            %following 15 lines based on
+            %H. Nobach ￜ M. Honkanen (2005)
+            %Two-dimensional Gaussian regression for sub-pixel displacement
+            %estimation in particle image velocimetry or particle position
+            %estimation in particle tracking velocimetry
+            %Experiments in Fluids (2005) 38: 511ￜ515
+            c10(j+2,i+2)=i*log(result_conv(y+j, x+i));
+            c01(j+2,i+2)=j*log(result_conv(y+j, x+i));
+            c11(j+2,i+2)=i*j*log(result_conv(y+j, x+i));
+            c20(j+2,i+2)=(3*i^2-2)*log(result_conv(y+j, x+i));
+            c02(j+2,i+2)=(3*j^2-2)*log(result_conv(y+j, x+i));
+        end
+    end
+    c10=(1/6)*sum(sum(c10));
+    c01=(1/6)*sum(sum(c01));
+    c11=(1/4)*sum(sum(c11));
+    c20=(1/6)*sum(sum(c20));
+    c02=(1/6)*sum(sum(c02)); 
+    deltax=(c11*c01-2*c10*c02)/(4*c20*c02-c11^2);
+    deltay=(c11*c10-2*c01*c20)/(4*c20*c02-c11^2);
+    if abs(deltax)<1
+        peakx=x+deltax;
+    end
+    if abs(deltay)<1
+        peaky=y+deltay;
+    end
+end
+vector=[peakx-floor(npx/2)-1 peaky-floor(npy/2)-1];