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

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

ima2concentration added

File size: 6.7 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%% request input parameters
23if 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
28end
29cpath=which('uvmat');
30addpath(fullfile(fileparts(cpath),'transform_field'))% define path for phys_polar.m
31DataOut_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
48XmlData.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
57DataIn.Action.RUN=1;% avoid input menu in phys_polar
58DataOut=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);
89dX=(DataOut.Coord_x(2)-DataOut.Coord_x(1))/(npr-1);
90r_edge=XmlData.LIFCalib.RefLineRadius'*ones(1,npr);% radial position of the reference line extended as a matrix (npx,npy)
91A_ref=XmlData.LIFCalib.RefLineLum'*ones(1,npr);% luminosity on the reference line at the edge of the box,extended as a matrix (npx,npy)
92R=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
103gamma_coeff=XmlData.LIFCalib.DecayRate;
104
105
106DataOut.A(R<r_edge)=0;
107DataOut.A=double(DataOut.A)./A_ref;
108I=(r_edge-dX*gamma_coeff.*cumsum(R.*DataOut.A,2))./R;% expected laser intensity along the line
109DataOut.A=DataOut.A./I;%concentration
110DataOut.A(I<=0)=0;% eliminate values obtained with I<=0
111
112RangeX=DataIn.Coord_x-XmlData.LIFCalib.LightOrigin(1);
113RangeY=DataIn.Coord_y-XmlData.LIFCalib.LightOrigin(2);
114%
115DataOut=polar2phys(DataOut,RangeX,RangeY);
116DataOut.A=uint16(DataOut.A);
117DataOut.Coord_x=DataOut.Coord_x+XmlData.LIFCalib.LightOrigin(1);
118DataOut.Coord_y=DataOut.Coord_y+XmlData.LIFCalib.LightOrigin(2);
119
120
121
122function DataOut=polar2phys(DataIn,RangeX,RangeY)
123%%%%%%%%%%%%%%%%%%%%
124DataOut=DataIn; %fdefault
125[npy,npx]=size(DataIn.A);
126dx=(DataIn.Coord_x(2)-DataIn.Coord_x(1))/(npx-1);
127dy=(DataIn.Coord_y(2)-DataIn.Coord_y(1))/(npy-1);%mesh
128rcorner=[DataIn.Coord_x(1) DataIn.Coord_x(2) DataIn.Coord_x(1) DataIn.Coord_x(2)];% radius of the corners
129ycorner=[DataIn.Coord_y(2) DataIn.Coord_y(2) DataIn.Coord_y(1) DataIn.Coord_y(1)];% azimuth of the corners
130thetacorner=pi*ycorner/180;% azimuth in radians
131[Xcorner,Ycorner] = pol2cart(thetacorner,rcorner);% cartesian coordinates of the corners (with respect to lser source)
132if ~exist('RangeX','var')
133RangeX(1)=min(Xcorner);
134RangeX(2)=max(Xcorner);
135end
136if ~exist('RangeY','var')
137RangeY(2)=min(Ycorner);
138RangeY(1)=max(Ycorner);
139end
140%Rangx=[-100 100];%bounds of the initial box
141%Rangy=[75 -150];
142% Rangy(1)=min(Ycorner);
143% Rangy(2)=max(Ycorner);
144x=linspace(RangeX(1),RangeX(2),npx);%coordinates of the new pixels
145y=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
149Theta=Theta*180/pi;
150%Theta=1+round((Theta-DataIn.Coord_y(1))/dy); %index along y (dy negative)
151Theta=1-round((Theta-DataIn.Coord_y(2))/dy); %index along y (dy negative)
152R=1+round((R-DataIn.Coord_x(1))/dx); %index along x
153R=reshape(R,1,npx*npy);%indices reorganized in 'line'
154Theta=reshape(Theta,1,npx*npy);
155flagin=R>=1 & R<=npx & Theta >=1 & Theta<=npy;%flagin=1 inside the original image
156vec_A=reshape(DataIn.A,1,npx*npy);%put the original image in line
157ind_in=find(flagin);
158ind_out=find(~flagin);
159ICOMB=((R-1)*npy+(npy+1-Theta));
160ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
161vec_B(ind_in)=vec_A(ICOMB);
162vec_B(ind_out)=zeros(size(ind_out));
163DataOut.A=flipdim(reshape(vec_B,npy,npx),1);%new image in real coordinates
164
165     %Rangx=Rangx-radius_ref;
166DataOut.Coord_x=RangeX;
167DataOut.Coord_y=RangeY; 
168
Note: See TracBrowser for help on using the repository browser.