[566]  1  %transform LIF images to concentration images


 2 


 3  function [DataOut,DataOut_1,DataMask]=concentration(Data,XmlData,Data_1,XmlData_1,Ref)


 4  cpath=which('uvmat');


 5  addpath(fullfile(fileparts(cpath),'transform_field'))% define path for phys_polar.m


 6  DataOut_1=[];


 7 


 8  %% for use in uvmat


 9  num_level=Data.ZIndex;


 10  if ~exist('Ref','var')


 11  huvmat=findobj(allchild(0),'tag','uvmat');


 12  hhuvmat=guidata(huvmat);


 13  RootPath=get(hhuvmat.RootPath,'String');


 14 


 15  %reference file


 16  RootPath=fullfile(RootPath,'LIF_REF');


 17  file_ref=fullfile(RootPath,['lif_ref_' num2str(num_level) '.nc']);


 18  Ref=nc2struct(file_ref);


 19  end


 20 


 21  %% Parameters


 22  XmlData.GeometryCalib.PolarCentre=Ref.IlluminationOrigin;%[515 175]; %position of the laser origin [x, y]


 23  XmlData_1.GeometryCalib.PolarCentre=Ref.IlluminationOrigin;%[515 175]; %position of the laser origin [x, y]


 24  ImageOffset=Ref.ImageOffset; %237;% image value for black background


 25  nfilt=64;


 26 


 27  %% concentration image


 28  Data.A(Ref.CoverIndex:end,:)=Ref.CoverCoeff*(double(Data.A(Ref.CoverIndex:end,:))ImageOffset(1))+ImageOffset(1);% COMPENSATION OF BRIGHTNESS UNDER THE COVER


 29  [DataOut,DataOut_1]=phys_polar(Data,XmlData,Data_1,XmlData_1);


 30  A=Ref.Aref;%default


 31  ind_good=find(Ref.Aref~=0);


 32  ind_bad=find(Ref.Aref==0);


 33  A(ind_good)=double(DataOut.A(ind_good))ImageOffset(1)0.07*(double(DataOut_1.A(ind_good))ImageOffset(2));%substract PIV image information for removing particles


 34  %filtering and decimate


 35  Afilt=filter2(ones(nfilt,nfilt),A);


 36  Mask=filter2(ones(nfilt,nfilt),double(Ref.Aref~=0));


 37  B=Afilt./Mask;


 38  A(ind_bad)=B(ind_bad);


 39  [npy,npx]=size(A);


 40  DataMask=DataOut;


 41  DataMask.A=2*ones(npy,npx);%mask=2 for good data


 42 


 43  DataMask.A(Ref.Aref==0)=1;%mask=0 for undefined data


 44 


 45 


 46 


 47  C=filter2(ones(nfilt,nfilt),Ref.Aref);


 48  D=C./Mask;


 49  Ref.Aref(ind_bad)=D(ind_bad);


 50  DataOut_1=[];


[782]  51  Coord_x=DataOut.Coord_x;


 52  Coord_y=DataOut.Coord_y;


[566]  53 


[782]  54  dX=(Coord_x(2)Coord_x(1))/(npx1);


 55  dY=(Coord_y(1)Coord_y(2))/(npy1);%mesh of new pixels


 56  [R,Y]=meshgrid(linspace(Coord_x(1),Coord_x(2),npx),linspace(Coord_y(1),Coord_y(2),npy));


 57  r=Coord_x(1)+[0:npx1]*dX;%distance from laser


[566]  58  %A(ind_good)=(A(ind_good)>=0).*A(ind_good); %replaces negative values by zeros


 59  A=A./Ref.Aref;% luminosity normalised by the reference (value at the edge of the box)


 60 


 61  %% Interpolation


 62  % [Rindex,Yindex]=meshgrid(linspace(0.5,npx0.5,npx),linspace(npy0.5,0.5,npy));


 63  % Rgood=Rindex(ind_good);


 64  % Ygood=Yindex(ind_good);


 65  %F=TriScatteredInterp(Rgood,Ygood,A(ind_good));


 66  %A=F(Rindex,Yindex);


 67 


 68 


 69  DataMask.A(isnan(A)isinf(A)A>1.5)=0;% mask=1 for interpolated data


 70  r_edge=Ref.r_edge*ones(1,npx);


 71  Edge_ind=find((abs(Rr_edge)/dX)<=1 & DataMask.A~=0);%indies of positions close to r_edge, values greater than 1 are not expected


 72  yedge=min(min(Y(Edge_ind)));


[782]  73  jmax=round((yedgeCoord_y(1))/dY+1);


[566]  74  DataMask.A(jmax:end,:)=0;


 75 


 76  A(isnan(A)isinf(A))=0;


 77 


 78  % radius along the reference line


[782]  79  Theta=(linspace(Coord_y(1),Coord_y(2),npy)*pi/180)'*ones(1,npx);%theta in radians


[566]  80 


 81  gamma_coeff=Ref.GammaCoeff*ones(1,npx);


 82 


 83  A(R<r_edge)=0;


 84  I=(r_edgedX*gamma_coeff.*cumsum(R.*A,2))./R;% expected laser intensity along the line


 85 


 86  DataOut.A=A./I;%concentration


 87  DataOut.A(I<=0)=0;% eliminate values obtained with I<=0


 88  DataOut.A(jmax:end,:)=0;%put to zeros points for which the e laser ray is not visible from the edge


 89  RangeX=Ref.RangeXXmlData.GeometryCalib.PolarCentre(1);


 90  RangeY=Ref.RangeYXmlData.GeometryCalib.PolarCentre(2);


 91 


 92  DataOut=polar2phys(DataOut,RangeX,RangeY);


[782]  93  DataOut.Coord_x=DataOut.Coord_x+XmlData.GeometryCalib.PolarCentre(1);


 94  DataOut.Coord_y=DataOut.Coord_y+XmlData.GeometryCalib.PolarCentre(2);


[566]  95  DataMask=polar2phys(DataMask,RangeX,RangeY);


[782]  96  DataMask.Coord_x=DataMask.Coord_x+XmlData.GeometryCalib.PolarCentre(1);


 97  DataMask.Coord_y=DataMask.Coord_y+XmlData.GeometryCalib.PolarCentre(2);


[566]  98 


 99 


 100  function DataOut=polar2phys(DataIn,RangeX,RangeY)


 101  %%%%%%%%%%%%%%%%%%%%


 102  DataOut=DataIn; %fdefault


 103  [npy,npx]=size(DataIn.A);


[782]  104  dx=(DataIn.Coord_x(2)DataIn.Coord_x(1))/(npx1);


 105  dy=(DataIn.Coord_y(2)DataIn.Coord_y(1))/(npy1);%mesh


 106  rcorner=[DataIn.Coord_x(1) DataIn.Coord_x(2) DataIn.Coord_x(1) DataIn.Coord_x(2)];% radius of the corners


 107  ycorner=[DataIn.Coord_y(2) DataIn.Coord_y(2) DataIn.Coord_y(1) DataIn.Coord_y(1)];% azimuth of the corners


[566]  108  thetacorner=pi*ycorner/180;% azimuth in radians


 109  [Xcorner,Ycorner] = pol2cart(thetacorner,rcorner);% cartesian coordinates of the corners (with respect to lser source)


 110  if ~exist('RangeX','var')


 111  RangeX(1)=min(Xcorner);


 112  RangeX(2)=max(Xcorner);


 113  end


 114  if ~exist('RangeY','var')


 115  RangeY(2)=min(Ycorner);


 116  RangeY(1)=max(Ycorner);


 117  end


 118  %Rangx=[100 100];%bounds of the initial box


 119  %Rangy=[75 150];


 120  % Rangy(1)=min(Ycorner);


 121  % Rangy(2)=max(Ycorner);


 122  x=linspace(RangeX(1),RangeX(2),npx);%coordinates of the new pixels


 123  y=linspace(RangeY(2),RangeY(1),npy);


 124  [X,Y]=meshgrid(x,y);%grid for new pixels in cartesian coordiantes


 125 


 126  [Theta,R] = cart2pol(X,Y);%corresponding polar coordiantes


 127  Theta=Theta*180/pi;


[782]  128  %Theta=1+round((ThetaDataIn.Coord_y(1))/dy); %index along y (dy negative)


 129  Theta=1round((ThetaDataIn.Coord_y(2))/dy); %index along y (dy negative)


 130  R=1+round((RDataIn.Coord_x(1))/dx); %index along x


[566]  131  R=reshape(R,1,npx*npy);%indices reorganized in 'line'


 132  Theta=reshape(Theta,1,npx*npy);


 133  flagin=R>=1 & R<=npx & Theta >=1 & Theta<=npy;%flagin=1 inside the original image


 134  vec_A=reshape(DataIn.A,1,npx*npy);%put the original image in line


 135  ind_in=find(flagin);


 136  ind_out=find(~flagin);


 137  ICOMB=((R1)*npy+(npy+1Theta));


 138  ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A


 139  vec_B(ind_in)=vec_A(ICOMB);


 140  vec_B(ind_out)=zeros(size(ind_out));


 141  DataOut.A=flipdim(reshape(vec_B,npy,npx),1);%new image in real coordinates


 142 


 143  %Rangx=Rangxradius_ref;


[782]  144  DataOut.Coord_x=RangeX;


 145  DataOut.Coord_y=RangeY;


[566]  146 

