[566]  1  %transform LIF images to concentration images


 2 


[809]  3  %=======================================================================


 4  % Copyright 20082014, LEGI UMR 5519 / CNRS UJF GINP, Grenoble, France


 5  % http://www.legi.grenobleinp.fr


 6  % Joel.Sommeria  Joel.Sommeria (A) legi.cnrs.fr


 7  %


 8  % This file is part of the toolbox UVMAT.


 9  %


 10  % UVMAT is free software; you can redistribute it and/or modify


 11  % it under the terms of the GNU General Public License as published


 12  % by the Free Software Foundation; either version 2 of the license,


 13  % or (at your option) any later version.


 14  %


 15  % UVMAT is distributed in the hope that it will be useful,


 16  % but WITHOUT ANY WARRANTY; without even the implied warranty of


 17  % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the


 18  % GNU General Public License (see LICENSE.txt) for more details.


 19  %=======================================================================


 20 


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


 22  cpath=which('uvmat');


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


 24  DataOut_1=[];


 25 


 26  %% for use in uvmat


 27  num_level=Data.ZIndex;


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


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


 30  hhuvmat=guidata(huvmat);


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


 32 


 33  %reference file


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


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


 36  Ref=nc2struct(file_ref);


 37  end


 38 


 39  %% Parameters


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


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


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


 43  nfilt=64;


 44 


 45  %% concentration image


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


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


 48  A=Ref.Aref;%default


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


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


 51  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


 52  %filtering and decimate


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


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


 55  B=Afilt./Mask;


 56  A(ind_bad)=B(ind_bad);


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


 58  DataMask=DataOut;


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


 60 


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


 62 


 63 


 64 


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


 66  D=C./Mask;


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


 68  DataOut_1=[];


[782]  69  Coord_x=DataOut.Coord_x;


 70  Coord_y=DataOut.Coord_y;


[566]  71 


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


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


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


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


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


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


 78 


 79  %% Interpolation


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


 81  % Rgood=Rindex(ind_good);


 82  % Ygood=Yindex(ind_good);


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


 84  %A=F(Rindex,Yindex);


 85 


 86 


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


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


 89  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


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


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


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


 93 


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


 95 


 96  % radius along the reference line


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


[566]  98 


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


 100 


 101  A(R<r_edge)=0;


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


 103 


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


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


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


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


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


 109 


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


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


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


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


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


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


[566]  116 


 117 


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


 119  %%%%%%%%%%%%%%%%%%%%


 120  DataOut=DataIn; %fdefault


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


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


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


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


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


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


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


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


 129  RangeX(1)=min(Xcorner);


 130  RangeX(2)=max(Xcorner);


 131  end


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


 133  RangeY(2)=min(Ycorner);


 134  RangeY(1)=max(Ycorner);


 135  end


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


 137  %Rangy=[75 150];


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


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


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


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


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


 143 


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


 145  Theta=Theta*180/pi;


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


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


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


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


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


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


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


 153  ind_in=find(flagin);


 154  ind_out=find(~flagin);


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


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


 157  vec_B(ind_in)=vec_A(ICOMB);


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


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


 160 


 161  %Rangx=Rangxradius_ref;


[782]  162  DataOut.Coord_x=RangeX;


 163  DataOut.Coord_y=RangeY;


[566]  164 

