Changeset 396
 Timestamp:
 Apr 25, 2012, 2:16:56 PM (9 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/src/transform_field/phys.m
r250 r396 4 4 5 5 % OUTPUT: 6 % DataOut: structure representing the modified field7 % DataOut_1: structure representing the second modified field6 % DataOut: structure representing the first field in phys coordinates 7 % DataOut_1: structure representing the second field in phys coordinates 8 8 9 9 %INPUT: … … 14 14 % .CoordType='phys' or 'px', The function ACTS ONLY IF .CoordType='px' 15 15 % CalibData: structure containing calibration parameters or a subtree Calib.GeometryCalib =calibration data (tsai parameters) 16 % Data_1, CalibData_1: same as Data, CalibData for the second field. 16 17 17 18 function [DataOut,DataOut_1]=phys(varargin) … … 24 25 % 'D_i' if '.Role='vector_x,...', 25 26 % 'scalar', else (thenno change except scale factor) 26 Calib{1}=[]; 27 if nargin==2nargin==4 % nargin =nbre of input variables 28 Data=varargin{1}; 29 DataOut=Data;%default 30 DataOut_1=[];%default 31 CalibData=varargin{2}; 32 if isfield(CalibData,'GeometryCalib') 33 Calib{1}=CalibData.GeometryCalib; 34 end 35 Calib{2}=Calib{1}; 36 else 37 DataOut.Txt='wrong input: need two or four structures'; 38 end 39 test_1=0; 40 if nargin==4 41 test_1=1; 42 Data_1=varargin{3}; 43 DataOut_1=Data_1;%default 44 CalibData_1=varargin{4}; 45 if isfield(CalibData_1,'GeometryCalib') 46 Calib{2}=CalibData_1.GeometryCalib; 47 end 48 end 49 iscalar=0; 50 if ~isempty(Calib{1}) 51 DataOut=phys_1(Data,Calib{1}); 52 %case of images or scalar: in case of two input fields, we need to project the transform of on the same regular grid 53 if isfield(Data,'A') && isfield(Data,'AX') && ~isempty(Data.AX) && isfield(Data,'AY')&&... 54 ~isempty(Data.AY) && length(Data.A)>1 55 iscalar=1; 56 A{1}=Data.A; 57 end 58 end 59 %transform of X,Y coordinates for vector fields 60 if isfield(Data,'ZIndex')&&~isempty(Data.ZIndex)&&~isnan(Data.ZIndex) 61 ZIndex=Data.ZIndex; 27 %% analyse input and set default output 28 DataOut=varargin{1};%default first output field 29 DataOut_1=[]; %default second output field 30 if nargin>=2 % nargin =nbre of input variables 31 if isfield(varargin{2},'GeometryCalib') 32 Calib{1}=varargin{2}.GeometryCalib; 33 else 34 Calib{1}=[]; 35 end 36 if nargin>=3 %two input fields 37 DataOut_1=varargin{3};%default second output field 38 if nargin>=4 && isfield(varargin{4},'GeometryCalib') 39 Calib{2}=CalibData_1.GeometryCalib; 40 else 41 Calib{2}=Calib{1}; 42 end 43 end 44 end 45 46 %% get the z index defining the section plane 47 if isfield(varargin{1},'ZIndex')&&~isempty(varargin{1}.ZIndex)&&~isnan(varargin{1}.ZIndex) 48 ZIndex=varargin{1}.ZIndex; 62 49 else 63 50 ZIndex=1; 64 51 end 65 if test_1 66 DataOut_1=phys_1(Data_1,Calib{2}); 67 if isfield(Data_1,'A')&&isfield(Data_1,'AX')&&~isempty(Data_1.AX) && isfield(Data_1,'AY')&&... 68 ~isempty(Data_1.AY)&&length(Data_1.A)>1 69 iscalar=iscalar+1; 70 Calib{iscalar}=Calib{2}; 71 A{iscalar}=Data_1.A; 72 if isfield(Data_1,'ZIndex') && ~isequal(Data_1.ZIndex,ZIndex) 73 DataOut.Txt='inconsistent plane indexes in the two input fields'; 74 end 75 if iscalar==1% case for which only the second field is a scalar 76 [A,AX,AY]=phys_Ima(A,Calib,ZIndex); 77 DataOut_1.A=A{1}; 78 DataOut_1.AX=AX; 79 DataOut_1.AY=AY; 80 return 81 end 82 end 83 end 52 53 %% transform first field 54 iscalar=0;% counter of scalar fields 55 if ~isempty(Calib{1}) 56 if ~isfield(Calib{1},'CalibrationType')~isfield(Calib{1},'CoordUnit') 57 return %bad calib parameter input 58 end 59 if ~(isfield(varargin{1},'CoordUnit')&& strcmp(varargin{1}.CoordUnit,'pixel')) 60 return % transform only fields in pixel coordinates 61 end 62 DataOut=phys_1(varargin{1},Calib{1},ZIndex);% transform coordiantes and velocity components 63 %case of images or scalar: in case of two input fields, we need to project the transform on the same regular grid 64 if isfield(varargin{1},'A') && isfield(varargin{1},'AX') && ~isempty(varargin{1}.AX) && isfield(varargin{1},'AY')&&... 65 ~isempty(varargin{1}.AY) && length(varargin{1}.A)>1 66 iscalar=1; 67 A{1}=varargin{1}.A; 68 end 69 end 70 71 %% document the selected plane position and angle if relevant 72 if isfield(Calib{1},'SliceCoord')&&size(Calib{1}.SliceCoord,1)>=ZIndex 73 DataOut.PlaneCoord=Calib{1}.SliceCoord(ZIndex,:);% transfer the slice position corresponding to index ZIndex 74 if isfield(Calib{1},'SliceAngle') % transfer the slice rotation angles 75 if isequal(size(Calib{1}.SliceAngle,1),1)% case of a unique angle 76 DataOut.PlaneAngle=Calib{1}.SliceAngle; 77 else % case of multiple planes with different angles: select the plane with index ZIndex 78 DataOut.PlaneAngle=Calib{1}.SliceAngle(ZIndex,:); 79 end 80 end 81 end 82 83 %% transform second field if relevant 84 if ~isempty(DataOut_1) 85 if isfield(varargin{3},'ZIndex') && ~isequal(varargin{3}.ZIndex,ZIndex) 86 DataOut_1.Txt='different plane indices for the two input fields'; 87 return 88 end 89 if ~isfield(Calib{2},'CalibrationType')~isfield(Calib{2},'CoordUnit') 90 return %bad calib parameter input 91 end 92 if ~(isfield(varargin{3},'CoordUnit')&& strcmp(varargin{3}.CoordUnit,'pixel')) 93 return % transform only fields in pixel coordinates 94 end 95 DataOut_1=phys_1(Data_1,Calib{2},ZIndex); 96 if isfield(Calib{1},'SliceCoord') 97 if ~(isfield(Calib{2},'SliceCoord') && isequal(Calib{2}.SliceCoord,Calib{1}.SliceCoord)) 98 DataOut_1.Txt='different plane positions for the two input fields'; 99 return 100 end 101 DataOut_1.PlaneCoord=DataOut.PlaneCoord;% same plane position for the two input fields 102 if isfield(Calib{1},'SliceAngle') 103 if ~(isfield(Calib{2},'SliceAngle') && isequal(Calib{2}.SliceAngle,Calib{1}.SliceAngle)) 104 DataOut_1.Txt='different plane angles for the two input fields'; 105 return 106 end 107 DataOut_1.PlaneAngle=DataOut.PlaneAngle; % same plane angle for the two input fields 108 end 109 end 110 if isfield(varargin{3},'A')&&isfield(varargin{3},'AX')&&~isempty(varargin{3}.AX) && isfield(varargin{3},'AY')&&... 111 ~isempty(varargin{3}.AY)&&length(varargin{3}.A)>1 112 iscalar=iscalar+1; 113 Calib{iscalar}=Calib{2}; 114 A{iscalar}=varargin{3}.A; 115 end 116 end 117 118 %% transform the scalar(s) or image(s) 84 119 if iscalar~=0 85 120 [A,AX,AY]=phys_Ima(A,Calib,ZIndex);%TODO : introduire interp2_uvmat ds phys_ima 86 DataOut.A=A{1}; 87 DataOut.AX=AX; 88 DataOut.AY=AY; 121 if iscalar==1 && ~isempty(DataOut_1) % case for which only the second field is a scalar 122 DataOut_1.A=A{1}; 123 DataOut_1.AX=AX; 124 DataOut_1.AY=AY; 125 else 126 DataOut.A=A{1}; 127 DataOut.AX=AX; 128 DataOut.AY=AY; 129 end 89 130 if iscalar==2 90 131 DataOut_1.A=A{2}; … … 94 135 end 95 136 96 % DataOut.VarDimName{2}97 % DataOut.VarDimName{3}98 % DataOut.VarDimName{4}99 % DataOut.VarDimName{5}100 % DataOut.VarDimName{6}101 % DataOut.VarDimName{7}102 % DataOut.VarAttribute{1}103 % DataOut.VarAttribute{2}104 % DataOut.VarAttribute{3}105 % DataOut.VarAttribute{4}106 % DataOut.VarAttribute{5}107 % DataOut.VarAttribute{6}108 % DataOut.VarAttribute{7}109 137 % 110 function DataOut=phys_1(Data,Calib) 111 % for icell=1:length(Data) 112 138 % transform a single field 139 function DataOut=phys_1(Data,Calib,ZIndex) 140 % 141 %% set default output 113 142 DataOut=Data;%default 114 % DataOut.CoordUnit=Calib.CoordUnit; %put flag for physical coordinates 115 if isfield(Calib,'SliceCoord') && isfield(Data,'ZIndex')&&~isempty(Data.ZIndex)&&~isnan(Data.ZIndex) 116 DataOut.PlaneCoord=Calib.SliceCoord(Data.ZIndex,:);% transfer the slice position 117 if isfield(Calib,'SliceAngle') % transfer the slice rotation angles 118 if isequal(size(Calib.SliceAngle,1),1) 119 DataOut.PlaneAngle=Calib.SliceAngle; 120 else 121 DataOut.PlaneAngle=Calib.SliceAngle(Data.ZIndex,:); 122 end 123 end 124 end 125 % The transform ACTS ONLY IF .CoordType='px'and Calib defined 126 if isfield(Data,'CoordUnit')%&& isequal(Data.CoordType,'px')&& ~isempty(Calib) 127 if isfield(Calib,'CoordUnit') 128 DataOut.CoordUnit=Calib.CoordUnit; 129 else 130 DataOut.CoordUnit='cm'; %default 131 end 132 DataOut.TimeUnit='s'; 133 %transform of X,Y coordinates for vector fields 134 test_z=0; 135 if isfield(Data,'ZIndex') && ~isempty(Data.ZIndex)&&~isnan(Data.ZIndex) 136 Z=Data.ZIndex; 137 test_z=1; 138 else 139 Z=0; 140 end 141 if isfield(Data,'X') &&isfield(Data,'Y')&&~isempty(Data.X) && ~isempty(Data.Y) 142 [DataOut.X,DataOut.Y,DataOut.Z]=phys_XYZ(Calib,Data.X,Data.Y,Z); 143 if test_z 144 DataOut.ListVarName=[DataOut.ListVarName(1:2) {'Z'} DataOut.ListVarName(3:end)]; 145 DataOut.VarDimName=[DataOut.VarDimName(1:2) DataOut.VarDimName(1) DataOut.VarDimName(3:end)]; 146 ZAttribute{1}.Role='coord_z'; 147 DataOut.VarAttribute=[DataOut.VarAttribute(1:2) ZAttribute DataOut.VarAttribute(3:end)]; 148 end 149 if isfield(Data,'U')&&isfield(Data,'V')&&~isempty(Data.U) && ~isempty(Data.V)&& isfield(Data,'dt') 150 if ~isempty(Data.dt) 151 [XOut_1,YOut_1]=phys_XYZ(Calib,Data.XData.U/2,Data.YData.V/2,Z); 152 [XOut_2,YOut_2]=phys_XYZ(Calib,Data.X+Data.U/2,Data.Y+Data.V/2,Z); 153 DataOut.U=(XOut_2XOut_1)/Data.dt; 154 DataOut.V=(YOut_2YOut_1)/Data.dt; 155 end 156 end 157 end 158 %transform of an image or scalar: done in phys_ima 159 160 %transform of spatial derivatives 161 if isfield(Data,'X') && ~isempty(Data.X) && isfield(Data,'DjUi') && ~isempty(Data.DjUi)... 162 && isfield(Data,'dt') 163 if ~isempty(Data.dt) 164 % estimate the Jacobian matrix DXpx/DXphys 165 for ip=1:length(Data.X) 166 [Xp1,Yp1]=phys_XYZ(Calib,Data.X(ip)+0.5,Data.Y(ip),Z); 167 [Xm1,Ym1]=phys_XYZ(Calib,Data.X(ip)0.5,Data.Y(ip),Z); 168 [Xp2,Yp2]=phys_XYZ(Calib,Data.X(ip),Data.Y(ip)+0.5,Z); 169 [Xm2,Ym2]=phys_XYZ(Calib,Data.X(ip),Data.Y(ip)0.5,Z); 170 %Jacobian matrix DXpphys/DXpx 171 DjXi(1,1)=(Xp1Xm1); 172 DjXi(2,1)=(Yp1Ym1); 173 DjXi(1,2)=(Xp2Xm2); 174 DjXi(2,2)=(Yp2Ym2); 175 DjUi(:,:)=Data.DjUi(ip,:,:); 176 DjUi=(DjXi*DjUi')/DjXi;% =J1*M*J , curvature effects (derivatives of J) neglected 177 DataOut.DjUi(ip,:,:)=DjUi'; 178 end 179 DataOut.DjUi = DataOut.DjUi/Data.dt; % min(Data.DjUi(:,1,1))=DUDX 180 end 181 end 182 end 143 DataOut.CoordUnit=Calib.CoordUnit;% the output coord unit is set by the calibration parameters 144 145 %% transform X,Y coordinates for velocity fields (transform of an image or scalar done in phys_ima) 146 if isfield(Data,'X') &&isfield(Data,'Y')&&~isempty(Data.X) && ~isempty(Data.Y) 147 [DataOut.X,DataOut.Y]=phys_XYZ(Calib,Data.X,Data.Y,ZIndex); 148 Dt=1; %default 149 if isfield(Data,'dt')&&~isempty(Data.dt) 150 Dt=Data.dt; 151 end 152 if isfield(Data,'Dt')&&~isempty(Data.Dt) 153 Dt=Data.Dt; 154 end 155 if isfield(Data,'U')&&isfield(Data,'V')&&~isempty(Data.U) && ~isempty(Data.V) 156 [XOut_1,YOut_1]=phys_XYZ(Calib,Data.XData.U/2,Data.YData.V/2,ZIndex); 157 [XOut_2,YOut_2]=phys_XYZ(Calib,Data.X+Data.U/2,Data.Y+Data.V/2,ZIndex); 158 DataOut.U=(XOut_2XOut_1)/Dt; 159 DataOut.V=(YOut_2YOut_1)/Dt; 160 end 161 if ~strcmp(Calib.CalibrationType,'rescale') && isfield(Data,'X_tps') && isfield(Data,'Y_tps') 162 [DataOut.X_tps,DataOut.Y_tps]=phys_XYZ(Calib,Data.X,Data.Y,ZIndex); 163 end 164 end 165 166 %% transform of spatial derivatives: TODO check the case with plane angles 167 if isfield(Data,'X') && ~isempty(Data.X) && isfield(Data,'DjUi') && ~isempty(Data.DjUi)... 168 && isfield(Data,'dt') 169 if ~isempty(Data.dt) 170 % estimate the Jacobian matrix DXpx/DXphys 171 for ip=1:length(Data.X) 172 [Xp1,Yp1]=phys_XYZ(Calib,Data.X(ip)+0.5,Data.Y(ip),ZIndex); 173 [Xm1,Ym1]=phys_XYZ(Calib,Data.X(ip)0.5,Data.Y(ip),ZIndex); 174 [Xp2,Yp2]=phys_XYZ(Calib,Data.X(ip),Data.Y(ip)+0.5,ZIndex); 175 [Xm2,Ym2]=phys_XYZ(Calib,Data.X(ip),Data.Y(ip)0.5,ZIndex); 176 %Jacobian matrix DXpphys/DXpx 177 DjXi(1,1)=(Xp1Xm1); 178 DjXi(2,1)=(Yp1Ym1); 179 DjXi(1,2)=(Xp2Xm2); 180 DjXi(2,2)=(Yp2Ym2); 181 DjUi(:,:)=Data.DjUi(ip,:,:); 182 DjUi=(DjXi*DjUi')/DjXi;% =J1*M*J , curvature effects (derivatives of J) neglected 183 DataOut.DjUi(ip,:,:)=DjUi'; 184 end 185 DataOut.DjUi = DataOut.DjUi/Dt; % min(Data.DjUi(:,1,1))=DUDX 186 end 187 end 188 183 189 184 190 %%%%%%%%%%%%%%%%%%%% … … 217 223 for icell=1:length(A) 218 224 Calib=CalibIn{icell}; 219 if isfield(Calib,'R')  isfield(Calib,'kc') test_multi ~isequal(Calib,CalibIn{1})% the image needs to be interpolated to the new coordinates 225 % rescaling of the image coordinates without change of the image array 226 if strcmp(Calib.CalibrationType,'rescale') && isequal(Calib,CalibIn{1}) 227 A_out{icell}=A{icell};%no transform 228 Rangx=[0.5 npx0.5];%image coordiantes of corners 229 Rangy=[npy0.5 0.5]; 230 [Rangx]=phys_XYZ(Calib,Rangx,[0.5 0.5],ZIndex);%case of translations without rotation and quadratic deformation 231 [xx,Rangy]=phys_XYZ(Calib,[0.5 0.5],Rangy,ZIndex); 232 else 233 % the image needs to be interpolated to the new coordinates 220 234 zphys=0; %default 221 235 if isfield(Calib,'SliceCoord') %.Z= index of plane … … 263 277 if testuint16 264 278 A_out{icell}=uint16(A_out{icell}); 265 end 266 else% 267 A_out{icell}=A{icell};%no transform 268 Rangx=[0.5 npx0.5];%image coordiantes of corners 269 Rangy=[npy0.5 0.5]; 270 [Rangx]=phys_XYZ(Calib,Rangx,[0.5 0.5],ZIndex);%case of translations without rotation and quadratic deformation 271 [xx,Rangy]=phys_XYZ(Calib,[0.5 0.5],Rangy,ZIndex); 272 end 273 end 274 279 end 280 end 281 end 282
Note: See TracChangeset
for help on using the changeset viewer.