source: trunk/src/transform_field/ima_find_particles.m @ 671

Last change on this file since 671 was 671, checked in by sommeria, 10 years ago

geometry_calib corrected for translation and rotation + cleaning

File size: 4.1 KB
Line 
1% 'ima_remove_particles': removes particles from an image (keeping the local minimum)
2% requires the Matlab image processing toolbox
3
4%------------------------------------------------------------------------
5%%%%  Use the general syntax for transform fields with a single input %%%%
6% OUTPUT:
7% DataOut:   output field structure
8
9%INPUT:
10% DataIn:  first input field structure
11%------------------------------------------------------------------------
12function DataOut=ima_find_particles(DataIn)
13%------------------------------------------------------------------------
14DataOut=DataIn;  %default  output field
15if strcmp(DataIn,'*')
16    return
17end
18
19%parameters
20AbsThreshold=150;
21SizePart=3;
22%---------------------------------------------------------
23   %A=double(DataIn.A(:,:,3));% take the blue component
24% if ndims(DataIn.A)==3;%color images
25    A=sum(double(DataIn.A),3);% take the sum of color components
26% end
27%%  mask to reduce the  working area (optional)
28Mask=ones(size(A));
29Mask(1:SizePart,:)=0;
30Mask(end-SizePart:end,:)=0;
31Mask(:,1:SizePart)=0;
32Mask(:,end-SizePart:end)=0;
33[Js,Is]=find(A<AbsThreshold &abs(double(DataIn.A(:,:,1))-double(DataIn.A(:,:,3)))<20 & Mask==1);%indices (I,J) of dark pixels
34X=zeros(size(Is));
35Y=zeros(size(Js));
36F=zeros(size(Js));
37for ipart=1:numel(Is)
38    if Mask(Js(ipart),Is(ipart))==1
39        subimage=A(Js(ipart)-SizePart:Js(ipart)+SizePart,Is(ipart)-SizePart:Is(ipart)+SizePart);
40        subimage=max(max(subimage))-subimage;%take negative of the image
41        [vector,F(ipart)] = SUBPIX2DGAUSS (subimage,SizePart+1,SizePart+1);
42        %             X0(ipart)=Is(ipart);%TEST
43        %             Y0(ipart)=Js(ipart);%TEST
44        X(ipart)=Is(ipart)+vector(1);%corrected position
45        Y(ipart)=Js(ipart)+vector(2);
46        Xround=round(X(ipart));
47        Xlow=max(1,Xround-SizePart);
48        Xhigh=min(size(A,2),Xround+SizePart);
49        Yround=round(Y(ipart));
50        Ylow=max(1,Yround-SizePart);
51        Yhigh=min(size(A,1),Yround+SizePart);
52        Mask(Ylow:Yhigh,Xlow:Xhigh)=0;% mask the subregion already treated to
53        % avoid double counting
54    end
55end
56X=X(X>0);
57Y=Y(Y>0);
58huvmat=findobj(allchild(0),'Tag','uvmat');
59if ~isempty(huvmat)
60    haxes=findobj(huvmat,'Tag','PlotAxes');
61    set(haxes,'NextPlot','add')
62    % hold on
63    axes(haxes)
64    plot(X-0.5,size(A,1)-Y+0.5,'+')
65    set(haxes,'NextPlot','replace')
66end
67% hold off
68hmovie=findobj(allchild(0),'Tag','movieaxes');
69if ~isempty(hmovie)
70    set(hmovie,'NextPlot','add')
71    axes(hmovie)
72    plot(X-0.5,size(A,1)-Y+0.5,'+')
73    set(hmovies,'NextPlot','replace')
74end
75
76       
77        %------------------------------------------------------------------------
78% --- Find the maximum of the correlation function after interpolation
79function [vector,F] = SUBPIX2DGAUSS (result_conv,x,y)
80%------------------------------------------------------------------------
81vector=[0 0]; %default
82F=-2;
83peaky=y;
84peakx=x;
85[npy,npx]=size(result_conv);
86if (x <= npx-1) && (y <= npy-1) && (x >= 1) && (y >= 1)
87    F=0;
88    for i=-1:1
89        for j=-1:1
90            %following 15 lines based on
91            %H. Nobach ï¿œ M. Honkanen (2005)
92            %Two-dimensional Gaussian regression for sub-pixel displacement
93            %estimation in particle image velocimetry or particle position
94            %estimation in particle tracking velocimetry
95            %Experiments in Fluids (2005) 38: 511ï¿œ515
96            c10(j+2,i+2)=i*log(result_conv(y+j, x+i));
97            c01(j+2,i+2)=j*log(result_conv(y+j, x+i));
98            c11(j+2,i+2)=i*j*log(result_conv(y+j, x+i));
99            c20(j+2,i+2)=(3*i^2-2)*log(result_conv(y+j, x+i));
100            c02(j+2,i+2)=(3*j^2-2)*log(result_conv(y+j, x+i));
101        end
102    end
103    c10=(1/6)*sum(sum(c10));
104    c01=(1/6)*sum(sum(c01));
105    c11=(1/4)*sum(sum(c11));
106    c20=(1/6)*sum(sum(c20));
107    c02=(1/6)*sum(sum(c02));
108    deltax=(c11*c01-2*c10*c02)/(4*c20*c02-c11^2);
109    deltay=(c11*c10-2*c01*c20)/(4*c20*c02-c11^2);
110    if abs(deltax)<1
111        peakx=x+deltax;
112    end
113    if abs(deltay)<1
114        peaky=y+deltay;
115    end
116end
117vector=[peakx-floor(npx/2)-1 peaky-floor(npy/2)-1];
Note: See TracBrowser for help on using the repository browser.