source: trunk/src/concentration.m @ 983

Last change on this file since 983 was 977, checked in by g7moreau, 4 years ago
  • Update Copyright 2008-2017 notice
File size: 6.6 KB
Line 
1%transform LIF images to concentration images
2
3%=======================================================================
4% Copyright 2008-2017, 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,DataOut_1,DataMask]=concentration(Data,XmlData,Data_1,XmlData_1,Ref)
22cpath=which('uvmat');
23addpath(fullfile(fileparts(cpath),'transform_field'))% define path for phys_polar.m
24DataOut_1=[];
25
26%%  for use in uvmat
27num_level=Data.ZIndex;
28if ~exist('Ref','var')
29    huvmat=findobj(allchild(0),'tag','uvmat');
30    hhuvmat=guidata(huvmat);
31    RootPath=get(hhuvmat.RootPath,'String');
32   
33    %reference file
34    RootPath=fullfile(RootPath,'LIF_REF');
35    file_ref=fullfile(RootPath,['lif_ref_' num2str(num_level) '.nc']);
36    Ref=nc2struct(file_ref);
37end
38
39%% Parameters
40XmlData.GeometryCalib.PolarCentre=Ref.IlluminationOrigin;%[-515 -175]; %position of the laser origin [x, y]
41XmlData_1.GeometryCalib.PolarCentre=Ref.IlluminationOrigin;%[-515 -175]; %position of the laser origin [x, y]
42ImageOffset=Ref.ImageOffset; %237;% image value for black background
43nfilt=64;
44
45%% concentration image
46Data.A(Ref.CoverIndex:end,:)=Ref.CoverCoeff*(double(Data.A(Ref.CoverIndex:end,:))-ImageOffset(1))+ImageOffset(1);% COMPENSATION OF BRIGHTNESS UNDER THE COVER
47[DataOut,DataOut_1]=phys_polar(Data,XmlData,Data_1,XmlData_1);
48A=Ref.Aref;%default
49ind_good=find(Ref.Aref~=0);
50ind_bad=find(Ref.Aref==0);
51A(ind_good)=double(DataOut.A(ind_good))-ImageOffset(1)-0.07*(double(DataOut_1.A(ind_good))-ImageOffset(2));%substract PIV image information for removing particles
52%filtering and decimate
53Afilt=filter2(ones(nfilt,nfilt),A);
54Mask=filter2(ones(nfilt,nfilt),double(Ref.Aref~=0));
55B=Afilt./Mask;
56A(ind_bad)=B(ind_bad);
57[npy,npx]=size(A);
58DataMask=DataOut;
59DataMask.A=2*ones(npy,npx);%mask=2 for good data
60
61DataMask.A(Ref.Aref==0)=1;%mask=0 for undefined data
62
63
64
65C=filter2(ones(nfilt,nfilt),Ref.Aref);
66D=C./Mask;
67Ref.Aref(ind_bad)=D(ind_bad);
68DataOut_1=[];
69Coord_x=DataOut.Coord_x;
70Coord_y=DataOut.Coord_y;
71
72dX=(Coord_x(2)-Coord_x(1))/(npx-1);
73dY=(Coord_y(1)-Coord_y(2))/(npy-1);%mesh of new pixels
74[R,Y]=meshgrid(linspace(Coord_x(1),Coord_x(2),npx),linspace(Coord_y(1),Coord_y(2),npy));
75r=Coord_x(1)+[0:npx-1]*dX;%distance from laser
76%A(ind_good)=(A(ind_good)>=0).*A(ind_good); %replaces negative values  by zeros
77A=A./Ref.Aref;% luminosity normalised by the reference (value at the edge of the box)
78
79%% Interpolation
80% [Rindex,Yindex]=meshgrid(linspace(0.5,npx-0.5,npx),linspace(npy-0.5,0.5,npy));
81% Rgood=Rindex(ind_good);
82% Ygood=Yindex(ind_good);
83%F=TriScatteredInterp(Rgood,Ygood,A(ind_good));
84%A=F(Rindex,Yindex);
85
86
87DataMask.A(isnan(A)|isinf(A)|A>1.5)=0;% mask=1 for interpolated data
88r_edge=Ref.r_edge*ones(1,npx);
89Edge_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
90yedge=min(min(Y(Edge_ind)));
91jmax=round(-(yedge-Coord_y(1))/dY+1);
92DataMask.A(jmax:end,:)=0;
93
94A(isnan(A)|isinf(A))=0;
95
96% radius along the reference line
97Theta=(linspace(Coord_y(1),Coord_y(2),npy)*pi/180)'*ones(1,npx);%theta in radians
98
99gamma_coeff=Ref.GammaCoeff*ones(1,npx);
100
101A(R<r_edge)=0;
102I=(r_edge-dX*gamma_coeff.*cumsum(R.*A,2))./R;% expected laser intensity along the line
103
104DataOut.A=A./I;%concentration
105DataOut.A(I<=0)=0;% eliminate values obtained with I<=0
106DataOut.A(jmax:end,:)=0;%put to zeros points for which the e laser ray is not visible from the edge
107RangeX=Ref.RangeX-XmlData.GeometryCalib.PolarCentre(1);
108RangeY=Ref.RangeY-XmlData.GeometryCalib.PolarCentre(2);
109
110DataOut=polar2phys(DataOut,RangeX,RangeY);
111DataOut.Coord_x=DataOut.Coord_x+XmlData.GeometryCalib.PolarCentre(1);
112DataOut.Coord_y=DataOut.Coord_y+XmlData.GeometryCalib.PolarCentre(2);
113DataMask=polar2phys(DataMask,RangeX,RangeY);
114DataMask.Coord_x=DataMask.Coord_x+XmlData.GeometryCalib.PolarCentre(1);
115DataMask.Coord_y=DataMask.Coord_y+XmlData.GeometryCalib.PolarCentre(2);
116
117
118function DataOut=polar2phys(DataIn,RangeX,RangeY)
119%%%%%%%%%%%%%%%%%%%%
120DataOut=DataIn; %fdefault
121[npy,npx]=size(DataIn.A);
122dx=(DataIn.Coord_x(2)-DataIn.Coord_x(1))/(npx-1);
123dy=(DataIn.Coord_y(2)-DataIn.Coord_y(1))/(npy-1);%mesh
124rcorner=[DataIn.Coord_x(1) DataIn.Coord_x(2) DataIn.Coord_x(1) DataIn.Coord_x(2)];% radius of the corners
125ycorner=[DataIn.Coord_y(2) DataIn.Coord_y(2) DataIn.Coord_y(1) DataIn.Coord_y(1)];% azimuth of the corners
126thetacorner=pi*ycorner/180;% azimuth in radians
127[Xcorner,Ycorner] = pol2cart(thetacorner,rcorner);% cartesian coordinates of the corners (with respect to lser source)
128if ~exist('RangeX','var')
129RangeX(1)=min(Xcorner);
130RangeX(2)=max(Xcorner);
131end
132if ~exist('RangeY','var')
133RangeY(2)=min(Ycorner);
134RangeY(1)=max(Ycorner);
135end
136%Rangx=[-100 100];%bounds of the initial box
137%Rangy=[75 -150];
138% Rangy(1)=min(Ycorner);
139% Rangy(2)=max(Ycorner);
140x=linspace(RangeX(1),RangeX(2),npx);%coordinates of the new pixels
141y=linspace(RangeY(2),RangeY(1),npy);
142[X,Y]=meshgrid(x,y);%grid for new pixels in cartesian coordiantes
143
144[Theta,R] = cart2pol(X,Y);%corresponding polar coordiantes
145Theta=Theta*180/pi;
146%Theta=1+round((Theta-DataIn.Coord_y(1))/dy); %index along y (dy negative)
147Theta=1-round((Theta-DataIn.Coord_y(2))/dy); %index along y (dy negative)
148R=1+round((R-DataIn.Coord_x(1))/dx); %index along x
149R=reshape(R,1,npx*npy);%indices reorganized in 'line'
150Theta=reshape(Theta,1,npx*npy);
151flagin=R>=1 & R<=npx & Theta >=1 & Theta<=npy;%flagin=1 inside the original image
152vec_A=reshape(DataIn.A,1,npx*npy);%put the original image in line
153ind_in=find(flagin);
154ind_out=find(~flagin);
155ICOMB=((R-1)*npy+(npy+1-Theta));
156ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
157vec_B(ind_in)=vec_A(ICOMB);
158vec_B(ind_out)=zeros(size(ind_out));
159DataOut.A=flipdim(reshape(vec_B,npy,npx),1);%new image in real coordinates
160
161     %Rangx=Rangx-radius_ref;
162DataOut.Coord_x=RangeX;
163DataOut.Coord_y=RangeY; 
164
Note: See TracBrowser for help on using the repository browser.