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

Last change on this file since 924 was 924, checked in by g7moreau, 9 years ago
  • Update Copyright Copyright 2008-2016, LEGI UMR 5519 / CNRS UGA G-INP, Grenoble, France
File size: 4.9 KB
RevLine 
[671]1% 'ima_remove_particles': removes particles from an image (keeping the local minimum)
2% requires the Matlab image processing toolbox
3%------------------------------------------------------------------------
4%%%%  Use the general syntax for transform fields with a single input %%%%
5% OUTPUT:
6% DataOut:   output field structure
[810]7%
[671]8%INPUT:
9% DataIn:  first input field structure
[810]10
11%=======================================================================
[924]12% Copyright 2008-2016, LEGI UMR 5519 / CNRS UGA G-INP, Grenoble, France
[810]13%   http://www.legi.grenoble-inp.fr
14%   Joel.Sommeria - Joel.Sommeria (A) legi.cnrs.fr
15%
16%     This file is part of the toolbox UVMAT.
17%
18%     UVMAT is free software; you can redistribute it and/or modify
19%     it under the terms of the GNU General Public License as published
20%     by the Free Software Foundation; either version 2 of the license,
21%     or (at your option) any later version.
22%
23%     UVMAT is distributed in the hope that it will be useful,
24%     but WITHOUT ANY WARRANTY; without even the implied warranty of
25%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
26%     GNU General Public License (see LICENSE.txt) for more details.
27%=======================================================================
28
[671]29function DataOut=ima_find_particles(DataIn)
30%------------------------------------------------------------------------
31DataOut=DataIn;  %default  output field
32if strcmp(DataIn,'*')
33    return
34end
35
36%parameters
37AbsThreshold=150;
38SizePart=3;
39%---------------------------------------------------------
40   %A=double(DataIn.A(:,:,3));% take the blue component
41% if ndims(DataIn.A)==3;%color images
42    A=sum(double(DataIn.A),3);% take the sum of color components
43% end
44%%  mask to reduce the  working area (optional)
45Mask=ones(size(A));
46Mask(1:SizePart,:)=0;
47Mask(end-SizePart:end,:)=0;
48Mask(:,1:SizePart)=0;
49Mask(:,end-SizePart:end)=0;
50[Js,Is]=find(A<AbsThreshold &abs(double(DataIn.A(:,:,1))-double(DataIn.A(:,:,3)))<20 & Mask==1);%indices (I,J) of dark pixels
51X=zeros(size(Is));
52Y=zeros(size(Js));
53F=zeros(size(Js));
54for ipart=1:numel(Is)
55    if Mask(Js(ipart),Is(ipart))==1
56        subimage=A(Js(ipart)-SizePart:Js(ipart)+SizePart,Is(ipart)-SizePart:Is(ipart)+SizePart);
57        subimage=max(max(subimage))-subimage;%take negative of the image
58        [vector,F(ipart)] = SUBPIX2DGAUSS (subimage,SizePart+1,SizePart+1);
59        %             X0(ipart)=Is(ipart);%TEST
60        %             Y0(ipart)=Js(ipart);%TEST
61        X(ipart)=Is(ipart)+vector(1);%corrected position
62        Y(ipart)=Js(ipart)+vector(2);
63        Xround=round(X(ipart));
64        Xlow=max(1,Xround-SizePart);
65        Xhigh=min(size(A,2),Xround+SizePart);
66        Yround=round(Y(ipart));
67        Ylow=max(1,Yround-SizePart);
68        Yhigh=min(size(A,1),Yround+SizePart);
69        Mask(Ylow:Yhigh,Xlow:Xhigh)=0;% mask the subregion already treated to
70        % avoid double counting
71    end
72end
73X=X(X>0);
74Y=Y(Y>0);
75huvmat=findobj(allchild(0),'Tag','uvmat');
76if ~isempty(huvmat)
77    haxes=findobj(huvmat,'Tag','PlotAxes');
78    set(haxes,'NextPlot','add')
79    % hold on
80    axes(haxes)
81    plot(X-0.5,size(A,1)-Y+0.5,'+')
82    set(haxes,'NextPlot','replace')
83end
84% hold off
85hmovie=findobj(allchild(0),'Tag','movieaxes');
86if ~isempty(hmovie)
87    set(hmovie,'NextPlot','add')
88    axes(hmovie)
89    plot(X-0.5,size(A,1)-Y+0.5,'+')
90    set(hmovies,'NextPlot','replace')
91end
92
93       
94        %------------------------------------------------------------------------
95% --- Find the maximum of the correlation function after interpolation
96function [vector,F] = SUBPIX2DGAUSS (result_conv,x,y)
97%------------------------------------------------------------------------
98vector=[0 0]; %default
99F=-2;
100peaky=y;
101peakx=x;
102[npy,npx]=size(result_conv);
103if (x <= npx-1) && (y <= npy-1) && (x >= 1) && (y >= 1)
104    F=0;
105    for i=-1:1
106        for j=-1:1
107            %following 15 lines based on
108            %H. Nobach ï¿œ M. Honkanen (2005)
109            %Two-dimensional Gaussian regression for sub-pixel displacement
110            %estimation in particle image velocimetry or particle position
111            %estimation in particle tracking velocimetry
112            %Experiments in Fluids (2005) 38: 511ï¿œ515
113            c10(j+2,i+2)=i*log(result_conv(y+j, x+i));
114            c01(j+2,i+2)=j*log(result_conv(y+j, x+i));
115            c11(j+2,i+2)=i*j*log(result_conv(y+j, x+i));
116            c20(j+2,i+2)=(3*i^2-2)*log(result_conv(y+j, x+i));
117            c02(j+2,i+2)=(3*j^2-2)*log(result_conv(y+j, x+i));
118        end
119    end
120    c10=(1/6)*sum(sum(c10));
121    c01=(1/6)*sum(sum(c01));
122    c11=(1/4)*sum(sum(c11));
123    c20=(1/6)*sum(sum(c20));
124    c02=(1/6)*sum(sum(c02));
125    deltax=(c11*c01-2*c10*c02)/(4*c20*c02-c11^2);
126    deltay=(c11*c10-2*c01*c20)/(4*c20*c02-c11^2);
127    if abs(deltax)<1
128        peakx=x+deltax;
129    end
130    if abs(deltay)<1
131        peaky=y+deltay;
132    end
133end
134vector=[peakx-floor(npx/2)-1 peaky-floor(npy/2)-1];
Note: See TracBrowser for help on using the repository browser.