source: trunk/src/concentration.m @ 643

Last change on this file since 643 was 566, checked in by sommeria, 11 years ago

introduce LIF, first step

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