source: trunk/src/transform_field/ima2concentration.m @ 1077

Last change on this file since 1077 was 1077, checked in by sommeria, 4 years ago

various bug repair

File size: 5.2 KB
Line 
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
21function [DataOut]=ima2concentration(DataIn,XmlData)
22
23%% request input parameters
24DataOut=[];
25if (isfield(DataIn,'Action') && isfield(DataIn.Action,'RUN') && isequal(DataIn.Action.RUN,0))
26    return
27end
28if ~isfield(XmlData,'LIFCalib')
29        msgbox_uvmat('ERROR','no LIF calibration data available, first run LIFCalib in uvmat')
30    return
31end
32cpath=which('uvmat');
33addpath(fullfile(fileparts(cpath),'transform_field'))% define path for phys_polar.m
34
35%% rescale the image
36[nby,nbx]=size(DataIn.A);
37x=linspace(DataIn.Coord_x(1),DataIn.Coord_x(2),nbx)-nbx/2;
38y=linspace(DataIn.Coord_y(1),DataIn.Coord_y(2),nby)-nby/2;
39[X,Y]=meshgrid(x,y);
40coeff_quad=0.15*4/(nbx*nbx);% image luminosity reduced by 10% at the edge
41DataIn.A=double(DataIn.A).*(1+coeff_quad*(X.*X+Y.*Y));
42
43%% Transform images to polar coordinates with origin at the light source position
44XmlData.TransformInput.PolarCentre=XmlData.LIFCalib.LightOrigin; %position of the laser origin [x, y]
45DataIn.Action.RUN=1;% avoid input menu in phys_polar
46DataOut=phys_polar(DataIn,XmlData);
47[npangle,npr]=size(DataOut.A);%size of the image in polar coordinates
48dX=(DataOut.Coord_x(2)-DataOut.Coord_x(1))/(npr-1);% radial step
49
50%% introduce the reference line where the laser enters the fluid region
51r_edge=XmlData.LIFCalib.RefLineRadius'*ones(1,npr);% radial position of the reference line extended as a matrix (npx,npy)
52A_ref=XmlData.LIFCalib.RefLineLum'*ones(1,npr);% luminosity on the reference line extended as a matrix (npx,npy)
53R=ones(npangle,1)*linspace(DataOut.Coord_x(1), DataOut.Coord_x(2),npr);%radial coordinate extended as a matrix (npx,npy)
54
55%gamma_coeff=XmlData.LIFCalib.DecayRate;
56DataOut.A(R<r_edge)=0;
57DataOut.A=double(DataOut.A)./A_ref;% renormalize the luminosity with the reference luminosity at the same azimuth on the reference line
58I=(r_edge-dX*XmlData.LIFCalib.DecayRate.*cumsum(R.*DataOut.A,2))./R;% expected laser intensity along the line
59DataOut.A=DataOut.A./I;%concentration normalized by the uniform concentration assumed in the ref image used for calibration
60DataOut.A(I<=0)=0;% eliminate values obtained with I<=0
61
62DataOut=polar2phys(DataOut);% back to phys cartesian coordinates with origin at the light source
63DataOut.A=uint16(1000*DataOut.A);% concentration multiplied by 1000 to get an image
64DataOut.Coord_x=DataOut.Coord_x+XmlData.LIFCalib.LightOrigin(1);%shift to original cartesian coordinates
65DataOut.Coord_y=DataOut.Coord_y+XmlData.LIFCalib.LightOrigin(2);
66
67
68function DataOut=polar2phys(DataIn)
69%%%%%%%%%%%%%%%%%%%%
70DataOut=DataIn; %default
71[npy,npx]=size(DataIn.A);
72dx=(DataIn.Coord_x(2)-DataIn.Coord_x(1))/(npx-1); %mesh along radius
73dy=(DataIn.Coord_y(2)-DataIn.Coord_y(1))/(npy-1);%mesh along azimuth
74
75%% create cartesian coordinates in the domain defined by the four image corners
76rcorner=[DataIn.Coord_x(1) DataIn.Coord_x(2) DataIn.Coord_x(1) DataIn.Coord_x(2)];% radius of the corners
77ycorner=[DataIn.Coord_y(2) DataIn.Coord_y(2) DataIn.Coord_y(1) DataIn.Coord_y(1)];% azimuth of the corners
78thetacorner=pi*ycorner/180;% azimuth in radians
79[Xcorner,Ycorner] = pol2cart(thetacorner,rcorner);% cartesian coordinates of the corners (with respect to lser source)
80RangeX(1)=min(Xcorner);
81RangeX(2)=max(Xcorner);
82RangeY(2)=min(Ycorner);
83RangeY(1)=max(Ycorner);
84x=linspace(RangeX(1),RangeX(2),npx);%coordinates of the new pixels
85y=linspace(RangeY(2),RangeY(1),npy);
86[X,Y]=meshgrid(x,y);%grid for new pixels in cartesian coordinates
87
88%% image indices corresponding to the cartesian grid
89[Theta,R] = cart2pol(X,Y);%corresponding polar coordiantes
90Theta=180*Theta/pi;%angles in degrees
91Theta=1-round((Theta-DataIn.Coord_y(2))/dy); %angular index along y (dy negative)
92R=1+round((R-DataIn.Coord_x(1))/dx); %angular index along x
93R=reshape(R,1,npx*npy);%indices reorganized in 'line'
94Theta=reshape(Theta,1,npx*npy);
95flagin=R>=1 & R<=npx & Theta >=1 & Theta<=npy;%flagin=1 inside the original image
96vec_A=reshape(DataIn.A,1,npx*npy);%put the original image in line
97ind_in=find(flagin);
98ind_out=find(~flagin);
99ICOMB=((R-1)*npy+(npy+1-Theta));
100ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
101vec_B(ind_in)=vec_A(ICOMB);
102vec_B(ind_out)=zeros(size(ind_out));
103DataOut.A=flipdim(reshape(vec_B,npy,npx),1);%new image in real coordinates
104DataOut.Coord_x=RangeX;
105DataOut.Coord_y=RangeY; 
106
Note: See TracBrowser for help on using the repository browser.