[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  %=======================================================================


[1071]  12  % Copyright 20082020, LEGI UMR 5519 / CNRS UGA GINP, Grenoble, France


[810]  13  % http://www.legi.grenobleinp.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(endSizePart:end,:)=0;


 48  Mask(:,1:SizePart)=0;


 49  Mask(:,endSizePart: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,XroundSizePart);


 65  Xhigh=min(size(A,2),Xround+SizePart);


 66  Yround=round(Y(ipart));


 67  Ylow=max(1,YroundSizePart);


 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(X0.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(X0.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 <= npx1) && (y <= npy1) && (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  %Twodimensional Gaussian regression for subpixel 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^22)*log(result_conv(y+j, x+i));


 117  c02(j+2,i+2)=(3*j^22)*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*c012*c10*c02)/(4*c20*c02c11^2);


 126  deltay=(c11*c102*c01*c20)/(4*c20*c02c11^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=[peakxfloor(npx/2)1 peakyfloor(npy/2)1];

