[1073] | 1 | %transform LIF images to concentration images
|
---|
| 2 |
|
---|
| 3 | %=======================================================================
|
---|
| 4 | % Copyright 2008-2019, LEGI UMR 5519 / CNRS UGA G-INP, Grenoble, France
|
---|
| 5 | % http://www.legi.grenoble-inp.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 |
|
---|
[1074] | 35 | %% Transform images to polar coordinates with origin at the light source position
|
---|
[1073] | 36 | XmlData.TransformInput.PolarCentre=XmlData.LIFCalib.LightOrigin; %position of the laser origin [x, y]
|
---|
| 37 | DataIn.Action.RUN=1;% avoid input menu in phys_polar
|
---|
| 38 | DataOut=phys_polar(DataIn,XmlData);
|
---|
[1074] | 39 | [npangle,npr]=size(DataOut.A);%size of the image in polar coordinates
|
---|
| 40 | dX=(DataOut.Coord_x(2)-DataOut.Coord_x(1))/(npr-1);% radial step
|
---|
[1073] | 41 |
|
---|
[1074] | 42 | %% introduce the reference line where the laser enters the fluid region
|
---|
[1073] | 43 | r_edge=XmlData.LIFCalib.RefLineRadius'*ones(1,npr);% radial position of the reference line extended as a matrix (npx,npy)
|
---|
[1074] | 44 | A_ref=XmlData.LIFCalib.RefLineLum'*ones(1,npr);% luminosity on the reference line extended as a matrix (npx,npy)
|
---|
[1073] | 45 | R=ones(npangle,1)*linspace(DataOut.Coord_x(1), DataOut.Coord_x(2),npr);%radial coordinate extended as a matrix (npx,npy)
|
---|
| 46 |
|
---|
[1074] | 47 | %gamma_coeff=XmlData.LIFCalib.DecayRate;
|
---|
[1073] | 48 | DataOut.A(R<r_edge)=0;
|
---|
[1074] | 49 | DataOut.A=double(DataOut.A)./A_ref;% renormalize the luminosity with the reference luminosity at the same azimuth on the reference line
|
---|
| 50 | I=(r_edge-dX*XmlData.LIFCalib.DecayRate.*cumsum(R.*DataOut.A,2))./R;% expected laser intensity along the line
|
---|
| 51 | DataOut.A=DataOut.A./I;%concentration normalized by the uniform concentration assumed in the ref image used for calibration
|
---|
[1073] | 52 | DataOut.A(I<=0)=0;% eliminate values obtained with I<=0
|
---|
| 53 |
|
---|
[1074] | 54 | DataOut=polar2phys(DataOut);% back to phys cartesian coordinates with origin at the light source
|
---|
| 55 | DataOut.A=uint16(1000*DataOut.A);% concentration multiplied by 1000 to get an image
|
---|
| 56 | DataOut.Coord_x=DataOut.Coord_x+XmlData.LIFCalib.LightOrigin(1);%shift to original cartesian coordinates
|
---|
[1073] | 57 | DataOut.Coord_y=DataOut.Coord_y+XmlData.LIFCalib.LightOrigin(2);
|
---|
| 58 |
|
---|
| 59 |
|
---|
[1074] | 60 | function DataOut=polar2phys(DataIn)
|
---|
[1073] | 61 | %%%%%%%%%%%%%%%%%%%%
|
---|
[1074] | 62 | DataOut=DataIn; %default
|
---|
[1073] | 63 | [npy,npx]=size(DataIn.A);
|
---|
[1074] | 64 | dx=(DataIn.Coord_x(2)-DataIn.Coord_x(1))/(npx-1); %mesh along radius
|
---|
| 65 | dy=(DataIn.Coord_y(2)-DataIn.Coord_y(1))/(npy-1);%mesh along azimuth
|
---|
| 66 |
|
---|
| 67 | %% create cartesian coordinates in the domain defined by the four image corners
|
---|
[1073] | 68 | rcorner=[DataIn.Coord_x(1) DataIn.Coord_x(2) DataIn.Coord_x(1) DataIn.Coord_x(2)];% radius of the corners
|
---|
| 69 | ycorner=[DataIn.Coord_y(2) DataIn.Coord_y(2) DataIn.Coord_y(1) DataIn.Coord_y(1)];% azimuth of the corners
|
---|
| 70 | thetacorner=pi*ycorner/180;% azimuth in radians
|
---|
| 71 | [Xcorner,Ycorner] = pol2cart(thetacorner,rcorner);% cartesian coordinates of the corners (with respect to lser source)
|
---|
| 72 | RangeX(1)=min(Xcorner);
|
---|
| 73 | RangeX(2)=max(Xcorner);
|
---|
| 74 | RangeY(2)=min(Ycorner);
|
---|
| 75 | RangeY(1)=max(Ycorner);
|
---|
| 76 | x=linspace(RangeX(1),RangeX(2),npx);%coordinates of the new pixels
|
---|
| 77 | y=linspace(RangeY(2),RangeY(1),npy);
|
---|
[1074] | 78 | [X,Y]=meshgrid(x,y);%grid for new pixels in cartesian coordinates
|
---|
[1073] | 79 |
|
---|
[1074] | 80 | %% image indices corresponding to the cartesian grid
|
---|
[1073] | 81 | [Theta,R] = cart2pol(X,Y);%corresponding polar coordiantes
|
---|
[1074] | 82 | Theta=180*Theta/pi;%angles in degrees
|
---|
| 83 | Theta=1-round((Theta-DataIn.Coord_y(2))/dy); %angular index along y (dy negative)
|
---|
| 84 | R=1+round((R-DataIn.Coord_x(1))/dx); %angular index along x
|
---|
[1073] | 85 | R=reshape(R,1,npx*npy);%indices reorganized in 'line'
|
---|
| 86 | Theta=reshape(Theta,1,npx*npy);
|
---|
| 87 | flagin=R>=1 & R<=npx & Theta >=1 & Theta<=npy;%flagin=1 inside the original image
|
---|
| 88 | vec_A=reshape(DataIn.A,1,npx*npy);%put the original image in line
|
---|
| 89 | ind_in=find(flagin);
|
---|
| 90 | ind_out=find(~flagin);
|
---|
| 91 | ICOMB=((R-1)*npy+(npy+1-Theta));
|
---|
| 92 | ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
|
---|
| 93 | vec_B(ind_in)=vec_A(ICOMB);
|
---|
| 94 | vec_B(ind_out)=zeros(size(ind_out));
|
---|
| 95 | DataOut.A=flipdim(reshape(vec_B,npy,npx),1);%new image in real coordinates
|
---|
| 96 | DataOut.Coord_x=RangeX;
|
---|
| 97 | DataOut.Coord_y=RangeY;
|
---|
| 98 |
|
---|