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
|
---|
7 | %
|
---|
8 | %INPUT:
|
---|
9 | % DataIn: first input field structure
|
---|
10 |
|
---|
11 | %=======================================================================
|
---|
12 | % Copyright 2008-2016, LEGI UMR 5519 / CNRS UGA G-INP, Grenoble, France
|
---|
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 |
|
---|
29 | function DataOut=ima_find_particles(DataIn)
|
---|
30 | %------------------------------------------------------------------------
|
---|
31 | DataOut=DataIn; %default output field
|
---|
32 | if strcmp(DataIn,'*')
|
---|
33 | return
|
---|
34 | end
|
---|
35 |
|
---|
36 | %parameters
|
---|
37 | AbsThreshold=150;
|
---|
38 | SizePart=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)
|
---|
45 | Mask=ones(size(A));
|
---|
46 | Mask(1:SizePart,:)=0;
|
---|
47 | Mask(end-SizePart:end,:)=0;
|
---|
48 | Mask(:,1:SizePart)=0;
|
---|
49 | Mask(:,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
|
---|
51 | X=zeros(size(Is));
|
---|
52 | Y=zeros(size(Js));
|
---|
53 | F=zeros(size(Js));
|
---|
54 | for 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
|
---|
72 | end
|
---|
73 | X=X(X>0);
|
---|
74 | Y=Y(Y>0);
|
---|
75 | huvmat=findobj(allchild(0),'Tag','uvmat');
|
---|
76 | if ~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')
|
---|
83 | end
|
---|
84 | % hold off
|
---|
85 | hmovie=findobj(allchild(0),'Tag','movieaxes');
|
---|
86 | if ~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')
|
---|
91 | end
|
---|
92 |
|
---|
93 |
|
---|
94 | %------------------------------------------------------------------------
|
---|
95 | % --- Find the maximum of the correlation function after interpolation
|
---|
96 | function [vector,F] = SUBPIX2DGAUSS (result_conv,x,y)
|
---|
97 | %------------------------------------------------------------------------
|
---|
98 | vector=[0 0]; %default
|
---|
99 | F=-2;
|
---|
100 | peaky=y;
|
---|
101 | peakx=x;
|
---|
102 | [npy,npx]=size(result_conv);
|
---|
103 | if (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
|
---|
133 | end
|
---|
134 | vector=[peakx-floor(npx/2)-1 peaky-floor(npy/2)-1];
|
---|