[1073]  1  %transform LIF images to concentration images


 2 


 3  %=======================================================================


 4  % Copyright 20082019, LEGI UMR 5519 / CNRS UGA 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 


 21  function [DataOut]=ima2concentration(DataIn,XmlData)


[1074]  22 


[1073]  23  %% request input parameters


[1075]  24  DataOut=[];


 25  if (isfield(DataIn,'Action') && isfield(DataIn.Action,'RUN') && isequal(DataIn.Action.RUN,0))


 26  return


 27  end


 28  if ~isfield(XmlData,'LIFCalib')


[1073]  29  msgbox_uvmat('ERROR','no LIF calibration data available, first run LIFCalib in uvmat')


 30  return


 31  end


 32  cpath=which('uvmat');


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


 34 


[1077]  35  %% rescale the image


 36  [nby,nbx]=size(DataIn.A);


 37  x=linspace(DataIn.Coord_x(1),DataIn.Coord_x(2),nbx)nbx/2;


 38  y=linspace(DataIn.Coord_y(1),DataIn.Coord_y(2),nby)nby/2;


 39  [X,Y]=meshgrid(x,y);


 40  coeff_quad=0.15*4/(nbx*nbx);% image luminosity reduced by 10% at the edge


 41  DataIn.A=double(DataIn.A).*(1+coeff_quad*(X.*X+Y.*Y));


 42 


[1074]  43  %% Transform images to polar coordinates with origin at the light source position


[1073]  44  XmlData.TransformInput.PolarCentre=XmlData.LIFCalib.LightOrigin; %position of the laser origin [x, y]


 45  DataIn.Action.RUN=1;% avoid input menu in phys_polar


 46  DataOut=phys_polar(DataIn,XmlData);


[1074]  47  [npangle,npr]=size(DataOut.A);%size of the image in polar coordinates


 48  dX=(DataOut.Coord_x(2)DataOut.Coord_x(1))/(npr1);% radial step


[1073]  49 


[1074]  50  %% introduce the reference line where the laser enters the fluid region


[1073]  51  r_edge=XmlData.LIFCalib.RefLineRadius'*ones(1,npr);% radial position of the reference line extended as a matrix (npx,npy)


[1074]  52  A_ref=XmlData.LIFCalib.RefLineLum'*ones(1,npr);% luminosity on the reference line extended as a matrix (npx,npy)


[1073]  53  R=ones(npangle,1)*linspace(DataOut.Coord_x(1), DataOut.Coord_x(2),npr);%radial coordinate extended as a matrix (npx,npy)


 54 


[1074]  55  %gamma_coeff=XmlData.LIFCalib.DecayRate;


[1073]  56  DataOut.A(R<r_edge)=0;


[1074]  57  DataOut.A=double(DataOut.A)./A_ref;% renormalize the luminosity with the reference luminosity at the same azimuth on the reference line


 58  I=(r_edgedX*XmlData.LIFCalib.DecayRate.*cumsum(R.*DataOut.A,2))./R;% expected laser intensity along the line


 59  DataOut.A=DataOut.A./I;%concentration normalized by the uniform concentration assumed in the ref image used for calibration


[1073]  60  DataOut.A(I<=0)=0;% eliminate values obtained with I<=0


 61 


[1074]  62  DataOut=polar2phys(DataOut);% back to phys cartesian coordinates with origin at the light source


 63  DataOut.A=uint16(1000*DataOut.A);% concentration multiplied by 1000 to get an image


 64  DataOut.Coord_x=DataOut.Coord_x+XmlData.LIFCalib.LightOrigin(1);%shift to original cartesian coordinates


[1073]  65  DataOut.Coord_y=DataOut.Coord_y+XmlData.LIFCalib.LightOrigin(2);


 66 


 67 


[1074]  68  function DataOut=polar2phys(DataIn)


[1073]  69  %%%%%%%%%%%%%%%%%%%%


[1074]  70  DataOut=DataIn; %default


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


[1074]  72  dx=(DataIn.Coord_x(2)DataIn.Coord_x(1))/(npx1); %mesh along radius


 73  dy=(DataIn.Coord_y(2)DataIn.Coord_y(1))/(npy1);%mesh along azimuth


 74 


 75  %% create cartesian coordinates in the domain defined by the four image corners


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


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


 78  thetacorner=pi*ycorner/180;% azimuth in radians


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


 80  RangeX(1)=min(Xcorner);


 81  RangeX(2)=max(Xcorner);


 82  RangeY(2)=min(Ycorner);


 83  RangeY(1)=max(Ycorner);


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


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


[1074]  86  [X,Y]=meshgrid(x,y);%grid for new pixels in cartesian coordinates


[1073]  87 


[1074]  88  %% image indices corresponding to the cartesian grid


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


[1074]  90  Theta=180*Theta/pi;%angles in degrees


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


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


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


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


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


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


 97  ind_in=find(flagin);


 98  ind_out=find(~flagin);


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


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


 101  vec_B(ind_in)=vec_A(ICOMB);


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


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


 104  DataOut.Coord_x=RangeX;


 105  DataOut.Coord_y=RangeY;


 106 

