[566] | 1 | %transform LIF images to concentration images
|
---|
| 2 |
|
---|
| 3 | function [DataOut,DataOut_1,DataMask]=concentration(Data,XmlData,Data_1,XmlData_1,Ref)
|
---|
| 4 | cpath=which('uvmat');
|
---|
| 5 | addpath(fullfile(fileparts(cpath),'transform_field'))% define path for phys_polar.m
|
---|
| 6 | DataOut_1=[];
|
---|
| 7 |
|
---|
| 8 | %% for use in uvmat
|
---|
| 9 | num_level=Data.ZIndex;
|
---|
| 10 | if ~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);
|
---|
| 19 | end
|
---|
| 20 |
|
---|
| 21 | %% Parameters
|
---|
| 22 | XmlData.GeometryCalib.PolarCentre=Ref.IlluminationOrigin;%[-515 -175]; %position of the laser origin [x, y]
|
---|
| 23 | XmlData_1.GeometryCalib.PolarCentre=Ref.IlluminationOrigin;%[-515 -175]; %position of the laser origin [x, y]
|
---|
| 24 | ImageOffset=Ref.ImageOffset; %237;% image value for black background
|
---|
| 25 | nfilt=64;
|
---|
| 26 |
|
---|
| 27 | %% concentration image
|
---|
| 28 | Data.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);
|
---|
| 30 | A=Ref.Aref;%default
|
---|
| 31 | ind_good=find(Ref.Aref~=0);
|
---|
| 32 | ind_bad=find(Ref.Aref==0);
|
---|
| 33 | A(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
|
---|
| 35 | Afilt=filter2(ones(nfilt,nfilt),A);
|
---|
| 36 | Mask=filter2(ones(nfilt,nfilt),double(Ref.Aref~=0));
|
---|
| 37 | B=Afilt./Mask;
|
---|
| 38 | A(ind_bad)=B(ind_bad);
|
---|
| 39 | [npy,npx]=size(A);
|
---|
| 40 | DataMask=DataOut;
|
---|
| 41 | DataMask.A=2*ones(npy,npx);%mask=2 for good data
|
---|
| 42 |
|
---|
| 43 | DataMask.A(Ref.Aref==0)=1;%mask=0 for undefined data
|
---|
| 44 |
|
---|
| 45 |
|
---|
| 46 |
|
---|
| 47 | C=filter2(ones(nfilt,nfilt),Ref.Aref);
|
---|
| 48 | D=C./Mask;
|
---|
| 49 | Ref.Aref(ind_bad)=D(ind_bad);
|
---|
| 50 | DataOut_1=[];
|
---|
[782] | 51 | Coord_x=DataOut.Coord_x;
|
---|
| 52 | Coord_y=DataOut.Coord_y;
|
---|
[566] | 53 |
|
---|
[782] | 54 | dX=(Coord_x(2)-Coord_x(1))/(npx-1);
|
---|
| 55 | dY=(Coord_y(1)-Coord_y(2))/(npy-1);%mesh of new pixels
|
---|
| 56 | [R,Y]=meshgrid(linspace(Coord_x(1),Coord_x(2),npx),linspace(Coord_y(1),Coord_y(2),npy));
|
---|
| 57 | r=Coord_x(1)+[0:npx-1]*dX;%distance from laser
|
---|
[566] | 58 | %A(ind_good)=(A(ind_good)>=0).*A(ind_good); %replaces negative values by zeros
|
---|
| 59 | A=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 |
|
---|
| 69 | DataMask.A(isnan(A)|isinf(A)|A>1.5)=0;% mask=1 for interpolated data
|
---|
| 70 | r_edge=Ref.r_edge*ones(1,npx);
|
---|
| 71 | 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
|
---|
| 72 | yedge=min(min(Y(Edge_ind)));
|
---|
[782] | 73 | jmax=round(-(yedge-Coord_y(1))/dY+1);
|
---|
[566] | 74 | DataMask.A(jmax:end,:)=0;
|
---|
| 75 |
|
---|
| 76 | A(isnan(A)|isinf(A))=0;
|
---|
| 77 |
|
---|
| 78 | % radius along the reference line
|
---|
[782] | 79 | Theta=(linspace(Coord_y(1),Coord_y(2),npy)*pi/180)'*ones(1,npx);%theta in radians
|
---|
[566] | 80 |
|
---|
| 81 | gamma_coeff=Ref.GammaCoeff*ones(1,npx);
|
---|
| 82 |
|
---|
| 83 | A(R<r_edge)=0;
|
---|
| 84 | I=(r_edge-dX*gamma_coeff.*cumsum(R.*A,2))./R;% expected laser intensity along the line
|
---|
| 85 |
|
---|
| 86 | DataOut.A=A./I;%concentration
|
---|
| 87 | DataOut.A(I<=0)=0;% eliminate values obtained with I<=0
|
---|
| 88 | DataOut.A(jmax:end,:)=0;%put to zeros points for which the e laser ray is not visible from the edge
|
---|
| 89 | RangeX=Ref.RangeX-XmlData.GeometryCalib.PolarCentre(1);
|
---|
| 90 | RangeY=Ref.RangeY-XmlData.GeometryCalib.PolarCentre(2);
|
---|
| 91 |
|
---|
| 92 | DataOut=polar2phys(DataOut,RangeX,RangeY);
|
---|
[782] | 93 | DataOut.Coord_x=DataOut.Coord_x+XmlData.GeometryCalib.PolarCentre(1);
|
---|
| 94 | DataOut.Coord_y=DataOut.Coord_y+XmlData.GeometryCalib.PolarCentre(2);
|
---|
[566] | 95 | DataMask=polar2phys(DataMask,RangeX,RangeY);
|
---|
[782] | 96 | DataMask.Coord_x=DataMask.Coord_x+XmlData.GeometryCalib.PolarCentre(1);
|
---|
| 97 | DataMask.Coord_y=DataMask.Coord_y+XmlData.GeometryCalib.PolarCentre(2);
|
---|
[566] | 98 |
|
---|
| 99 |
|
---|
| 100 | function DataOut=polar2phys(DataIn,RangeX,RangeY)
|
---|
| 101 | %%%%%%%%%%%%%%%%%%%%
|
---|
| 102 | DataOut=DataIn; %fdefault
|
---|
| 103 | [npy,npx]=size(DataIn.A);
|
---|
[782] | 104 | dx=(DataIn.Coord_x(2)-DataIn.Coord_x(1))/(npx-1);
|
---|
| 105 | dy=(DataIn.Coord_y(2)-DataIn.Coord_y(1))/(npy-1);%mesh
|
---|
| 106 | rcorner=[DataIn.Coord_x(1) DataIn.Coord_x(2) DataIn.Coord_x(1) DataIn.Coord_x(2)];% radius of the corners
|
---|
| 107 | ycorner=[DataIn.Coord_y(2) DataIn.Coord_y(2) DataIn.Coord_y(1) DataIn.Coord_y(1)];% azimuth of the corners
|
---|
[566] | 108 | thetacorner=pi*ycorner/180;% azimuth in radians
|
---|
| 109 | [Xcorner,Ycorner] = pol2cart(thetacorner,rcorner);% cartesian coordinates of the corners (with respect to lser source)
|
---|
| 110 | if ~exist('RangeX','var')
|
---|
| 111 | RangeX(1)=min(Xcorner);
|
---|
| 112 | RangeX(2)=max(Xcorner);
|
---|
| 113 | end
|
---|
| 114 | if ~exist('RangeY','var')
|
---|
| 115 | RangeY(2)=min(Ycorner);
|
---|
| 116 | RangeY(1)=max(Ycorner);
|
---|
| 117 | end
|
---|
| 118 | %Rangx=[-100 100];%bounds of the initial box
|
---|
| 119 | %Rangy=[75 -150];
|
---|
| 120 | % Rangy(1)=min(Ycorner);
|
---|
| 121 | % Rangy(2)=max(Ycorner);
|
---|
| 122 | x=linspace(RangeX(1),RangeX(2),npx);%coordinates of the new pixels
|
---|
| 123 | y=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
|
---|
| 127 | Theta=Theta*180/pi;
|
---|
[782] | 128 | %Theta=1+round((Theta-DataIn.Coord_y(1))/dy); %index along y (dy negative)
|
---|
| 129 | Theta=1-round((Theta-DataIn.Coord_y(2))/dy); %index along y (dy negative)
|
---|
| 130 | R=1+round((R-DataIn.Coord_x(1))/dx); %index along x
|
---|
[566] | 131 | R=reshape(R,1,npx*npy);%indices reorganized in 'line'
|
---|
| 132 | Theta=reshape(Theta,1,npx*npy);
|
---|
| 133 | flagin=R>=1 & R<=npx & Theta >=1 & Theta<=npy;%flagin=1 inside the original image
|
---|
| 134 | vec_A=reshape(DataIn.A,1,npx*npy);%put the original image in line
|
---|
| 135 | ind_in=find(flagin);
|
---|
| 136 | ind_out=find(~flagin);
|
---|
| 137 | ICOMB=((R-1)*npy+(npy+1-Theta));
|
---|
| 138 | ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
|
---|
| 139 | vec_B(ind_in)=vec_A(ICOMB);
|
---|
| 140 | vec_B(ind_out)=zeros(size(ind_out));
|
---|
| 141 | DataOut.A=flipdim(reshape(vec_B,npy,npx),1);%new image in real coordinates
|
---|
| 142 |
|
---|
| 143 | %Rangx=Rangx-radius_ref;
|
---|
[782] | 144 | DataOut.Coord_x=RangeX;
|
---|
| 145 | DataOut.Coord_y=RangeY;
|
---|
[566] | 146 |
|
---|