[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] | 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];
|
---|