Ignore:
Timestamp:
Jan 18, 2020, 10:47:39 AM (5 years ago)
Author:
sommeria
Message:

LIF bug corrected

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/transform_field/ima2concentration.m

    r1073 r1074  
    2020
    2121function [DataOut]=ima2concentration(DataIn,XmlData)
     22
    2223%% request input parameters
    2324if isfield(DataIn,'Action') && isfield(DataIn.Action,'RUN') && isequal(DataIn.Action.RUN,0)
     
    2930cpath=which('uvmat');
    3031addpath(fullfile(fileparts(cpath),'transform_field'))% define path for phys_polar.m
    31 DataOut_1=[];
    3232
    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
     33%% Transform images to polar coordinates with origin at the light source position
    4834XmlData.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
    5735DataIn.Action.RUN=1;% avoid input menu in phys_polar
    5836DataOut=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)
     37[npangle,npr]=size(DataOut.A);%size of the image in polar coordinates
     38dX=(DataOut.Coord_x(2)-DataOut.Coord_x(1))/(npr-1);% radial step
    8739
    88 [npangle,npr]=size(DataOut.A);
    89 dX=(DataOut.Coord_x(2)-DataOut.Coord_x(1))/(npr-1);
     40%% introduce the reference line where the laser enters the fluid region
    9041r_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)
     42A_ref=XmlData.LIFCalib.RefLineLum'*ones(1,npr);% luminosity on the reference line extended as a matrix (npx,npy)
    9243R=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;
    9944
    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 
     45%gamma_coeff=XmlData.LIFCalib.DecayRate;
    10646DataOut.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
     47DataOut.A=double(DataOut.A)./A_ref;% renormalize the luminosity with the reference luminosity at the same azimuth on the reference line
     48I=(r_edge-dX*XmlData.LIFCalib.DecayRate.*cumsum(R.*DataOut.A,2))./R;% expected laser intensity along the line
     49DataOut.A=DataOut.A./I;%concentration normalized by the uniform concentration assumed in the ref image used for calibration
    11050DataOut.A(I<=0)=0;% eliminate values obtained with I<=0
    11151
    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);
     52DataOut=polar2phys(DataOut);% back to phys cartesian coordinates with origin at the light source
     53DataOut.A=uint16(1000*DataOut.A);% concentration multiplied by 1000 to get an image
     54DataOut.Coord_x=DataOut.Coord_x+XmlData.LIFCalib.LightOrigin(1);%shift to original cartesian coordinates
    11855DataOut.Coord_y=DataOut.Coord_y+XmlData.LIFCalib.LightOrigin(2);
    11956
    12057
     58function DataOut=polar2phys(DataIn)
     59%%%%%%%%%%%%%%%%%%%%
     60DataOut=DataIn; %default
     61[npy,npx]=size(DataIn.A);
     62dx=(DataIn.Coord_x(2)-DataIn.Coord_x(1))/(npx-1); %mesh along radius
     63dy=(DataIn.Coord_y(2)-DataIn.Coord_y(1))/(npy-1);%mesh along azimuth
    12164
    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
     65%% create cartesian coordinates in the domain defined by the four image corners
    12866rcorner=[DataIn.Coord_x(1) DataIn.Coord_x(2) DataIn.Coord_x(1) DataIn.Coord_x(2)];% radius of the corners
    12967ycorner=[DataIn.Coord_y(2) DataIn.Coord_y(2) DataIn.Coord_y(1) DataIn.Coord_y(1)];% azimuth of the corners
    13068thetacorner=pi*ycorner/180;% azimuth in radians
    13169[Xcorner,Ycorner] = pol2cart(thetacorner,rcorner);% cartesian coordinates of the corners (with respect to lser source)
    132 if ~exist('RangeX','var')
    13370RangeX(1)=min(Xcorner);
    13471RangeX(2)=max(Xcorner);
    135 end
    136 if ~exist('RangeY','var')
    13772RangeY(2)=min(Ycorner);
    13873RangeY(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);
    14474x=linspace(RangeX(1),RangeX(2),npx);%coordinates of the new pixels
    14575y=linspace(RangeY(2),RangeY(1),npy);
    146 [X,Y]=meshgrid(x,y);%grid for new pixels in cartesian coordiantes
     76[X,Y]=meshgrid(x,y);%grid for new pixels in cartesian coordinates
    14777
     78%% image indices corresponding to the cartesian grid
    14879[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
     80Theta=180*Theta/pi;%angles in degrees
     81Theta=1-round((Theta-DataIn.Coord_y(2))/dy); %angular index along y (dy negative)
     82R=1+round((R-DataIn.Coord_x(1))/dx); %angular index along x
    15383R=reshape(R,1,npx*npy);%indices reorganized in 'line'
    15484Theta=reshape(Theta,1,npx*npy);
     
    16292vec_B(ind_out)=zeros(size(ind_out));
    16393DataOut.A=flipdim(reshape(vec_B,npy,npx),1);%new image in real coordinates
    164 
    165      %Rangx=Rangx-radius_ref;
    16694DataOut.Coord_x=RangeX;
    16795DataOut.Coord_y=RangeY; 
Note: See TracChangeset for help on using the changeset viewer.