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)
|
---|
22 | %% request input parameters
|
---|
23 | if isfield(DataIn,'Action') && isfield(DataIn.Action,'RUN') && isequal(DataIn.Action.RUN,0)
|
---|
24 | if ~isfield(XmlData,'LIFCalib')
|
---|
25 | msgbox_uvmat('ERROR','no LIF calibration data available, first run LIFCalib in uvmat')
|
---|
26 | return
|
---|
27 | end
|
---|
28 | end
|
---|
29 | cpath=which('uvmat');
|
---|
30 | addpath(fullfile(fileparts(cpath),'transform_field'))% define path for phys_polar.m
|
---|
31 | DataOut_1=[];
|
---|
32 |
|
---|
33 |
|
---|
34 | %% for use in uvmat
|
---|
35 | % num_level=Data.ZIndex;
|
---|
36 | % if ~exist('Ref','var')
|
---|
37 | % huvmat=findobj(allchild(0),'tag','uvmat');
|
---|
38 | % hhuvmat=guidata(huvmat);
|
---|
39 | % RootPath=get(hhuvmat.RootPath,'String');
|
---|
40 | %
|
---|
41 | % %reference file
|
---|
42 | % RootPath=fullfile(RootPath,'LIF_REF');
|
---|
43 | % file_ref=fullfile(RootPath,['lif_ref_' num2str(num_level) '.nc']);
|
---|
44 | % Ref=nc2struct(file_ref);
|
---|
45 | % end
|
---|
46 |
|
---|
47 | %% Parameters
|
---|
48 | XmlData.TransformInput.PolarCentre=XmlData.LIFCalib.LightOrigin; %position of the laser origin [x, y]
|
---|
49 |
|
---|
50 | % if isfield(XmlData.TransformInput,'PolarReferenceRadius')
|
---|
51 | % def{2}=num2str(XmlData.TransformInput.PolarReferenceRadius);
|
---|
52 | % end
|
---|
53 | % if isfield(XmlData.TransformInput,'PolarReferenceAngle')
|
---|
54 | % def{3}=num2str(XmlData.TransformInput.PolarReferenceAngle);
|
---|
55 |
|
---|
56 | %% concentration image
|
---|
57 | DataIn.Action.RUN=1;% avoid input menu in phys_polar
|
---|
58 | DataOut=phys_polar(DataIn,XmlData);
|
---|
59 | % A=Ref.Aref;%default
|
---|
60 | % ind_good=find(Ref.Aref~=0);
|
---|
61 | % ind_bad=find(Ref.Aref==0);
|
---|
62 | % A(ind_good)=double(DataOut.A(ind_good))-ImageOffset(1)
|
---|
63 | % %filtering and decimate
|
---|
64 | % Afilt=filter2(ones(nfilt,nfilt),A);
|
---|
65 | % Mask=filter2(ones(nfilt,nfilt),double(Ref.Aref~=0));
|
---|
66 | % B=Afilt./Mask;
|
---|
67 | % A(ind_bad)=B(ind_bad);
|
---|
68 | % [npy,npx]=size(A);
|
---|
69 | % DataMask=DataOut;
|
---|
70 | % DataMask.A=2*ones(npy,npx);%mask=2 for good data
|
---|
71 | %
|
---|
72 | % DataMask.A(Ref.Aref==0)=1;%mask=0 for undefined data
|
---|
73 | %
|
---|
74 | % C=filter2(ones(nfilt,nfilt),Ref.Aref);
|
---|
75 | % D=C./Mask;
|
---|
76 | % Ref.Aref(ind_bad)=D(ind_bad);
|
---|
77 | % DataOut_1=[];
|
---|
78 | % Coord_x=DataOut.Coord_x;
|
---|
79 | % Coord_y=DataOut.Coord_y;
|
---|
80 | %
|
---|
81 | % dX=(Coord_x(2)-Coord_x(1))/(npx-1);
|
---|
82 | % dY=(Coord_y(1)-Coord_y(2))/(npy-1);%mesh of new pixels
|
---|
83 | % [R,Y]=meshgrid(linspace(Coord_x(1),Coord_x(2),npx),linspace(Coord_y(1),Coord_y(2),npy));
|
---|
84 | % r=Coord_x(1)+[0:npx-1]*dX;%distance from laser
|
---|
85 | % %A(ind_good)=(A(ind_good)>=0).*A(ind_good); %replaces negative values by zeros
|
---|
86 | % A=A./Ref.Aref;% luminosity normalised by the reference (value at the edge of the box)
|
---|
87 |
|
---|
88 | [npangle,npr]=size(DataOut.A);
|
---|
89 | dX=(DataOut.Coord_x(2)-DataOut.Coord_x(1))/(npr-1);
|
---|
90 | r_edge=XmlData.LIFCalib.RefLineRadius'*ones(1,npr);% radial position of the reference line extended as a matrix (npx,npy)
|
---|
91 | A_ref=XmlData.LIFCalib.RefLineLum'*ones(1,npr);% luminosity on the reference line at the edge of the box,extended as a matrix (npx,npy)
|
---|
92 | R=ones(npangle,1)*linspace(DataOut.Coord_x(1), DataOut.Coord_x(2),npr);%radial coordinate extended as a matrix (npx,npy)
|
---|
93 | %Edge_ind=find((abs(R-r_edge)/dX)<=1 & DataMask.A~=0);%indies of positions close to r_edge, values greater than 1 are not expected
|
---|
94 | %yedge=min(min(Y(Edge_ind)));
|
---|
95 | % jmax=round(-(yedge-Coord_y(1))/dY+1);
|
---|
96 | % DataMask.A(jmax:end,:)=0;
|
---|
97 | %
|
---|
98 | % A(isnan(A)|isinf(A))=0;
|
---|
99 |
|
---|
100 | % radius along the reference line
|
---|
101 | %Theta=(linspace(Coord_y(1),Coord_y(2),npy)*pi/180)'*ones(1,npx);%theta in radians
|
---|
102 |
|
---|
103 | gamma_coeff=XmlData.LIFCalib.DecayRate;
|
---|
104 |
|
---|
105 |
|
---|
106 | DataOut.A(R<r_edge)=0;
|
---|
107 | DataOut.A=double(DataOut.A)./A_ref;
|
---|
108 | I=(r_edge-dX*gamma_coeff.*cumsum(R.*DataOut.A,2))./R;% expected laser intensity along the line
|
---|
109 | DataOut.A=DataOut.A./I;%concentration
|
---|
110 | DataOut.A(I<=0)=0;% eliminate values obtained with I<=0
|
---|
111 |
|
---|
112 | RangeX=DataIn.Coord_x-XmlData.LIFCalib.LightOrigin(1);
|
---|
113 | RangeY=DataIn.Coord_y-XmlData.LIFCalib.LightOrigin(2);
|
---|
114 | %
|
---|
115 | DataOut=polar2phys(DataOut,RangeX,RangeY);
|
---|
116 | DataOut.A=uint16(DataOut.A);
|
---|
117 | DataOut.Coord_x=DataOut.Coord_x+XmlData.LIFCalib.LightOrigin(1);
|
---|
118 | DataOut.Coord_y=DataOut.Coord_y+XmlData.LIFCalib.LightOrigin(2);
|
---|
119 |
|
---|
120 |
|
---|
121 |
|
---|
122 | function DataOut=polar2phys(DataIn,RangeX,RangeY)
|
---|
123 | %%%%%%%%%%%%%%%%%%%%
|
---|
124 | DataOut=DataIn; %fdefault
|
---|
125 | [npy,npx]=size(DataIn.A);
|
---|
126 | dx=(DataIn.Coord_x(2)-DataIn.Coord_x(1))/(npx-1);
|
---|
127 | dy=(DataIn.Coord_y(2)-DataIn.Coord_y(1))/(npy-1);%mesh
|
---|
128 | rcorner=[DataIn.Coord_x(1) DataIn.Coord_x(2) DataIn.Coord_x(1) DataIn.Coord_x(2)];% radius of the corners
|
---|
129 | ycorner=[DataIn.Coord_y(2) DataIn.Coord_y(2) DataIn.Coord_y(1) DataIn.Coord_y(1)];% azimuth of the corners
|
---|
130 | thetacorner=pi*ycorner/180;% azimuth in radians
|
---|
131 | [Xcorner,Ycorner] = pol2cart(thetacorner,rcorner);% cartesian coordinates of the corners (with respect to lser source)
|
---|
132 | if ~exist('RangeX','var')
|
---|
133 | RangeX(1)=min(Xcorner);
|
---|
134 | RangeX(2)=max(Xcorner);
|
---|
135 | end
|
---|
136 | if ~exist('RangeY','var')
|
---|
137 | RangeY(2)=min(Ycorner);
|
---|
138 | RangeY(1)=max(Ycorner);
|
---|
139 | end
|
---|
140 | %Rangx=[-100 100];%bounds of the initial box
|
---|
141 | %Rangy=[75 -150];
|
---|
142 | % Rangy(1)=min(Ycorner);
|
---|
143 | % Rangy(2)=max(Ycorner);
|
---|
144 | x=linspace(RangeX(1),RangeX(2),npx);%coordinates of the new pixels
|
---|
145 | y=linspace(RangeY(2),RangeY(1),npy);
|
---|
146 | [X,Y]=meshgrid(x,y);%grid for new pixels in cartesian coordiantes
|
---|
147 |
|
---|
148 | [Theta,R] = cart2pol(X,Y);%corresponding polar coordiantes
|
---|
149 | Theta=Theta*180/pi;
|
---|
150 | %Theta=1+round((Theta-DataIn.Coord_y(1))/dy); %index along y (dy negative)
|
---|
151 | Theta=1-round((Theta-DataIn.Coord_y(2))/dy); %index along y (dy negative)
|
---|
152 | R=1+round((R-DataIn.Coord_x(1))/dx); %index along x
|
---|
153 | R=reshape(R,1,npx*npy);%indices reorganized in 'line'
|
---|
154 | Theta=reshape(Theta,1,npx*npy);
|
---|
155 | flagin=R>=1 & R<=npx & Theta >=1 & Theta<=npy;%flagin=1 inside the original image
|
---|
156 | vec_A=reshape(DataIn.A,1,npx*npy);%put the original image in line
|
---|
157 | ind_in=find(flagin);
|
---|
158 | ind_out=find(~flagin);
|
---|
159 | ICOMB=((R-1)*npy+(npy+1-Theta));
|
---|
160 | ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
|
---|
161 | vec_B(ind_in)=vec_A(ICOMB);
|
---|
162 | vec_B(ind_out)=zeros(size(ind_out));
|
---|
163 | DataOut.A=flipdim(reshape(vec_B,npy,npx),1);%new image in real coordinates
|
---|
164 |
|
---|
165 | %Rangx=Rangx-radius_ref;
|
---|
166 | DataOut.Coord_x=RangeX;
|
---|
167 | DataOut.Coord_y=RangeY;
|
---|
168 |
|
---|