[854] | 1 | %'proj_field': projects the field on a projection object |
---|
| 2 | %-------------------------------------------------------------------------- |
---|
[871] | 3 | % function [ProjData,errormsg]=proj_field(FieldData,ObjectData,VarMesh) |
---|
[854] | 4 | % |
---|
| 5 | % OUTPUT: |
---|
| 6 | % ProjData structure containing the fields of the input field FieldData, |
---|
| 7 | % transmitted or projected on the object, plus the additional fields |
---|
| 8 | % .UMax, .UMin, .VMax, .VMin: min and max of velocity components in a domain |
---|
| 9 | % .UMean,VMean: mean of the velocity components in a domain |
---|
| 10 | % .AMin, AMax: min and max of a scalar |
---|
| 11 | % .AMean: mean of a scalar in a domain |
---|
| 12 | % .NbPix; |
---|
| 13 | % .DimName= names of the matrix dimensions (matlab cell) |
---|
| 14 | % .VarName= names of the variables [ProjData.VarName {'A','AMean','AMin','AMax'}]; |
---|
| 15 | % .VarDimNameIndex= dimensions of the variables, indicated by indices in the list .DimName; |
---|
| 16 | % |
---|
| 17 | %INPUT |
---|
| 18 | % ObjectData: structure characterizing the projection object |
---|
| 19 | % .Type : type of projection object |
---|
| 20 | % .ProjMode=mode of projection ; |
---|
| 21 | % .CoordUnit: 'px', 'cm' units for the coordinates defining the object |
---|
| 22 | % .Angle ( angles of rotation (=[0 0 0] by default) |
---|
| 23 | % .ProjAngle=angle of projection; |
---|
| 24 | % .DX,.DY,.DZ=increments along each coordinate |
---|
| 25 | % .Coord(nbpoints,3): set of coordinates defining the object position; |
---|
| 26 | |
---|
| 27 | %FieldData: data of the field to be projected on the projection object, with optional fields |
---|
| 28 | % .Txt: error message, transmitted to the projection |
---|
| 29 | % .FieldList: cell array of strings representing the fields to calculate |
---|
| 30 | % .CoordMesh: typical distance between data points (used for mouse action or display), transmitted |
---|
| 31 | % .CoordUnit, .TimeUnit, .dt: transmitted |
---|
| 32 | % standardised description of fields, nc-formated Matlab structure with fields: |
---|
| 33 | % .ListGlobalAttribute: cell listing the names of the global attributes |
---|
| 34 | % .Att_1,Att_2... : values of the global attributes |
---|
| 35 | % .ListVarName: cell listing the names of the variables |
---|
| 36 | % .VarAttribute: cell of structures s containing names and values of variable attributes (s.name=value) for each variable of .ListVarName |
---|
| 37 | % .Var1, .Var2....: variables (Matlab arrays) with names listed in .ListVarName |
---|
| 38 | % The variables are grouped in 'fields', made of a set of variables with common dimensions (using the function find_field_cells) |
---|
| 39 | % The variable attribute 'Role' is used to define the role for plotting: |
---|
| 40 | % Role = 'scalar': (default) represents a scalar field |
---|
| 41 | % = 'coord': represents a set of unstructured coordinates, whose |
---|
| 42 | % space dimension is given by the last array dimension (called 'NbDim'). |
---|
| 43 | % = 'coord_x', 'coord_y', 'coord_z': represents a separate set of |
---|
| 44 | % unstructured coordinate x, y or z |
---|
| 45 | % = 'vector': represents a vector field whose number of components |
---|
| 46 | % is given by the last dimension (called 'NbDim') |
---|
| 47 | % = 'vector_x', 'vector_y', 'vector_z' :represents the x, y or z component of a vector |
---|
| 48 | % = 'warnflag' : provides a warning flag about the quality of data in a 'Field', default=0, no warning |
---|
| 49 | % = 'errorflag': provides an error flag marking false data, |
---|
| 50 | % default=0, no error. Different non zero values can represent different criteria of elimination. |
---|
| 51 | % |
---|
| 52 | % Default role of variables (by name) |
---|
| 53 | % vector field: |
---|
| 54 | % .X,.Y: position of the velocity vectors, projected on the object |
---|
| 55 | % .U, .V, .W: velocity components, projected on the object |
---|
| 56 | % .C, .CName: scalar associated to the vector |
---|
| 57 | % .F : equivalent to 'warnflag' |
---|
| 58 | % .FF: equivalent to 'errorflag' |
---|
| 59 | % scalar field or image: |
---|
| 60 | % .AName: name of a scalar (to be calculated from velocity fields after projection), transmitted |
---|
| 61 | % .A: scalar, projected on the object |
---|
| 62 | % .AX, .AY: positions for the scalar |
---|
| 63 | % case of a structured grid: A is a dim 2 matrix and .AX=[first last] (length 2 vector) represents the first and last abscissa of the grid |
---|
| 64 | % case of an unstructured scalar: A is a vector, AX and AY the corresponding coordinates |
---|
| 65 | |
---|
| 66 | %======================================================================= |
---|
[924] | 67 | % Copyright 2008-2016, LEGI UMR 5519 / CNRS UGA G-INP, Grenoble, France |
---|
[854] | 68 | % http://www.legi.grenoble-inp.fr |
---|
| 69 | % Joel.Sommeria - Joel.Sommeria (A) legi.cnrs.fr |
---|
| 70 | % |
---|
| 71 | % This file is part of the toolbox UVMAT. |
---|
| 72 | % |
---|
| 73 | % UVMAT is free software; you can redistribute it and/or modify |
---|
| 74 | % it under the terms of the GNU General Public License as published |
---|
| 75 | % by the Free Software Foundation; either version 2 of the license, |
---|
| 76 | % or (at your option) any later version. |
---|
| 77 | % |
---|
| 78 | % UVMAT is distributed in the hope that it will be useful, |
---|
| 79 | % but WITHOUT ANY WARRANTY; without even the implied warranty of |
---|
| 80 | % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
---|
| 81 | % GNU General Public License (see LICENSE.txt) for more details. |
---|
| 82 | %======================================================================= |
---|
| 83 | |
---|
[871] | 84 | function [ProjData,errormsg]=proj_field(FieldData,ObjectData,VarMesh) |
---|
[854] | 85 | errormsg='';%default |
---|
| 86 | ProjData=[]; |
---|
| 87 | |
---|
| 88 | %% check input projection object: type, projection mode and Coord: |
---|
| 89 | if ~isfield(ObjectData,'Type')||~isfield(ObjectData,'ProjMode') |
---|
| 90 | return |
---|
| 91 | end |
---|
| 92 | ListProjMode={'projection','interp_lin','interp_tps','inside','outside'};%list of effective projection modes |
---|
| 93 | if isempty(find(strcmp(ObjectData.ProjMode,ListProjMode), 1))% no projection in case |
---|
| 94 | return |
---|
| 95 | end |
---|
| 96 | if ~isfield(ObjectData,'Coord')||isempty(ObjectData.Coord) |
---|
| 97 | if strcmp(ObjectData.Type,'plane') |
---|
| 98 | ObjectData.Coord=[0 0];%default |
---|
| 99 | else |
---|
| 100 | return |
---|
| 101 | end |
---|
| 102 | end |
---|
| 103 | |
---|
| 104 | %% apply projection depending on the object type |
---|
| 105 | switch ObjectData.Type |
---|
| 106 | case 'points' |
---|
| 107 | [ProjData,errormsg]=proj_points(FieldData,ObjectData); |
---|
| 108 | case {'line','polyline'} |
---|
| 109 | [ProjData,errormsg] = proj_line(FieldData,ObjectData); |
---|
| 110 | case {'polygon','rectangle','ellipse'} |
---|
| 111 | if isequal(ObjectData.ProjMode,'inside')||isequal(ObjectData.ProjMode,'outside') |
---|
[871] | 112 | if ~exist('VarMesh','var') |
---|
| 113 | VarMesh=[]; |
---|
| 114 | end |
---|
| 115 | [ProjData,errormsg] = proj_patch(FieldData,ObjectData,VarMesh); |
---|
[854] | 116 | else |
---|
| 117 | [ProjData,errormsg] = proj_line(FieldData,ObjectData); |
---|
| 118 | end |
---|
| 119 | case 'plane' |
---|
| 120 | [ProjData,errormsg] = proj_plane(FieldData,ObjectData); |
---|
| 121 | case 'volume' |
---|
| 122 | [ProjData,errormsg] = proj_volume(FieldData,ObjectData); |
---|
| 123 | end |
---|
| 124 | |
---|
| 125 | %----------------------------------------------------------------- |
---|
| 126 | %project on a set of points |
---|
| 127 | function [ProjData,errormsg]=proj_points(FieldData,ObjectData)%% |
---|
| 128 | %------------------------------------------------------------------- |
---|
| 129 | |
---|
| 130 | siz=size(ObjectData.Coord); |
---|
| 131 | width=0; |
---|
| 132 | if isfield(ObjectData,'Range') |
---|
| 133 | width=ObjectData.Range(1,2); |
---|
| 134 | end |
---|
| 135 | if isfield(ObjectData,'RangeX')&&~isempty(ObjectData.RangeX) |
---|
| 136 | width=max(ObjectData.RangeX); |
---|
| 137 | end |
---|
| 138 | if isfield(ObjectData,'RangeY')&&~isempty(ObjectData.RangeY) |
---|
| 139 | width=max(width,max(ObjectData.RangeY)); |
---|
| 140 | end |
---|
| 141 | if isfield(ObjectData,'RangeZ')&&~isempty(ObjectData.RangeZ) |
---|
| 142 | width=max(width,max(ObjectData.RangeZ)); |
---|
| 143 | end |
---|
| 144 | if isequal(ObjectData.ProjMode,'projection') |
---|
| 145 | if width==0 |
---|
| 146 | errormsg='projection range around points needed'; |
---|
| 147 | return |
---|
| 148 | end |
---|
| 149 | elseif ~isequal(ObjectData.ProjMode,'interp_lin') |
---|
| 150 | errormsg=(['ProjMode option ' ObjectData.ProjMode ' not available in proj_field']); |
---|
| 151 | return |
---|
| 152 | end |
---|
| 153 | [ProjData,errormsg]=proj_heading(FieldData,ObjectData); |
---|
| 154 | if ~isempty(errormsg) |
---|
| 155 | return |
---|
| 156 | end |
---|
| 157 | ProjData.NbDim=0; |
---|
| 158 | [CellInfo,NbDimArray,errormsg]=find_field_cells(FieldData); |
---|
| 159 | if ~isempty(errormsg) |
---|
| 160 | errormsg=['error in proj_field/proj_points:' errormsg]; |
---|
| 161 | return |
---|
| 162 | end |
---|
| 163 | %LOOP ON GROUPS OF VARIABLES SHARING THE SAME DIMENSIONS |
---|
| 164 | for icell=1:length(CellInfo) |
---|
| 165 | if NbDimArray(icell)<=1 |
---|
| 166 | continue %projection only for multidimensional fields |
---|
| 167 | end |
---|
| 168 | VarIndex=CellInfo{icell}.VarIndex;% indices of the selected variables in the list FieldData.ListVarName |
---|
| 169 | ivar_X=CellInfo{icell}.CoordIndex(end); |
---|
| 170 | ivar_Y=CellInfo{icell}.CoordIndex(end-1); |
---|
| 171 | ivar_Z=[]; |
---|
| 172 | if NbDimArray(icell)==3 |
---|
| 173 | ivar_Z=CellInfo{icell}.CoordIndex(1); |
---|
| 174 | end |
---|
| 175 | ivar_FF=[]; |
---|
| 176 | if isfield(CellInfo{icell},'VarIndex_errorflag') |
---|
| 177 | ivar_FF=CellInfo{icell}.VarIndex_errorflag; |
---|
| 178 | if numel(ivar_FF)>1 |
---|
| 179 | errormsg='multiple error flag input'; |
---|
| 180 | return |
---|
| 181 | end |
---|
| 182 | end |
---|
| 183 | % select types of variables to be projected |
---|
| 184 | ListProj={'VarIndex_scalar','VarIndex_image','VarIndex_color','VarIndex_vector_x','VarIndex_vector_y'}; |
---|
| 185 | check_proj=false(size(FieldData.ListVarName)); |
---|
| 186 | for ilist=1:numel(ListProj) |
---|
| 187 | if isfield(CellInfo{icell},ListProj{ilist}) |
---|
| 188 | check_proj(CellInfo{icell}.(ListProj{ilist}))=1; |
---|
| 189 | end |
---|
| 190 | end |
---|
| 191 | VarIndex=find(check_proj); |
---|
| 192 | ProjData.ListVarName={'Y','X','NbVal'}; |
---|
| 193 | ProjData.VarDimName={'nb_points','nb_points','nb_points'}; |
---|
| 194 | ProjData.VarAttribute{1}.Role='ancillary'; |
---|
| 195 | ProjData.VarAttribute{2}.Role='ancillary'; |
---|
| 196 | ProjData.VarAttribute{3}.Role='ancillary'; |
---|
| 197 | for ivar=VarIndex |
---|
| 198 | VarName=FieldData.ListVarName{ivar}; |
---|
| 199 | ProjData.ListVarName=[ProjData.ListVarName {VarName}];% add the current variable to the list of projected variables |
---|
| 200 | ProjData.VarDimName=[ProjData.VarDimName {'nb_points'}]; % projected VarName has a single dimension called 'nb_points' (set of projection points) |
---|
| 201 | |
---|
| 202 | end |
---|
| 203 | if strcmp( CellInfo{icell}.CoordType,'scattered') |
---|
| 204 | coord_x=FieldData.(FieldData.ListVarName{ivar_X}); |
---|
| 205 | coord_y=FieldData.(FieldData.ListVarName{ivar_Y}); |
---|
| 206 | test3D=0;% TEST 3D CASE : NOT COMPLETED , 3D CASE : NOT COMPLETED |
---|
| 207 | if length(ivar_Z)==1 |
---|
| 208 | coord_z=FieldData.(FieldData.ListVarName{ivar_Z}); |
---|
| 209 | test3D=1; |
---|
| 210 | end |
---|
| 211 | |
---|
| 212 | for ipoint=1:siz(1) |
---|
| 213 | Xpoint=ObjectData.Coord(ipoint,:); |
---|
| 214 | distX=coord_x-Xpoint(1); |
---|
| 215 | distY=coord_y-Xpoint(2); |
---|
| 216 | dist=distX.*distX+distY.*distY; |
---|
| 217 | indsel=find(dist<width*width); |
---|
| 218 | ProjData.X(ipoint,1)=Xpoint(1); |
---|
| 219 | ProjData.Y(ipoint,1)=Xpoint(2); |
---|
| 220 | if isequal(length(ivar_FF),1) |
---|
| 221 | FFName=FieldData.ListVarName{ivar_FF}; |
---|
| 222 | FF=FieldData.(FFName)(indsel); |
---|
| 223 | indsel=indsel(~FF); |
---|
| 224 | end |
---|
| 225 | ProjData.NbVal(ipoint,1)=length(indsel); |
---|
| 226 | for ivar=VarIndex |
---|
| 227 | VarName=FieldData.ListVarName{ivar}; |
---|
| 228 | if isempty(indsel) |
---|
| 229 | ProjData.(VarName)(ipoint,1)=NaN; |
---|
| 230 | else |
---|
| 231 | Var=FieldData.(VarName)(indsel); |
---|
| 232 | ProjData.(VarName)(ipoint,1)=mean(Var); |
---|
| 233 | if isequal(ObjectData.ProjMode,'interp_lin') |
---|
| 234 | ProjData.(VarName)(ipoint,1)=griddata_uvmat(coord_x(indsel),coord_y(indsel),Var,Xpoint(1),Xpoint(2)); |
---|
| 235 | end |
---|
| 236 | end |
---|
| 237 | end |
---|
| 238 | end |
---|
| 239 | else %case of structured coordinates |
---|
| 240 | if strcmp( CellInfo{icell}.CoordType,'grid') |
---|
| 241 | AYName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(end-1)}; |
---|
| 242 | AXName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)}; |
---|
| 243 | eval(['AX=FieldData.' AXName ';']);% set of x positions |
---|
| 244 | eval(['AY=FieldData.' AYName ';']);% set of y positions |
---|
| 245 | AName=FieldData.ListVarName{VarIndex(1)};% a single variable assumed in the current cell |
---|
| 246 | eval(['A=FieldData.' AName ';']);% scalar |
---|
| 247 | npxy=size(A); |
---|
| 248 | %update VarDimName in case of components (non coordinate dimensions e;g. color components) |
---|
| 249 | if numel(npxy)>NbDimArray(icell) |
---|
| 250 | ProjData.VarDimName{end}={'nb_points','component'}; |
---|
| 251 | end |
---|
| 252 | for idim=1:NbDimArray(icell) %loop on space dimensions |
---|
| 253 | test_interp(idim)=0;%test for coordiate interpolation (non regular grid), =0 by default |
---|
| 254 | test_coord(idim)=0;%test for defined coordinates, =0 by default |
---|
| 255 | ivar=CellInfo{icell}.CoordIndex(idim); |
---|
| 256 | Coord{idim}=FieldData.(FieldData.ListVarName{ivar}); % position for the first index |
---|
| 257 | if numel(Coord{idim})==2 |
---|
| 258 | DCoord_min(idim)= (Coord{idim}(2)-Coord{idim}(1))/(npxy(idim)-1); |
---|
| 259 | test_direct(idim)=DCoord_min(idim)>0;% =1 for increasing values, 0 otherwise |
---|
| 260 | else |
---|
| 261 | DCoord=diff(Coord{idim}); |
---|
| 262 | DCoord_min(idim)=min(DCoord); |
---|
| 263 | DCoord_max=max(DCoord); |
---|
| 264 | test_direct(idim)=DCoord_max>0;% =1 for increasing values, 0 otherwise |
---|
| 265 | test_direct_min=DCoord_min(idim)>0;% =1 for increasing values, 0 otherwise |
---|
| 266 | if ~isequal(test_direct(idim),test_direct_min) |
---|
| 267 | errormsg=['non monotonic dimension variable # ' num2str(idim) ' in proj_field.m']; |
---|
| 268 | return |
---|
| 269 | end |
---|
| 270 | test_interp(idim)=(DCoord_max-DCoord_min(idim))> 0.0001*abs(DCoord_max);% test grid regularity |
---|
| 271 | test_coord(idim)=1; |
---|
| 272 | end |
---|
| 273 | end |
---|
| 274 | DX=DCoord_min(2); |
---|
| 275 | DY=DCoord_min(1); |
---|
| 276 | for ipoint=1:siz(1) |
---|
| 277 | xwidth=width/(abs(DX)); |
---|
| 278 | ywidth=width/(abs(DY)); |
---|
| 279 | i_min=round((ObjectData.Coord(ipoint,1)-Coord{2}(1))/DX+0.5-xwidth); %minimum index of the selected region |
---|
| 280 | i_min=max(1,i_min);%restrict to field limit |
---|
| 281 | i_plus=round((ObjectData.Coord(ipoint,1)-Coord{2}(1))/DX+0.5+xwidth); |
---|
| 282 | i_plus=min(npxy(2),i_plus); %restrict to field limit |
---|
| 283 | j_min=round((ObjectData.Coord(ipoint,2)-Coord{1}(1))/DY-ywidth+0.5); |
---|
| 284 | j_min=max(1,j_min); |
---|
| 285 | j_plus=round((ObjectData.Coord(ipoint,2)-Coord{1}(1))/DY+ywidth+0.5); |
---|
| 286 | j_plus=min(npxy(1),j_plus); |
---|
| 287 | ProjData.X(ipoint,1)=ObjectData.Coord(ipoint,1); |
---|
| 288 | ProjData.Y(ipoint,1)=ObjectData.Coord(ipoint,2); |
---|
| 289 | i_int=(i_min:i_plus); |
---|
| 290 | j_int=(j_min:j_plus); |
---|
| 291 | ProjData.NbVal(ipoint,1)=length(j_int)*length(i_int); |
---|
| 292 | if isempty(i_int) || isempty(j_int) |
---|
| 293 | for ivar=VarIndex |
---|
| 294 | eval(['ProjData.' FieldData.ListVarName{ivar} '(ipoint,:)=NaN;']); |
---|
| 295 | end |
---|
| 296 | errormsg=['no data points in the selected projection range ' num2str(width) ]; |
---|
| 297 | else |
---|
| 298 | %TODO: introduce circle in the selected subregion |
---|
| 299 | %[I,J]=meshgrid([1:j_int],[1:i_int]); |
---|
| 300 | for ivar=VarIndex |
---|
| 301 | Avalue=FieldData.(FieldData.ListVarName{ivar})(j_int,i_int,:); |
---|
| 302 | ProjData.(FieldData.ListVarName{ivar})(ipoint,:)=mean(mean(Avalue)); |
---|
| 303 | end |
---|
| 304 | end |
---|
| 305 | end |
---|
| 306 | end |
---|
| 307 | end |
---|
| 308 | end |
---|
| 309 | |
---|
| 310 | %----------------------------------------------------------------- |
---|
| 311 | %project in a patch |
---|
[871] | 312 | function [ProjData,errormsg]=proj_patch(FieldData,ObjectData,VarMesh)%% |
---|
[854] | 313 | %------------------------------------------------------------------- |
---|
| 314 | [ProjData,errormsg]=proj_heading(FieldData,ObjectData); |
---|
| 315 | if ~isempty(errormsg) |
---|
| 316 | return |
---|
| 317 | end |
---|
| 318 | %objectfield=fieldnames(ObjectData); |
---|
| 319 | widthx=0; |
---|
| 320 | widthy=0; |
---|
| 321 | if isfield(ObjectData,'RangeX') && ~isempty(ObjectData.RangeX) |
---|
| 322 | widthx=max(ObjectData.RangeX); |
---|
| 323 | end |
---|
| 324 | if isfield(ObjectData,'RangeY') && ~isempty(ObjectData.RangeY) |
---|
| 325 | widthy=max(ObjectData.RangeY); |
---|
| 326 | end |
---|
| 327 | |
---|
| 328 | %A REVOIR, GENERALISER: UTILISER proj_line |
---|
| 329 | ProjData.NbDim=1; |
---|
| 330 | ProjData.ListVarName={}; |
---|
| 331 | ProjData.VarDimName={}; |
---|
| 332 | ProjData.VarAttribute={}; |
---|
| 333 | |
---|
| 334 | CoordMesh=zeros(1,numel(FieldData.ListVarName)); |
---|
| 335 | if isfield (FieldData,'VarAttribute') |
---|
| 336 | for iattr=1:length(FieldData.VarAttribute)%initialization of variable attribute values |
---|
| 337 | if isfield(FieldData.VarAttribute{iattr},'Unit') |
---|
| 338 | unit{iattr}=FieldData.VarAttribute{iattr}.Unit; |
---|
| 339 | end |
---|
| 340 | if isfield(FieldData.VarAttribute{iattr},'CoordMesh') |
---|
| 341 | CoordMesh(iattr)=FieldData.VarAttribute{iattr}.CoordMesh; |
---|
| 342 | end |
---|
| 343 | end |
---|
| 344 | end |
---|
| 345 | |
---|
| 346 | %group the variables (fields of 'FieldData') in cells of variables with the same dimensions |
---|
| 347 | [CellInfo,NbDim,errormsg]=find_field_cells(FieldData); |
---|
| 348 | if ~isempty(errormsg) |
---|
| 349 | errormsg=['error in proj_field/proj_patch:' errormsg]; |
---|
| 350 | return |
---|
| 351 | end |
---|
| 352 | |
---|
| 353 | %LOOP ON GROUPS OF VARIABLES SHARING THE SAME DIMENSIONS |
---|
| 354 | for icell=1:length(CellInfo) |
---|
| 355 | CoordType=CellInfo{icell}.CoordType; |
---|
| 356 | test_Amat=0; |
---|
| 357 | if NbDim(icell)~=2% proj_patch acts only on fields of space dimension 2 |
---|
| 358 | continue |
---|
| 359 | end |
---|
| 360 | ivar_FF=[]; |
---|
| 361 | testfalse=isfield(CellInfo{icell},'VarIndex_errorflag'); |
---|
| 362 | if testfalse |
---|
| 363 | ivar_FF=CellInfo{icell}.VarIndex_errorflag; |
---|
| 364 | FFName=FieldData.ListVarName{ivar_FF}; |
---|
| 365 | errorflag=FieldData.(FFName); |
---|
| 366 | end |
---|
| 367 | % select types of variables to be projected |
---|
| 368 | ListProj={'VarIndex_scalar','VarIndex_image','VarIndex_color','VarIndex_vector_x','VarIndex_vector_y'}; |
---|
| 369 | check_proj=false(size(FieldData.ListVarName)); |
---|
| 370 | for ilist=1:numel(ListProj) |
---|
| 371 | if isfield(CellInfo{icell},ListProj{ilist}) |
---|
| 372 | check_proj(CellInfo{icell}.(ListProj{ilist}))=1; |
---|
| 373 | end |
---|
| 374 | end |
---|
| 375 | VarIndex=find(check_proj); |
---|
| 376 | |
---|
| 377 | ivar_X=CellInfo{icell}.CoordIndex(end); |
---|
| 378 | ivar_Y=CellInfo{icell}.CoordIndex(end-1); |
---|
| 379 | ivar_Z=[]; |
---|
| 380 | if NbDim(icell)==3 |
---|
| 381 | ivar_Z=CellInfo{icell}.CoordIndex(1); |
---|
| 382 | end |
---|
| 383 | switch CellInfo{icell}.CoordType |
---|
| 384 | case 'scattered' %case of unstructured coordinates |
---|
| 385 | for ivar=[VarIndex ivar_X ivar_Y ivar_FF] |
---|
| 386 | VarName=FieldData.ListVarName{ivar}; |
---|
| 387 | FieldData.(VarName)=reshape(FieldData.(VarName),[],1); |
---|
| 388 | end |
---|
| 389 | XName=FieldData.ListVarName{ivar_X}; |
---|
| 390 | YName=FieldData.ListVarName{ivar_Y}; |
---|
[866] | 391 | coord_x=FieldData.(FieldData.ListVarName{ivar_X}); |
---|
| 392 | coord_y=FieldData.(FieldData.ListVarName{ivar_Y}); |
---|
[854] | 393 | % image or 2D matrix |
---|
| 394 | case 'grid' %case of structured coordinates |
---|
| 395 | test_Amat=1;% test for image or 2D matrix |
---|
| 396 | AYName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(end-1)}; |
---|
| 397 | AXName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)}; |
---|
| 398 | AX=FieldData.(AXName);% x coordinate |
---|
| 399 | AY=FieldData.(AYName);% y coordinate |
---|
| 400 | VarName=FieldData.ListVarName{VarIndex(1)}; |
---|
| 401 | DimValue=size(FieldData.(VarName)); |
---|
| 402 | if length(AX)==2 |
---|
| 403 | AX=linspace(AX(1),AX(end),DimValue(2)); |
---|
| 404 | end |
---|
| 405 | if length(AY)==2 |
---|
| 406 | AY=linspace(AY(1),AY(end),DimValue(1)); |
---|
| 407 | end |
---|
| 408 | if length(DimValue)==3 |
---|
| 409 | testcolor=1; |
---|
| 410 | npxy(3)=3; |
---|
| 411 | else |
---|
| 412 | testcolor=0; |
---|
| 413 | npxy(3)=1; |
---|
| 414 | end |
---|
| 415 | [Xi,Yi]=meshgrid(AX,AY); |
---|
| 416 | npxy(1)=length(AY); |
---|
| 417 | npxy(2)=length(AX); |
---|
| 418 | Xi=reshape(Xi,npxy(1)*npxy(2),1); |
---|
| 419 | Yi=reshape(Yi,npxy(1)*npxy(2),1); |
---|
| 420 | for ivar=1:length(VarIndex) |
---|
| 421 | VarName=FieldData.ListVarName{VarIndex(ivar)}; |
---|
| 422 | FieldData.(VarName)=reshape(FieldData.(VarName),npxy(1)*npxy(2),npxy(3)); % keep only non false vectors |
---|
| 423 | end |
---|
| 424 | end |
---|
| 425 | %select the indices in the range of action |
---|
| 426 | testin=[];%default |
---|
| 427 | if isequal(ObjectData.Type,'rectangle') |
---|
| 428 | if strcmp(CellInfo{icell}.CoordType,'scattered') |
---|
| 429 | distX=abs(coord_x-ObjectData.Coord(1,1)); |
---|
| 430 | distY=abs(coord_y-ObjectData.Coord(1,2)); |
---|
| 431 | testin=distX<widthx & distY<widthy; |
---|
| 432 | elseif test_Amat |
---|
| 433 | distX=abs(Xi-ObjectData.Coord(1,1)); |
---|
| 434 | distY=abs(Yi-ObjectData.Coord(1,2)); |
---|
| 435 | testin=distX<widthx & distY<widthy; |
---|
| 436 | end |
---|
| 437 | elseif isequal(ObjectData.Type,'polygon') |
---|
| 438 | if strcmp(CoordType,'scattered') |
---|
| 439 | testin=inpolygon(coord_x,coord_y,ObjectData.Coord(:,1),ObjectData.Coord(:,2)); |
---|
| 440 | elseif strcmp(CoordType,'grid') |
---|
| 441 | testin=inpolygon(Xi,Yi,ObjectData.Coord(:,1),ObjectData.Coord(:,2)); |
---|
| 442 | else%calculate the scalar |
---|
| 443 | testin=[]; %A REVOIR |
---|
| 444 | end |
---|
| 445 | elseif isequal(ObjectData.Type,'ellipse') |
---|
| 446 | X2Max=widthx*widthx; |
---|
| 447 | Y2Max=(widthy)*(widthy); |
---|
| 448 | if strcmp(CoordType,'scattered') |
---|
| 449 | distX=(coord_x-ObjectData.Coord(1,1)); |
---|
| 450 | distY=(coord_y-ObjectData.Coord(1,2)); |
---|
| 451 | testin=(distX.*distX/X2Max+distY.*distY/Y2Max)<1; |
---|
| 452 | elseif strcmp(CoordType,'grid') %case of usual 2x2 matrix |
---|
| 453 | distX=(Xi-ObjectData.Coord(1,1)); |
---|
| 454 | distY=(Yi-ObjectData.Coord(1,2)); |
---|
| 455 | testin=(distX.*distX/X2Max+distY.*distY/Y2Max)<1; |
---|
| 456 | end |
---|
| 457 | end |
---|
| 458 | %selected indices |
---|
| 459 | if isequal(ObjectData.ProjMode,'outside') |
---|
| 460 | testin=~testin; |
---|
| 461 | end |
---|
| 462 | if testfalse |
---|
| 463 | testin=testin & (errorflag==0); % keep only non false vectors |
---|
| 464 | end |
---|
| 465 | indsel=find(testin); |
---|
[871] | 466 | nbvar=0; |
---|
| 467 | VarSize=zeros(size(VarIndex)); |
---|
[854] | 468 | for ivar=VarIndex |
---|
| 469 | VarName=FieldData.ListVarName{ivar}; |
---|
[871] | 470 | ProjData.([VarName 'Mean'])=mean(double(FieldData.(VarName)(indsel,:))); % take the mean in the selected region, for each co |
---|
[854] | 471 | ProjData.([VarName 'Min'])=min(double(FieldData.(VarName)(indsel,:))); % take the min in the selected region , for each color component |
---|
[871] | 472 | ProjData.([VarName 'Max'])=max(double(FieldData.(VarName)(indsel,:))); % take the max in the selected region , for each color co |
---|
| 473 | nbvar=nbvar+1; |
---|
| 474 | VarSize(nbvar)=mean((ProjData.([VarName 'Max'])-ProjData.([VarName 'Min']))/100); |
---|
| 475 | end |
---|
[880] | 476 | if isempty(VarMesh)% || isnan(VarMesh) % mesh not specified as input, estimate from the bounds |
---|
[871] | 477 | VarMesh=mean(VarSize); |
---|
| 478 | ord=10^(floor(log10(VarMesh)));%order of magnitude |
---|
| 479 | if VarMesh/ord >=5 |
---|
| 480 | VarMesh=5*ord; |
---|
| 481 | elseif VarMesh/ord >=2 |
---|
| 482 | VarMesh=2*ord; |
---|
[854] | 483 | else |
---|
[871] | 484 | VarMesh=ord; |
---|
[854] | 485 | end |
---|
[871] | 486 | end |
---|
| 487 | for ivar=VarIndex |
---|
| 488 | VarName=FieldData.ListVarName{ivar}; |
---|
| 489 | LowBound=VarMesh*ceil(ProjData.([VarName 'Min'])/VarMesh); |
---|
| 490 | UpperBound=VarMesh*floor(ProjData.([VarName 'Max'])/VarMesh); |
---|
[880] | 491 | if numel(indsel)<=1 |
---|
| 492 | errormsg='only one data point or less for histogram'; |
---|
| 493 | return |
---|
| 494 | elseif isequal(LowBound,UpperBound) |
---|
| 495 | errormsg='attempt histogram of uniform field: low bound = high bound'; |
---|
| 496 | return |
---|
| 497 | end |
---|
[871] | 498 | ProjData.(VarName)=LowBound:VarMesh:UpperBound; % list of bin values |
---|
| 499 | ProjData.([VarName 'Histo'])=hist(double(FieldData.(VarName)(indsel,:)),ProjData.(VarName)); % histogram at predefined bin positions |
---|
[854] | 500 | ProjData.ListVarName=[ProjData.ListVarName {VarName} {[VarName 'Histo']} {[VarName 'Mean']} {[VarName 'Min']} {[VarName 'Max']}]; |
---|
| 501 | if test_Amat && testcolor |
---|
| 502 | ProjData.VarDimName=[ProjData.VarDimName {VarName} {{VarName,'rgb'}} {'rgb'} {'rgb'} {'rgb'}];%{{'nb_point','rgb'}}; |
---|
| 503 | else |
---|
| 504 | ProjData.VarDimName=[ProjData.VarDimName {VarName} {VarName} {'one'} {'one'} {'one'}]; |
---|
| 505 | end |
---|
[871] | 506 | VarAttribute_var=[]; |
---|
[854] | 507 | if isfield(FieldData,'VarAttribute')&& numel(FieldData.VarAttribute)>=ivar |
---|
[871] | 508 | VarAttribute_var=FieldData.VarAttribute{ivar}; |
---|
[854] | 509 | end |
---|
[871] | 510 | % VarAttribute_var.Role='coord_x';% the variable is now used as an absissa |
---|
| 511 | VarAttribute_histo.Role='histo'; |
---|
| 512 | ProjData.VarAttribute=[ProjData.VarAttribute {VarAttribute_var} {VarAttribute_histo} {[]} {[]} {[]}]; |
---|
[854] | 513 | end |
---|
| 514 | end |
---|
| 515 | |
---|
| 516 | %----------------------------------------------------------------- |
---|
| 517 | %project on a line |
---|
| 518 | % AJOUTER flux,circul,error |
---|
| 519 | % OUTPUT: |
---|
| 520 | % ProjData: projected field |
---|
| 521 | % |
---|
| 522 | function [ProjData,errormsg] = proj_line(FieldData, ObjectData) |
---|
| 523 | %----------------------------------------------------------------- |
---|
[905] | 524 | |
---|
| 525 | %% prepare heading for the projected field |
---|
[854] | 526 | [ProjData,errormsg]=proj_heading(FieldData,ObjectData);%transfer global attributes |
---|
| 527 | if ~isempty(errormsg) |
---|
| 528 | return |
---|
| 529 | end |
---|
| 530 | ProjData.NbDim=1; |
---|
| 531 | %initialisation of the input parameters and defaultoutput |
---|
| 532 | ProjMode=ObjectData.ProjMode; %rmq: ProjMode always defined from input={'projection','interp_lin','interp_tps'} |
---|
| 533 | % ProjAngle=90; %90 degrees projection by default |
---|
| 534 | width=0; |
---|
| 535 | if isfield(ObjectData,'RangeY') |
---|
| 536 | width=max(ObjectData.RangeY);%Rangey needed bfor mode 'projection' |
---|
| 537 | end |
---|
| 538 | % default output |
---|
| 539 | errormsg='';%default |
---|
| 540 | Xline=[]; |
---|
| 541 | flux=0; |
---|
| 542 | circul=0; |
---|
| 543 | liny=ObjectData.Coord(:,2); |
---|
| 544 | NbPoints=size(ObjectData.Coord,1); |
---|
| 545 | testfalse=0; |
---|
| 546 | ListIndex={}; |
---|
| 547 | |
---|
[905] | 548 | %% group the variables (fields of 'FieldData') in cells of variables with the same dimensions |
---|
| 549 | [CellInfo,NbDim,errormsg]=find_field_cells(FieldData); |
---|
| 550 | if ~isempty(errormsg) |
---|
| 551 | errormsg=['error in proj_field/proj_line:' errormsg]; |
---|
| 552 | return |
---|
| 553 | end |
---|
| 554 | CellInfo=CellInfo(NbDim==2); %keep only the 2D cells |
---|
| 555 | %%%%%% TODO: treat 1D fields: project as identity so that P o P=P for projection operation |
---|
| 556 | cell_select=true(size(CellInfo)); |
---|
[854] | 557 | |
---|
[905] | 558 | for icell=1:length(CellInfo) |
---|
| 559 | if isfield(CellInfo{icell},'ProjModeRequest') |
---|
| 560 | if ~strcmp(CellInfo{icell}.ProjModeRequest, ProjMode) |
---|
| 561 | cell_select(icell)=0; |
---|
| 562 | end |
---|
| 563 | if strcmp(ProjMode,'interp_tps')&& ~strcmp(CellInfo{icell}.CoordType,'tps') |
---|
| 564 | cell_select(icell)=0; |
---|
| 565 | end |
---|
| 566 | end |
---|
| 567 | end |
---|
| 568 | if isempty(find(cell_select)) |
---|
| 569 | errormsg=[' invalid projection mode ''' ProjMode ''': use ''interp_tps'' to interpolate spatial derivatives']; |
---|
| 570 | return |
---|
| 571 | end |
---|
| 572 | CellInfo=CellInfo(cell_select); |
---|
| 573 | |
---|
[854] | 574 | %% projection line: object types selected from proj_field='line','polyline','polygon','rectangle','ellipse': |
---|
| 575 | LineCoord=ObjectData.Coord; |
---|
| 576 | switch ObjectData.Type |
---|
| 577 | case 'ellipse' |
---|
| 578 | LineLength=2*pi*ObjectData.RangeX*ObjectData.RangeY; |
---|
| 579 | NbSegment=0; |
---|
| 580 | case 'rectangle' |
---|
| 581 | LineCoord([1 4],1)=ObjectData.Coord(1,1)-ObjectData.RangeX; |
---|
| 582 | LineCoord([1 2],2)=ObjectData.Coord(1,2)-ObjectData.RangeY; |
---|
| 583 | LineCoord([2 3],1)=ObjectData.Coord(1,1)+ObjectData.RangeX; |
---|
| 584 | LineCoord([4 1],2)=ObjectData.Coord(1,2)+ObjectData.RangeY; |
---|
| 585 | case 'polygon' |
---|
| 586 | LineCoord(NbPoints+1)=LineCoord(1); |
---|
| 587 | end |
---|
| 588 | if ~strcmp(ObjectData.Type,'ellipse') |
---|
| 589 | if ~strcmp(ObjectData.Type,'rectangle') && NbPoints<2 |
---|
| 590 | return% line needs at least 2 points to be defined |
---|
| 591 | end |
---|
| 592 | dlinx=diff(LineCoord(:,1)); |
---|
| 593 | dliny=diff(LineCoord(:,2)); |
---|
| 594 | [theta,dlength]=cart2pol(dlinx,dliny);%angle and length of each segment |
---|
| 595 | LineLength=sum(dlength); |
---|
| 596 | NbSegment=numel(LineLength); |
---|
| 597 | end |
---|
| 598 | CheckClosedLine=~isempty(find(strcmp(ObjectData.Type,{'rectangle','ellipse','polygon'}))); |
---|
| 599 | |
---|
| 600 | % x = a \ \cosh \mu \ \cos \nu |
---|
| 601 | % |
---|
| 602 | % y = a \ \sinh \mu \ \sin \nu |
---|
| 603 | |
---|
| 604 | %% angles of the polyline and boundaries of action for mode 'projection' |
---|
| 605 | |
---|
| 606 | % determine a rectangles at +-width from the line (only used for the ProjMode='projection or 'interp_tps') |
---|
| 607 | xsup=zeros(1,NbPoints); xinf=zeros(1,NbPoints); ysup=zeros(1,NbPoints); yinf=zeros(1,NbPoints); |
---|
| 608 | if isequal(ProjMode,'projection') |
---|
| 609 | if strcmp(ObjectData.Type,'line') |
---|
| 610 | xsup=ObjectData.Coord(:,1)-width*sin(theta); |
---|
| 611 | xinf=ObjectData.Coord(:,1)+width*sin(theta); |
---|
| 612 | ysup=ObjectData.Coord(:,2)+width*cos(theta); |
---|
| 613 | yinf=ObjectData.Coord(:,2)-width*cos(theta); |
---|
| 614 | else |
---|
| 615 | errormsg='mode projection only available for simple line, use interpolation otherwise'; |
---|
| 616 | return |
---|
| 617 | end |
---|
| 618 | else % need to define the set of interpolation points |
---|
| 619 | if isfield(ObjectData,'DX') && ~isempty(ObjectData.DX) |
---|
| 620 | DX=abs(ObjectData.DX);%mesh of interpolation points along the line |
---|
| 621 | if CheckClosedLine |
---|
| 622 | NbPoint=ceil(LineLength/DX); |
---|
| 623 | DX=LineLength/NbPoint;%adjust DX to get an integer nbre of intervals in a closed line |
---|
| 624 | DX_edge=DX/2; |
---|
| 625 | else |
---|
| 626 | DX_edge=(LineLength-DX*floor(LineLength/DX))/2;%margin from the first point and first interpolation point, the same for the end point |
---|
| 627 | end |
---|
| 628 | XI=[]; |
---|
| 629 | YI=[]; |
---|
| 630 | ThetaI=[]; |
---|
| 631 | dlengthI=[]; |
---|
| 632 | if strcmp(ObjectData.Type,'ellipse') |
---|
| 633 | phi=(DX_edge:DX:LineLength)*2*pi/LineLength; |
---|
| 634 | XI=ObjectData.RangeX*cos(phi); |
---|
| 635 | YI=ObjectData.RangeY*sin(phi); |
---|
| 636 | dphi=2*pi*DX/LineLength; |
---|
| 637 | [ThetaI,dlengthI]=cart2pol(-ObjectData.RangeX*sin(phi)*dphi,ObjectData.RangeY*cos(phi)*dphi); |
---|
| 638 | else |
---|
| 639 | for isegment=1:NbSegment |
---|
| 640 | costheta=cos(theta(isegment)); |
---|
| 641 | sintheta=sin(theta(isegment)); |
---|
[866] | 642 | % XIsegment=LineCoord(isegment,1)+DX_edge*costheta:DX*costheta:LineCoord(isegment+1,1)); |
---|
| 643 | % YIsegment=(LineCoord(isegment,2)+DX_edge*sintheta:DX*sintheta:LineCoord(isegment+1,2)); |
---|
[854] | 644 | NbInterval=floor((dlength(isegment)-DX_edge)/DX); |
---|
| 645 | LastX=DX_edge+DX*NbInterval; |
---|
| 646 | NbPoint=NbInterval+1; |
---|
| 647 | XIsegment=linspace(LineCoord(isegment,1)+DX_edge*costheta,LineCoord(isegment,1)+LastX*costheta,NbPoint); |
---|
| 648 | YIsegment=linspace(LineCoord(isegment,2)+DX_edge*sintheta,LineCoord(isegment,2)+LastX*sintheta,NbPoint); |
---|
| 649 | XI=[XI XIsegment]; |
---|
| 650 | YI=[YI YIsegment]; |
---|
| 651 | ThetaI=[ThetaI theta(isegment)*ones(1,numel(XIsegment))]; |
---|
| 652 | dlengthI=[dlengthI DX*ones(1,numel(XIsegment))]; |
---|
| 653 | DX_edge=DX-(dlength(isegment)-LastX);%edge for the next segment set to keep DX=DX_end+DX_edge between two segments |
---|
| 654 | end |
---|
| 655 | end |
---|
| 656 | Xproj=cumsum(dlengthI); |
---|
| 657 | else |
---|
| 658 | errormsg='abscissa mesh along line DX needed for interpolation'; |
---|
| 659 | return |
---|
| 660 | end |
---|
| 661 | end |
---|
| 662 | |
---|
| 663 | %% loop on variable cells with the same space dimension 2 |
---|
| 664 | ProjData.ListVarName={}; |
---|
| 665 | ProjData.VarDimName={}; |
---|
[866] | 666 | check_abscissa=0; |
---|
[854] | 667 | for icell=1:length(CellInfo) |
---|
| 668 | % list of variable types to be projected |
---|
| 669 | ListProj={'VarIndex_scalar','VarIndex_image','VarIndex_color','VarIndex_vector_x','VarIndex_vector_y'}; |
---|
| 670 | check_proj=false(size(FieldData.ListVarName)); |
---|
| 671 | for ilist=1:numel(ListProj) |
---|
| 672 | if isfield(CellInfo{icell},ListProj{ilist}) |
---|
| 673 | check_proj(CellInfo{icell}.(ListProj{ilist}))=1; |
---|
| 674 | end |
---|
| 675 | end |
---|
| 676 | VarIndex=find(check_proj);% indices of the variables to be projected |
---|
| 677 | |
---|
| 678 | %% identify vector components |
---|
[867] | 679 | %testU=isfield(CellInfo{icell},'VarIndex_vector_x') &&isfield(CellInfo{icell},'VarIndex_vector_y') ;% test for vectors |
---|
[866] | 680 | % if testU |
---|
| 681 | % UName=FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_x}; |
---|
| 682 | % VName=FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_y}; |
---|
| 683 | % vector_x=FieldData.(UName); |
---|
| 684 | % vector_y=FieldData.(VName); |
---|
| 685 | % end |
---|
[854] | 686 | %identify error flag |
---|
| 687 | errorflag=0; %default, no error flag |
---|
| 688 | if isfield(CellInfo{icell},'VarIndex_errorflag');% test for error flag |
---|
| 689 | FFName=FieldData.ListVarName{CellInfo{icell}.VarIndex_errorflag}; |
---|
| 690 | errorflag=FieldData.(FFName); |
---|
| 691 | end |
---|
| 692 | VarName=FieldData.ListVarName(VarIndex);% cell array of the names of variables to pje |
---|
| 693 | ivar_U=[]; |
---|
| 694 | ivar_V=[]; |
---|
| 695 | %% check needed object properties for unstructured positions (position given by the variables with role coord_x, coord_y |
---|
| 696 | |
---|
| 697 | % circul=0; |
---|
| 698 | % flux=0; |
---|
| 699 | %%%%%%% % A FAIRE CALCULER MEAN DES QUANTITES %%%%%% |
---|
| 700 | switch CellInfo{icell}.CoordType |
---|
| 701 | %case of unstructured coordinates |
---|
| 702 | case 'scattered' |
---|
[866] | 703 | % XName= FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)}; |
---|
| 704 | % YName= FieldData.ListVarName{CellInfo{icell}.CoordIndex(end-1)}; |
---|
| 705 | coord_x=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)}); |
---|
| 706 | coord_y=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(end-1)}); |
---|
| 707 | |
---|
[854] | 708 | if isequal(ProjMode,'projection') |
---|
| 709 | if width==0 |
---|
| 710 | errormsg='range of the projection object is missing'; |
---|
| 711 | return |
---|
| 712 | end |
---|
[867] | 713 | ProjData.ListVarName=[ProjData.ListVarName FieldData.ListVarName(CellInfo{icell}.CoordIndex(end))]; |
---|
| 714 | ProjData.VarDimName=[ProjData.VarDimName FieldData.ListVarName(CellInfo{icell}.CoordIndex(end))]; |
---|
[866] | 715 | nbvar=numel(ProjData.ListVarName); |
---|
| 716 | ProjData.VarAttribute{nbvar}.long_name='abscissa along line'; |
---|
[854] | 717 | % select the (non false) input data located in the band of projection |
---|
| 718 | flagsel=(errorflag==0) & ((coord_y -yinf(1))*(xinf(2)-xinf(1))>(coord_x-xinf(1))*(yinf(2)-yinf(1))) ... |
---|
| 719 | & ((coord_y -ysup(1))*(xsup(2)-xsup(1))<(coord_x-xsup(1))*(ysup(2)-ysup(1))) ... |
---|
| 720 | & ((coord_y -yinf(2))*(xsup(2)-xinf(2))>(coord_x-xinf(2))*(ysup(2)-yinf(2))) ... |
---|
| 721 | & ((coord_y -yinf(1))*(xsup(1)-xinf(1))<(coord_x-xinf(1))*(ysup(1)-yinf(1))); |
---|
| 722 | coord_x=coord_x(flagsel); |
---|
| 723 | coord_y=coord_y(flagsel); |
---|
| 724 | costheta=cos(theta); |
---|
| 725 | sintheta=sin(theta); |
---|
| 726 | Xproj=(coord_x-ObjectData.Coord(1,1))*costheta + (coord_y-ObjectData.Coord(1,2))*sintheta; %projection on the line |
---|
| 727 | [Xproj,indsort]=sort(Xproj);% sort points by increasing absissa along the projection line |
---|
[866] | 728 | ProjData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)})=Xproj; |
---|
[854] | 729 | for ivar=1:numel(VarIndex) |
---|
[867] | 730 | ProjData.(VarName{ivar})=FieldData.(VarName{ivar})(flagsel);% restrict variables to the projection band |
---|
[854] | 731 | ProjData.(VarName{ivar})=ProjData.(VarName{ivar})(indsort);% sort by absissa |
---|
| 732 | ProjData.ListVarName=[ProjData.ListVarName VarName{ivar}]; |
---|
[867] | 733 | ProjData.VarDimName=[ProjData.VarDimName FieldData.ListVarName(CellInfo{icell}.CoordIndex(end))]; |
---|
[854] | 734 | ProjData.VarAttribute{nbvar+ivar}=FieldData.VarAttribute{VarIndex(ivar)};%reproduce var attribute |
---|
| 735 | if isfield(ProjData.VarAttribute{nbvar+ivar},'Role') |
---|
| 736 | if strcmp(ProjData.VarAttribute{nbvar+ivar}.Role,'vector_x'); |
---|
[867] | 737 | ivar_U=nbvar+ivar; |
---|
[854] | 738 | elseif strcmp(ProjData.VarAttribute{nbvar+ivar}.Role,'vector_y'); |
---|
[867] | 739 | ivar_V=nbvar+ivar; |
---|
[854] | 740 | end |
---|
| 741 | end |
---|
| 742 | ProjData.VarAttribute{ivar+nbvar}.Role='discrete';% will promote plots of the profiles with continuous lines |
---|
| 743 | end |
---|
| 744 | elseif isequal(ProjMode,'interp_lin') %filtering %linear interpolation: |
---|
[866] | 745 | if ~check_abscissa |
---|
| 746 | XName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)}; |
---|
| 747 | ProjData.ListVarName=[ProjData.ListVarName {XName}]; |
---|
| 748 | ProjData.VarDimName=[ProjData.VarDimName {XName}]; |
---|
| 749 | nbvar=numel(ProjData.ListVarName); |
---|
| 750 | ProjData.VarAttribute{nbvar}.long_name='abscissa along line'; |
---|
| 751 | check_abscissa=1; % define abcissa only once |
---|
| 752 | end |
---|
[863] | 753 | if ~isequal(errorflag,0) |
---|
| 754 | VarName_FF=FieldData.ListVarName{CellInfo{icell}.VarIndex_errorflag}; |
---|
| 755 | indsel=find(FieldData.(VarName_FF)==0); |
---|
| 756 | coord_x=coord_x(indsel); |
---|
| 757 | coord_y=coord_y(indsel); |
---|
| 758 | for ivar=1:numel(CellInfo{icell}.VarIndex) |
---|
| 759 | VarName=FieldData.ListVarName{CellInfo{icell}.VarIndex(ivar)}; |
---|
| 760 | FieldData.(VarName)=FieldData.(VarName)(indsel); |
---|
| 761 | end |
---|
[866] | 762 | end |
---|
[854] | 763 | [ProjVar,ListFieldProj,VarAttribute,errormsg]=calc_field_interp([coord_x coord_y],FieldData,CellInfo{icell}.FieldName,XI,YI); |
---|
| 764 | ProjData.X=Xproj; |
---|
[866] | 765 | nbvar=numel(ProjData.ListVarName); |
---|
[854] | 766 | ProjData.ListVarName=[ProjData.ListVarName ListFieldProj]; |
---|
| 767 | ProjData.VarAttribute=[ProjData.VarAttribute VarAttribute]; |
---|
| 768 | for ivar=1:numel(VarAttribute) |
---|
| 769 | ProjData.VarDimName=[ProjData.VarDimName {XName}]; |
---|
| 770 | if isfield(VarAttribute{ivar},'Role') |
---|
| 771 | if strcmp(VarAttribute{ivar}.Role,'vector_x'); |
---|
[866] | 772 | ivar_U=ivar+nbvar; |
---|
[854] | 773 | elseif strcmp(VarAttribute{ivar}.Role,'vector_y'); |
---|
[866] | 774 | ivar_V=ivar+nbvar; |
---|
[854] | 775 | end |
---|
| 776 | end |
---|
| 777 | ProjData.VarAttribute{ivar+nbvar}.Role='continuous';% will promote plots of the profiles with continuous lines |
---|
| 778 | ProjData.(ListFieldProj{ivar})=ProjVar{ivar}; |
---|
| 779 | end |
---|
| 780 | end |
---|
| 781 | case 'tps'%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 782 | if strcmp(ProjMode,'interp_tps') |
---|
| 783 | Coord=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex}); |
---|
| 784 | NbCentres=FieldData.(FieldData.ListVarName{CellInfo{icell}.NbCentres_tps}); |
---|
| 785 | SubRange=FieldData.(FieldData.ListVarName{CellInfo{icell}.SubRange_tps}); |
---|
| 786 | if isfield(CellInfo{icell},'VarIndex_vector_x_tps')&&isfield(CellInfo{icell},'VarIndex_vector_y_tps') |
---|
| 787 | FieldVar=cat(3,FieldData.(FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_x_tps}),FieldData.(FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_y_tps})); |
---|
| 788 | end |
---|
| 789 | [DataOut,VarAttribute,errormsg]=calc_field_tps(Coord,NbCentres,SubRange,FieldVar,CellInfo{icell}.FieldName,cat(3,XI,YI)); |
---|
[867] | 790 | ProjData.ListVarName=[ProjData.ListVarName {'X'}]; |
---|
| 791 | ProjData.VarDimName=[ProjData.VarDimName {'X'}]; |
---|
[854] | 792 | ProjData.X=Xproj; |
---|
| 793 | nbvar=numel(ProjData.ListVarName); |
---|
| 794 | ProjData.VarAttribute{nbvar}.long_name='abscissa along line'; |
---|
| 795 | ProjVarName=(fieldnames(DataOut))'; |
---|
| 796 | ProjData.ListVarName=[ProjData.ListVarName ProjVarName]; |
---|
| 797 | ProjData.VarAttribute=[ProjData.VarAttribute VarAttribute]; |
---|
| 798 | for ivar=1:numel(VarAttribute) |
---|
[867] | 799 | ProjData.VarDimName=[ProjData.VarDimName {'X'}]; |
---|
[854] | 800 | if isfield(VarAttribute{ivar},'Role') |
---|
| 801 | if strcmp(VarAttribute{ivar}.Role,'vector_x'); |
---|
[866] | 802 | ivar_U=ivar+nbvar; |
---|
[854] | 803 | elseif strcmp(VarAttribute{ivar}.Role,'vector_y'); |
---|
[866] | 804 | ivar_V=ivar+nbvar; |
---|
[854] | 805 | end |
---|
| 806 | end |
---|
| 807 | ProjData.VarAttribute{ivar+nbvar}.Role='continuous';% will promote plots of the profiles with continuous lines |
---|
| 808 | ProjData.(ProjVarName{ivar})=DataOut.(ProjVarName{ivar}); |
---|
| 809 | end |
---|
| 810 | end |
---|
| 811 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 812 | |
---|
| 813 | case 'grid' %case of structured coordinates |
---|
| 814 | if ~isequal(ObjectData.Type,'line')% exclude polyline |
---|
| 815 | errormsg=['no projection available on ' ObjectData.Type 'for structured coordinates']; % |
---|
| 816 | else |
---|
| 817 | test_Amat=1;%image or 2D matrix |
---|
| 818 | test_interp2=0;%default |
---|
| 819 | AYName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(end-1)}; |
---|
| 820 | AXName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)}; |
---|
| 821 | eval(['AX=FieldData.' AXName ';']);% set of x positions |
---|
| 822 | eval(['AY=FieldData.' AYName ';']);% set of y positions |
---|
| 823 | AName=FieldData.ListVarName{VarIndex(1)}; |
---|
| 824 | eval(['A=FieldData.' AName ';']);% scalar |
---|
| 825 | npxy=size(A); |
---|
| 826 | npx=npxy(2); |
---|
| 827 | npy=npxy(1); |
---|
| 828 | if numel(AX)==2 |
---|
| 829 | DX=(AX(2)-AX(1))/(npx-1); |
---|
| 830 | else |
---|
| 831 | DX_vec=diff(AX); |
---|
| 832 | DX=max(DX_vec); |
---|
| 833 | DX_min=min(DX_vec); |
---|
| 834 | if (DX-DX_min)>0.0001*abs(DX) |
---|
| 835 | test_interp2=1; |
---|
| 836 | DX=DX_min; |
---|
| 837 | end |
---|
| 838 | end |
---|
| 839 | if numel(AY)==2 |
---|
| 840 | DY=(AY(2)-AY(1))/(npy-1); |
---|
| 841 | else |
---|
| 842 | DY_vec=diff(AY); |
---|
| 843 | DY=max(DY_vec); |
---|
| 844 | DY_min=min(DY_vec); |
---|
| 845 | if (DY-DY_min)>0.0001*abs(DY) |
---|
| 846 | test_interp2=1; |
---|
| 847 | DY=DY_min; |
---|
| 848 | end |
---|
| 849 | end |
---|
| 850 | AXI=linspace(AX(1),AX(end), npx);%set of x positions for the interpolated input data |
---|
| 851 | AYI=linspace(AY(1),AY(end), npy);%set of x positions for the interpolated input data |
---|
| 852 | if isfield(ObjectData,'DX') |
---|
| 853 | DXY_line=ObjectData.DX;%mesh on the projection line |
---|
| 854 | else |
---|
| 855 | DXY_line=sqrt(abs(DX*DY));% mesh on the projection line |
---|
| 856 | end |
---|
| 857 | dlinx=ObjectData.Coord(2,1)-ObjectData.Coord(1,1); |
---|
| 858 | dliny=ObjectData.Coord(2,2)-ObjectData.Coord(1,2); |
---|
| 859 | linelength=sqrt(dlinx*dlinx+dliny*dliny); |
---|
| 860 | theta=angle(dlinx+i*dliny);%angle of the line |
---|
| 861 | if isfield(FieldData,'RangeX') |
---|
| 862 | XMin=min(FieldData.RangeX);%shift of the origin on the line |
---|
| 863 | else |
---|
| 864 | XMin=0; |
---|
| 865 | end |
---|
| 866 | eval(['ProjData.' AXName '=linspace(XMin,XMin+linelength,linelength/DXY_line+1);'])%abscissa of the new pixels along the line |
---|
| 867 | y=linspace(-width,width,2*width/DXY_line+1);%ordintes of the new pixels (coordinate across the line) |
---|
| 868 | eval(['npX=length(ProjData.' AXName ');']) |
---|
| 869 | npY=length(y); %TODO: utiliser proj_grid |
---|
| 870 | eval(['[X,Y]=meshgrid(ProjData.' AXName ',y);'])%grid in the line coordinates |
---|
| 871 | XIMA=ObjectData.Coord(1,1)+(X-XMin)*cos(theta)-Y*sin(theta); |
---|
| 872 | YIMA=ObjectData.Coord(1,2)+(X-XMin)*sin(theta)+Y*cos(theta); |
---|
| 873 | XIMA=(XIMA-AX(1))/DX+1;% index of the original image along x |
---|
| 874 | YIMA=(YIMA-AY(1))/DY+1;% index of the original image along y |
---|
| 875 | XIMA=reshape(round(XIMA),1,npX*npY);%indices reorganized in 'line' |
---|
| 876 | YIMA=reshape(round(YIMA),1,npX*npY); |
---|
| 877 | flagin=XIMA>=1 & XIMA<=npx & YIMA >=1 & YIMA<=npy;%flagin=1 inside the original image |
---|
| 878 | ind_in=find(flagin); |
---|
| 879 | ind_out=find(~flagin); |
---|
| 880 | ICOMB=(XIMA-1)*npy+YIMA; |
---|
| 881 | ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A |
---|
| 882 | nbcolor=1; %color images |
---|
| 883 | if numel(npxy)==2 |
---|
| 884 | nbcolor=1; |
---|
| 885 | elseif length(npxy)==3 |
---|
| 886 | nbcolor=npxy(3); |
---|
| 887 | else |
---|
| 888 | errormsg='multicomponent field not projected'; |
---|
| 889 | display(errormsg) |
---|
| 890 | return |
---|
| 891 | end |
---|
| 892 | nbvar=length(ProjData.ListVarName);% number of var from previous cells |
---|
| 893 | ProjData.ListVarName=[ProjData.ListVarName {AXName}]; |
---|
| 894 | ProjData.VarDimName=[ProjData.VarDimName {AXName}]; |
---|
| 895 | for ivar=VarIndex |
---|
| 896 | %VarName{ivar}=FieldData.ListVarName{ivar}; |
---|
| 897 | if test_interp2% interpolate on new grid |
---|
| 898 | FieldData.(FieldData.ListVarName{ivar})=interp2(FieldData.(AXName),FieldData.(AYName),FieldData.(FieldData.ListVarName{ivar}),AXI,AYI);%TO TEST |
---|
| 899 | end |
---|
| 900 | vec_A=reshape(squeeze(FieldData.(FieldData.ListVarName{ivar})),npx*npy,nbcolor); %put the original image in colum |
---|
| 901 | if nbcolor==1 |
---|
| 902 | vec_B(ind_in)=vec_A(ICOMB); |
---|
| 903 | vec_B(ind_out)=zeros(size(ind_out)); |
---|
| 904 | A_out=reshape(vec_B,npY,npX); |
---|
| 905 | ProjData.(FieldData.ListVarName{ivar}) =sum(A_out,1)/npY; |
---|
| 906 | elseif nbcolor==3 |
---|
| 907 | vec_B(ind_in,1:3)=vec_A(ICOMB,:); |
---|
| 908 | vec_B(ind_out,1)=zeros(size(ind_out)); |
---|
| 909 | vec_B(ind_out,2)=zeros(size(ind_out)); |
---|
| 910 | vec_B(ind_out,3)=zeros(size(ind_out)); |
---|
| 911 | A_out=reshape(vec_B,npY,npX,nbcolor); |
---|
| 912 | ProjData.(FieldData.ListVarName{ivar})=squeeze(sum(A_out,1)/npY); |
---|
| 913 | end |
---|
| 914 | ProjData.ListVarName=[ProjData.ListVarName FieldData.ListVarName{ivar}]; |
---|
| 915 | ProjData.VarDimName=[ProjData.VarDimName {AXName}];%to generalize with the initial name of the x coordinate |
---|
| 916 | ProjData.VarAttribute{ivar}.Role='continuous';% for plot with continuous line |
---|
| 917 | end |
---|
| 918 | if nbcolor==3 |
---|
| 919 | ProjData.VarDimName{end}={AXName,'rgb'}; |
---|
| 920 | end |
---|
| 921 | end |
---|
[866] | 922 | end |
---|
[854] | 923 | if ~isempty(ivar_U) && ~isempty(ivar_V) |
---|
[866] | 924 | vector_x =ProjData.(ProjData.ListVarName{ivar_U}); |
---|
| 925 | ProjData.(ProjData.ListVarName{ivar_U}) =cos(theta)*vector_x+sin(theta)*ProjData.(ProjData.ListVarName{ivar_V}); |
---|
| 926 | ProjData.(ProjData.ListVarName{ivar_V}) =-sin(theta)*vector_x+cos(theta)*ProjData.(ProjData.ListVarName{ivar_V}); |
---|
[854] | 927 | end |
---|
| 928 | end |
---|
| 929 | |
---|
| 930 | % %shotarter case for horizontal or vertical line (A FAIRE |
---|
| 931 | % % Rangx=[0.5 npx-0.5];%image coordiantes of corners |
---|
| 932 | % % Rangy=[npy-0.5 0.5]; |
---|
| 933 | % % if isfield(Calib,'Pxcmx')&isfield(Calib,'Pxcmy')%old calib |
---|
| 934 | % % Rangx=Rangx/Calib.Pxcmx; |
---|
| 935 | % % Rangy=Rangy/Calib.Pxcmy; |
---|
| 936 | % % else |
---|
| 937 | % % [Rangx]=phys_XYZ(Calib,Rangx,[0.5 0.5],[0 0]);%case of translations without rotation and quadratic deformation |
---|
| 938 | % % [xx,Rangy]=phys_XYZ(Calib,[0.5 0.5],Rangy,[0 0]); |
---|
| 939 | % % end |
---|
| 940 | % |
---|
| 941 | % % test_scal=0;%default% 3- 'UserData':(get(handles.Tag,'UserData') |
---|
| 942 | |
---|
| 943 | |
---|
| 944 | %----------------------------------------------------------------- |
---|
| 945 | %project on a plane |
---|
| 946 | % AJOUTER flux,circul,error |
---|
| 947 | function [ProjData,errormsg] = proj_plane(FieldData, ObjectData) |
---|
| 948 | %----------------------------------------------------------------- |
---|
| 949 | |
---|
| 950 | %% rotation angles |
---|
| 951 | PlaneAngle=[0 0 0]; |
---|
| 952 | norm_plane=[0 0 1]; |
---|
[935] | 953 | %cos_om=1; |
---|
| 954 | %sin_om=0; |
---|
[854] | 955 | test90x=0;%=1 for 90 degree rotation alround x axis |
---|
| 956 | test90y=0;%=1 for 90 degree rotation alround y axis |
---|
| 957 | if isfield(ObjectData,'Angle')&& isequal(size(ObjectData.Angle),[1 3])&& ~isequal(ObjectData.Angle,[0 0 0]) |
---|
| 958 | test90y=isequal(ObjectData.Angle,[0 90 0]); |
---|
| 959 | PlaneAngle=(pi/180)*ObjectData.Angle; |
---|
| 960 | om=norm(PlaneAngle);%norm of rotation angle in radians |
---|
| 961 | OmAxis=PlaneAngle/om; %unit vector marking the rotation axis |
---|
| 962 | cos_om=cos(om); |
---|
| 963 | sin_om=sin(om); |
---|
| 964 | coeff=OmAxis(3)*(1-cos_om); |
---|
| 965 | %components of the unity vector norm_plane normal to the projection plane |
---|
| 966 | norm_plane(1)=OmAxis(1)*coeff+OmAxis(2)*sin_om; |
---|
| 967 | norm_plane(2)=OmAxis(2)*coeff-OmAxis(1)*sin_om; |
---|
| 968 | norm_plane(3)=OmAxis(3)*coeff+cos_om; |
---|
| 969 | end |
---|
| 970 | testangle=~isequal(PlaneAngle,[0 0 0])||~isequal(ObjectData.Coord(1:2),[0 0 ]) ;% && ~test90y && ~test90x;%=1 for slanted plane |
---|
| 971 | |
---|
| 972 | %% mesh sizes DX and DY |
---|
| 973 | DX=[]; |
---|
| 974 | DY=[];%default |
---|
| 975 | if isfield(ObjectData,'DX') && ~isempty(ObjectData.DX) |
---|
| 976 | DX=abs(ObjectData.DX);%mesh of interpolation points |
---|
| 977 | elseif isfield(FieldData,'CoordMesh') |
---|
| 978 | DX=FieldData.CoordMesh; |
---|
| 979 | end |
---|
| 980 | if isfield(ObjectData,'DY') && ~isempty(ObjectData.DY) |
---|
| 981 | DY=abs(ObjectData.DY);%mesh of interpolation points |
---|
| 982 | elseif isfield(FieldData,'CoordMesh') |
---|
| 983 | DY=FieldData.CoordMesh; |
---|
| 984 | end |
---|
| 985 | if ~strcmp(ObjectData.ProjMode,'projection') && (isempty(DX)||isempty(DY)) |
---|
| 986 | errormsg='DX or DY not defined'; |
---|
| 987 | return |
---|
| 988 | end |
---|
| 989 | |
---|
| 990 | %% extrema along each axis |
---|
| 991 | testXMin=0;% test if min of X coordinates defined on the projection object, =0 by default |
---|
| 992 | testXMax=0;% test if max of X coordinates defined on the projection object, =0 by default |
---|
| 993 | testYMin=0;% test if min of Y coordinates defined on the projection object, =0 by default |
---|
| 994 | testYMax=0;% test if max of Y coordinates defined on the projection object, =0 by default |
---|
| 995 | if isfield(ObjectData,'RangeX') % rangeX defined by the projection object |
---|
| 996 | XMin=min(ObjectData.RangeX); |
---|
| 997 | XMax=max(ObjectData.RangeX); |
---|
| 998 | testXMin=XMax>XMin;%=1 if XMin defined (i.e. RangeY has two distinct elements) |
---|
| 999 | testXMax=1;% max of X coordinates defined on the projection object |
---|
| 1000 | end |
---|
| 1001 | if isfield(ObjectData,'RangeY') % rangeY defined by the projection object |
---|
| 1002 | YMin=min(ObjectData.RangeY); |
---|
| 1003 | YMax=max(ObjectData.RangeY); |
---|
| 1004 | testYMin=YMax>YMin;%=1 if YMin defined (i.e. RangeY has tow distinct elements) |
---|
| 1005 | testYMax=1;% max of Y coordinates defined on the projection object |
---|
| 1006 | end |
---|
| 1007 | width=0;%default width of the projection band |
---|
| 1008 | if isfield(ObjectData,'RangeZ') |
---|
| 1009 | width=max(ObjectData.RangeZ); |
---|
| 1010 | end |
---|
| 1011 | |
---|
[912] | 1012 | %% interpolation range |
---|
| 1013 | thresh2=[]; |
---|
| 1014 | if isfield(ObjectData,'RangeInterp') |
---|
| 1015 | thresh2=ObjectData.RangeInterp*ObjectData.RangeInterp;%square of interpolation range (do not interpolate beyond this range) |
---|
| 1016 | end |
---|
| 1017 | |
---|
[854] | 1018 | %% initiate Matlab structure for physical field |
---|
| 1019 | [ProjData,errormsg]=proj_heading(FieldData,ObjectData); |
---|
| 1020 | if ~isempty(errormsg) |
---|
| 1021 | return |
---|
| 1022 | end |
---|
| 1023 | |
---|
| 1024 | %% reproduce initial plane position and angle |
---|
| 1025 | if isfield(FieldData,'PlaneCoord')&&length(FieldData.PlaneCoord)==3&& isfield(ProjData,'ProjObjectCoord') |
---|
| 1026 | if length(ProjData.ProjObjectCoord)==3% if the projection plane has a z coordinate |
---|
[864] | 1027 | if isfield(ProjData,'.PlaneCoord') && ~isequal(ProjData.PlaneCoord(3),ProjData.ProjObjectCoord) %check the consistency with the z coordinate of the field plane (set by calibration) |
---|
[854] | 1028 | errormsg='inconsistent z position for field and projection plane'; |
---|
| 1029 | return |
---|
| 1030 | end |
---|
| 1031 | else % the z coordinate is set only by the field plane (by calibration) |
---|
| 1032 | ProjData.ProjObjectCoord(3)=FieldData.PlaneCoord(3); |
---|
| 1033 | end |
---|
| 1034 | if isfield(FieldData,'PlaneAngle') |
---|
| 1035 | if isfield(ProjData,'ProjObjectAngle') |
---|
| 1036 | if ~isequal(FieldData.PlaneAngle,ProjData.ProjObjectAngle) %check the consistency with the z coordinate of the field plane (set by calibration) |
---|
| 1037 | errormsg='inconsistent plane angle for field and projection plane'; |
---|
| 1038 | return |
---|
| 1039 | end |
---|
| 1040 | else |
---|
| 1041 | ProjData.ProjObjectAngle=FieldData.PlaneAngle; |
---|
| 1042 | end |
---|
| 1043 | end |
---|
| 1044 | end |
---|
| 1045 | ProjData.NbDim=2; |
---|
| 1046 | ProjData.ListVarName={}; |
---|
| 1047 | ProjData.VarDimName={}; |
---|
| 1048 | ProjData.VarAttribute={}; |
---|
| 1049 | if ~isempty(DX) && ~isempty(DY) |
---|
| 1050 | ProjData.CoordMesh=sqrt(DX*DY);%define typical data mesh, useful for mouse selection in plots |
---|
| 1051 | elseif isfield(FieldData,'CoordMesh') |
---|
| 1052 | ProjData.CoordMesh=FieldData.CoordMesh; |
---|
| 1053 | end |
---|
[935] | 1054 | %error=0;%default |
---|
| 1055 | %flux=0; |
---|
| 1056 | %testfalse=0; |
---|
| 1057 | %ListIndex={}; |
---|
[854] | 1058 | |
---|
| 1059 | %% group the variables (fields of 'FieldData') in cells of variables with the same dimensions |
---|
| 1060 | [CellInfo,NbDimArray,errormsg]=find_field_cells(FieldData); |
---|
| 1061 | |
---|
| 1062 | if ~isempty(errormsg) |
---|
| 1063 | errormsg=['error in proj_field/proj_plane:' errormsg]; |
---|
| 1064 | return |
---|
| 1065 | end |
---|
| 1066 | check_grid=zeros(size(CellInfo));% =1 if a grid is needed , =0 otherwise, for each field cell |
---|
| 1067 | |
---|
| 1068 | ProjMode=cell(size(CellInfo)); |
---|
| 1069 | for icell=1:numel(CellInfo) |
---|
| 1070 | ProjMode{icell}=ObjectData.ProjMode;% projection mode of the plane object |
---|
| 1071 | end |
---|
| 1072 | icell_grid=[];% field cell index which defines the grid |
---|
| 1073 | if ~strcmp(ObjectData.ProjMode,'projection') |
---|
| 1074 | %% define the new coordinates in case of interpolation on a imposed grid |
---|
| 1075 | if ~testYMin |
---|
| 1076 | errormsg='min Y value not defined for the projection grid';return |
---|
| 1077 | end |
---|
| 1078 | if ~testYMax |
---|
| 1079 | errormsg='max Y value not defined for the projection grid';return |
---|
| 1080 | end |
---|
| 1081 | if ~testXMin |
---|
| 1082 | errormsg='min X value not defined for the projection grid';return |
---|
| 1083 | end |
---|
| 1084 | if ~testXMax |
---|
| 1085 | errormsg='max X value not defined for the projection grid';return |
---|
| 1086 | end |
---|
| 1087 | else |
---|
| 1088 | %% case of a grid requested by the input field |
---|
| 1089 | for icell=1:numel(CellInfo)% TODO: recalculate coordinates here to get the bounds in the rotated coordinates |
---|
| 1090 | if isfield(CellInfo{icell},'ProjModeRequest') |
---|
| 1091 | switch CellInfo{icell}.ProjModeRequest |
---|
| 1092 | case 'interp_lin' |
---|
| 1093 | ProjMode{icell}='interp_lin'; |
---|
| 1094 | case 'interp_tps' |
---|
| 1095 | ProjMode{icell}='interp_tps'; |
---|
| 1096 | end |
---|
| 1097 | end |
---|
| 1098 | if strcmp(ProjMode{icell},'interp_lin')||strcmp(ProjMode{icell},'interp_tps') |
---|
| 1099 | check_grid(icell)=1; |
---|
| 1100 | end |
---|
| 1101 | if strcmp(CellInfo{icell}.CoordType,'grid')&&NbDimArray(icell)>=2 |
---|
| 1102 | if ~testangle && isempty(icell_grid)% if the input gridded data is not modified, choose the first one in case of multiple gridded field cells |
---|
| 1103 | icell_grid=icell; |
---|
| 1104 | ProjMode{icell}='projection'; |
---|
| 1105 | end |
---|
| 1106 | check_grid(icell)=1; |
---|
| 1107 | end |
---|
| 1108 | end |
---|
| 1109 | if ~isempty(find(check_grid))% if a grid is requested by the input field |
---|
| 1110 | if isempty(icell_grid)% if the grid is not given by cell #icell_grid |
---|
| 1111 | if ~isfield(FieldData,'XMax') |
---|
| 1112 | FieldData=find_field_bounds(FieldData); |
---|
| 1113 | end |
---|
| 1114 | end |
---|
| 1115 | end |
---|
| 1116 | end |
---|
| 1117 | if ~isempty(find(check_grid))||~strcmp(ObjectData.ProjMode,'projection')%no existing gridded data used |
---|
| 1118 | if isempty(icell_grid)||~strcmp(ObjectData.ProjMode,'projection')%no existing gridded data used |
---|
| 1119 | AYName='coord_y'; |
---|
| 1120 | AXName='coord_x'; |
---|
| 1121 | if strcmp(ObjectData.ProjMode,'projection') |
---|
| 1122 | ProjData.coord_y=[FieldData.YMin FieldData.YMax];%note that if projection is done on a grid, the Min and Max along each direction must have been defined |
---|
| 1123 | ProjData.coord_x=[FieldData.XMin FieldData.XMax]; |
---|
| 1124 | coord_x_proj=FieldData.XMin:FieldData.CoordMesh:FieldData.XMax; |
---|
| 1125 | coord_y_proj=FieldData.YMin:FieldData.CoordMesh:FieldData.YMax; |
---|
| 1126 | else |
---|
| 1127 | ProjData.coord_y=[ObjectData.RangeY(1) ObjectData.RangeY(2)];%note that if projection is done on a grid, the Min and Max along each direction must have been defined |
---|
| 1128 | ProjData.coord_x=[ObjectData.RangeX(1) ObjectData.RangeX(2)]; |
---|
| 1129 | coord_x_proj=ObjectData.RangeX(1):ObjectData.DX:ObjectData.RangeX(2); |
---|
| 1130 | coord_y_proj=ObjectData.RangeY(1):ObjectData.DY:ObjectData.RangeY(2); |
---|
| 1131 | end |
---|
| 1132 | [XI,YI]=meshgrid(coord_x_proj,coord_y_proj);%grid in the new coordinates |
---|
[890] | 1133 | ProjData.VarDimName={AYName,AXName}; |
---|
[854] | 1134 | % XI=ObjectData.Coord(1,1)+(X)*cos(PlaneAngle(3))-YI*sin(PlaneAngle(3));%corresponding coordinates in the original system |
---|
| 1135 | % YI=ObjectData.Coord(1,2)+(X)*sin(PlaneAngle(3))+YI*cos(PlaneAngle(3)); |
---|
| 1136 | else% we use the existing grid from field cell #icell_grid |
---|
| 1137 | NbDim=NbDimArray(icell_grid); |
---|
| 1138 | AYName=FieldData.ListVarName{CellInfo{icell_grid}.CoordIndex(NbDim-1)};%name of input x coordinate (name preserved on projection) |
---|
| 1139 | AXName=FieldData.ListVarName{CellInfo{icell_grid}.CoordIndex(NbDim)};%name of input y coordinate (name preserved on projection) |
---|
[890] | 1140 | AYDimName=FieldData.VarDimName{CellInfo{icell_grid}.CoordIndex(NbDim-1)};% |
---|
| 1141 | AXDimName=FieldData.VarDimName{CellInfo{icell_grid}.CoordIndex(NbDim)};% |
---|
| 1142 | ProjData.VarDimName={AYDimName,AXDimName}; |
---|
[854] | 1143 | ProjData.(AYName)=FieldData.(AYName); % new (projected ) y coordinates |
---|
| 1144 | ProjData.(AXName)=FieldData.(AXName); % new (projected ) y coordinates |
---|
| 1145 | end |
---|
| 1146 | ProjData.ListVarName={AYName,AXName}; |
---|
[890] | 1147 | |
---|
[854] | 1148 | ProjData.VarAttribute={[],[]}; |
---|
| 1149 | end |
---|
| 1150 | |
---|
| 1151 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 1152 | % LOOP ON FIELD CELLS, PROJECT VARIABLES |
---|
| 1153 | % CellVarIndex=cells of variable index arrays |
---|
| 1154 | %ivar_new=0; % index of the current variable in the projected field |
---|
| 1155 | % icoord=0; |
---|
| 1156 | nbcoord=0;%number of added coordinate variables brought by projection |
---|
| 1157 | %nbvar=0; |
---|
| 1158 | vector_x_proj=[]; |
---|
| 1159 | vector_y_proj=[]; |
---|
| 1160 | for icell=1:length(CellInfo) |
---|
| 1161 | NbDim=NbDimArray(icell); |
---|
| 1162 | if NbDim<2 |
---|
| 1163 | continue % only cells represnting 2D or 3D fields are involved |
---|
| 1164 | end |
---|
| 1165 | VarIndex=CellInfo{icell}.VarIndex;% indices of the selected variables in the list FieldData.ListVarName |
---|
| 1166 | %dimensions |
---|
| 1167 | DimCell=FieldData.VarDimName{VarIndex(1)}; |
---|
| 1168 | if ischar(DimCell) |
---|
| 1169 | DimCell={DimCell};%name of dimensions |
---|
| 1170 | end |
---|
| 1171 | coord_z=0;%default |
---|
| 1172 | ListVarName={};% initiate list of projected variables for cell # icell |
---|
| 1173 | VarDimName={};% initiate coresponding list of dimensions for cell # icell |
---|
| 1174 | VarAttribute={};% initiate coresponding list of var attributes for cell # icell |
---|
| 1175 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 1176 | switch CellInfo{icell}.CoordType |
---|
| 1177 | |
---|
| 1178 | case 'scattered' |
---|
| 1179 | %% case of input fields with unstructured coordinates (applies for projMode ='projection' or 'interp_lin') |
---|
| 1180 | if strcmp(ProjMode{icell},'interp_tps') |
---|
| 1181 | continue %skip for next cell (needs tps field cell) |
---|
| 1182 | end |
---|
| 1183 | coord_x=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)});% initial x coordinates |
---|
| 1184 | coord_y=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(end-1)});% initial y coordinates |
---|
| 1185 | check3D=(numel(CellInfo{icell}.CoordIndex)==3); |
---|
| 1186 | if check3D |
---|
| 1187 | coord_z=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(1)}); |
---|
| 1188 | end |
---|
| 1189 | |
---|
| 1190 | % translate initial coordinates to account for the new origin |
---|
| 1191 | coord_x=coord_x-ObjectData.Coord(1,1); |
---|
| 1192 | coord_y=coord_y-ObjectData.Coord(1,2); |
---|
| 1193 | if check3D |
---|
| 1194 | coord_z=coord_z-ObjectData.Coord(1,3); |
---|
| 1195 | end |
---|
| 1196 | |
---|
| 1197 | % selection of the vectors in the projection range (3D case) |
---|
| 1198 | if check3D && width > 0 |
---|
| 1199 | %components of the unitiy vector normal to the projection plane |
---|
| 1200 | fieldZ=norm_plane(1)*coord_x + norm_plane(2)*coord_y+ norm_plane(3)*coord_z;% distance to the plane |
---|
| 1201 | indcut=find(abs(fieldZ) <= width); |
---|
| 1202 | for ivar=VarIndex |
---|
| 1203 | VarName=FieldData.ListVarName{ivar}; |
---|
| 1204 | FieldData.(VarName)=FieldData.(VarName)(indcut); |
---|
| 1205 | end |
---|
| 1206 | coord_x=coord_x(indcut); |
---|
| 1207 | coord_y=coord_y(indcut); |
---|
| 1208 | coord_z=coord_z(indcut); |
---|
| 1209 | end |
---|
| 1210 | |
---|
| 1211 | %rotate coordinates if needed: coord_X,coord_Y= = coordinates in the new plane |
---|
| 1212 | Psi=PlaneAngle(1); |
---|
| 1213 | Theta=PlaneAngle(2); |
---|
| 1214 | Phi=PlaneAngle(3); |
---|
| 1215 | if testangle && ~test90y && ~test90x;%=1 for slanted plane |
---|
| 1216 | coord_X=(coord_x *cos(Phi) + coord_y* sin(Phi)); |
---|
| 1217 | coord_Y=(-coord_x *sin(Phi) + coord_y *cos(Phi))*cos(Theta); |
---|
| 1218 | coord_Y=coord_Y+coord_z *sin(Theta); |
---|
| 1219 | coord_X=(coord_X *cos(Psi) - coord_Y* sin(Psi));%A VERIFIER |
---|
| 1220 | coord_Y=(coord_X *sin(Psi) + coord_Y* cos(Psi)); |
---|
| 1221 | else |
---|
| 1222 | coord_X=coord_x; |
---|
| 1223 | coord_Y=coord_y; |
---|
| 1224 | end |
---|
| 1225 | |
---|
| 1226 | %restriction to the range of X and Y if imposed by the projection object |
---|
| 1227 | testin=ones(size(coord_X)); %default |
---|
| 1228 | testbound=0; |
---|
| 1229 | if testXMin |
---|
| 1230 | testin=testin & (coord_X >= XMin); |
---|
| 1231 | testbound=1; |
---|
| 1232 | end |
---|
| 1233 | if testXMax |
---|
| 1234 | testin=testin & (coord_X <= XMax); |
---|
| 1235 | testbound=1; |
---|
| 1236 | end |
---|
| 1237 | if testYMin |
---|
| 1238 | testin=testin & (coord_Y >= YMin); |
---|
| 1239 | testbound=1; |
---|
| 1240 | end |
---|
| 1241 | if testYMin |
---|
| 1242 | testin=testin & (coord_Y <= YMax); |
---|
| 1243 | testbound=1; |
---|
| 1244 | end |
---|
| 1245 | if testbound |
---|
| 1246 | indcut=find(testin); |
---|
| 1247 | if isempty(indcut) |
---|
| 1248 | errormsg='data outside the bounds of the projection object'; |
---|
| 1249 | return |
---|
| 1250 | end |
---|
| 1251 | for ivar=VarIndex |
---|
| 1252 | VarName=FieldData.ListVarName{ivar}; |
---|
| 1253 | FieldData.(VarName)=FieldData.(VarName)(indcut); |
---|
| 1254 | end |
---|
| 1255 | coord_X=coord_X(indcut); |
---|
| 1256 | coord_Y=coord_Y(indcut); |
---|
| 1257 | if check3D |
---|
| 1258 | coord_Z=coord_Z(indcut); |
---|
| 1259 | end |
---|
| 1260 | end |
---|
| 1261 | |
---|
| 1262 | % two cases of projection for scattered coordinates |
---|
| 1263 | switch ProjMode{icell} |
---|
| 1264 | case 'projection' |
---|
| 1265 | nbvar=0; |
---|
| 1266 | %nbvar=numel(ProjData.ListVarName); |
---|
| 1267 | for ivar=VarIndex %transfer variables to the projection plane |
---|
| 1268 | VarName=FieldData.ListVarName{ivar}; |
---|
| 1269 | if ivar==CellInfo{icell}.CoordIndex(end) |
---|
| 1270 | ProjData.(VarName)=coord_X; |
---|
| 1271 | elseif ivar==CellInfo{icell}.CoordIndex(end-1) % y coordinate |
---|
| 1272 | ProjData.(VarName)=coord_Y; |
---|
| 1273 | elseif ~(check3D && ivar==CellInfo{icell}.CoordIndex(1)) % other variables (except Z coordinate wyhich is not reproduced) |
---|
| 1274 | ProjData.(VarName)=FieldData.(VarName); |
---|
| 1275 | end |
---|
| 1276 | if ~(check3D && ivar==CellInfo{icell}.CoordIndex(1)) |
---|
| 1277 | ListVarName=[ListVarName VarName]; |
---|
| 1278 | VarDimName=[VarDimName DimCell]; |
---|
| 1279 | nbvar=nbvar+1; |
---|
| 1280 | if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute) >=ivar |
---|
| 1281 | VarAttribute{nbvar}=FieldData.VarAttribute{ivar}; |
---|
| 1282 | end |
---|
| 1283 | end |
---|
| 1284 | end |
---|
| 1285 | case 'interp_lin'%interpolate data on a regular grid |
---|
| 1286 | if isfield(CellInfo{icell},'VarIndex_errorflag') |
---|
| 1287 | VarName_FF=FieldData.ListVarName{CellInfo{icell}.VarIndex_errorflag}; |
---|
| 1288 | indsel=find(FieldData.(VarName_FF)==0); |
---|
| 1289 | coord_X=coord_X(indsel); |
---|
| 1290 | coord_Y=coord_Y(indsel); |
---|
| 1291 | for ivar=1:numel(CellInfo{icell}.VarIndex) |
---|
| 1292 | VarName=FieldData.ListVarName{CellInfo{icell}.VarIndex(ivar)}; |
---|
| 1293 | FieldData.(VarName)=FieldData.(VarName)(indsel); |
---|
| 1294 | end |
---|
| 1295 | end |
---|
| 1296 | % interpolate and calculate field on the grid |
---|
[888] | 1297 | |
---|
[854] | 1298 | [VarVal,ListVarName,VarAttribute,errormsg]=calc_field_interp([coord_X coord_Y],FieldData,CellInfo{icell}.FieldName,XI,YI); |
---|
| 1299 | |
---|
[895] | 1300 | % set to NaN interpolation points which are too far from any initial data (more than 2 CoordMesh) |
---|
| 1301 | if exist('scatteredInterpolant','file')%recent Matlab versions |
---|
| 1302 | F=scatteredInterpolant(coord_X, coord_Y,coord_X,'nearest'); |
---|
| 1303 | G=scatteredInterpolant(coord_X, coord_Y,coord_Y,'nearest'); |
---|
| 1304 | else |
---|
| 1305 | F=TriScatteredInterp([coord_X coord_Y],coord_X,'nearest'); |
---|
| 1306 | G=TriScatteredInterp([coord_X coord_Y],coord_Y,'nearest'); |
---|
| 1307 | end |
---|
[888] | 1308 | Distx=F(XI,YI)-XI;% diff of x coordinates with the nearest measurement point |
---|
| 1309 | Disty=G(XI,YI)-YI;% diff of y coordinates with the nearest measurement point |
---|
| 1310 | Dist=Distx.*Distx+Disty.*Disty; |
---|
[912] | 1311 | if ~isempty(thresh2) |
---|
| 1312 | for ivar=1:numel(VarVal) |
---|
| 1313 | VarVal{ivar}(Dist>thresh2)=NaN;% % put to NaN interpolated positions further than 4 meshes from initial data |
---|
| 1314 | end |
---|
| 1315 | end |
---|
[854] | 1316 | if isfield(CellInfo{icell},'CheckSub') && CellInfo{icell}.CheckSub && ~isempty(vector_x_proj) |
---|
[866] | 1317 | ProjData.(FieldData.ListVarName{vector_x_proj})=ProjData.(FieldData.ListVarName{vector_x_proj})-VarVal{1}; |
---|
| 1318 | ProjData.(FieldData.ListVarName{vector_y_proj})=ProjData.(FieldData.ListVarName{vector_y_proj})-VarVal{2}; |
---|
| 1319 | ListVarName={};% no new variable |
---|
| 1320 | VarAttribute={}; |
---|
[854] | 1321 | else |
---|
| 1322 | VarDimName=cell(size(ListVarName)); |
---|
| 1323 | for ilist=1:numel(ListVarName)% reshape data, excluding coordinates (ilist=1-2), TODO: rationalise |
---|
| 1324 | ListVarName{ilist}=regexprep(ListVarName{ilist},'(.+',''); |
---|
| 1325 | if ~isempty(find(strcmp(ListVarName{ilist},ProjData.ListVarName))) |
---|
| 1326 | ListVarName{ilist}=[ListVarName{ilist} '_1']; |
---|
| 1327 | end |
---|
| 1328 | ProjData.(ListVarName{ilist})=VarVal{ilist}; |
---|
| 1329 | VarDimName{ilist}={'coord_y','coord_x'}; |
---|
| 1330 | end |
---|
| 1331 | end |
---|
[869] | 1332 | if isfield (CellInfo{icell},'VarIndex_vector_x')&& isfield (CellInfo{icell},'VarIndex_vector_y') |
---|
[866] | 1333 | vector_x_proj=CellInfo{icell}.VarIndex_vector_x; %preserve for next cell |
---|
| 1334 | vector_y_proj=CellInfo{icell}.VarIndex_vector_y; %preserve for next cell |
---|
[869] | 1335 | end |
---|
[854] | 1336 | end |
---|
| 1337 | |
---|
| 1338 | case 'tps' |
---|
| 1339 | %% case of tps data (applies only in interp_tps mode) |
---|
| 1340 | if strcmp(ProjMode{icell},'interp_tps') |
---|
| 1341 | Coord=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex}); |
---|
[935] | 1342 | |
---|
| 1343 | |
---|
| 1344 | coord_x=Coord(:,1,:);% initial x coordinates |
---|
| 1345 | coord_y=Coord(:,2,:);% initial y coordinates |
---|
| 1346 | check3D=(numel(CellInfo{icell}.CoordIndex)==3); |
---|
| 1347 | if check3D |
---|
| 1348 | coord_z=Coord(:,3,:); |
---|
| 1349 | end |
---|
| 1350 | |
---|
| 1351 | % translate initial coordinates to account for the new origin |
---|
| 1352 | coord_x=coord_x-ObjectData.Coord(1,1); |
---|
| 1353 | coord_y=coord_y-ObjectData.Coord(1,2); |
---|
| 1354 | if check3D |
---|
| 1355 | coord_z=coord_z-ObjectData.Coord(1,3); |
---|
| 1356 | end |
---|
| 1357 | |
---|
| 1358 | % selection of the vectors in the projection range (3D case) |
---|
| 1359 | if check3D && width > 0 |
---|
| 1360 | %components of the unitiy vector normal to the projection plane |
---|
| 1361 | fieldZ=norm_plane(1)*coord_x + norm_plane(2)*coord_y+ norm_plane(3)*coord_z;% distance to the plane |
---|
| 1362 | indcut=find(abs(fieldZ) <= width); |
---|
| 1363 | for ivar=VarIndex |
---|
| 1364 | VarName=FieldData.ListVarName{ivar}; |
---|
| 1365 | FieldData.(VarName)=FieldData.(VarName)(indcut); |
---|
| 1366 | end |
---|
| 1367 | coord_x=coord_x(indcut); |
---|
| 1368 | coord_y=coord_y(indcut); |
---|
| 1369 | coord_z=coord_z(indcut); |
---|
| 1370 | end |
---|
| 1371 | |
---|
| 1372 | %rotate coordinates if needed: coord_X,coord_Y= = coordinates in the new plane |
---|
| 1373 | Psi=PlaneAngle(1); |
---|
| 1374 | Theta=PlaneAngle(2); |
---|
| 1375 | Phi=PlaneAngle(3); |
---|
| 1376 | if testangle && ~test90y && ~test90x;%=1 for slanted plane |
---|
| 1377 | coord_X=(coord_x *cos(Phi) + coord_y* sin(Phi)); |
---|
| 1378 | coord_Y=(-coord_x *sin(Phi) + coord_y *cos(Phi))*cos(Theta); |
---|
| 1379 | coord_Y=coord_Y+coord_z *sin(Theta); |
---|
| 1380 | coord_X=(coord_X *cos(Psi) - coord_Y* sin(Psi));%A VERIFIER |
---|
| 1381 | coord_Y=(coord_X *sin(Psi) + coord_Y* cos(Psi)); |
---|
| 1382 | else |
---|
| 1383 | coord_X=coord_x; |
---|
| 1384 | coord_Y=coord_y; |
---|
| 1385 | end |
---|
| 1386 | |
---|
| 1387 | %restriction to the range of X and Y if imposed by the projection object |
---|
| 1388 | |
---|
| 1389 | |
---|
| 1390 | Coord(:,1,:)=coord_X;% initial x coordinates |
---|
| 1391 | Coord(:,2,:)=coord_Y;% initial y coordinates |
---|
| 1392 | |
---|
| 1393 | |
---|
[854] | 1394 | NbCentres=FieldData.(FieldData.ListVarName{CellInfo{icell}.NbCentres_tps}); |
---|
| 1395 | SubRange=FieldData.(FieldData.ListVarName{CellInfo{icell}.SubRange_tps}); |
---|
| 1396 | if isfield(CellInfo{icell},'VarIndex_vector_x_tps')&&isfield(CellInfo{icell},'VarIndex_vector_y_tps') |
---|
| 1397 | FieldVar=cat(3,FieldData.(FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_x_tps}),FieldData.(FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_y_tps})); |
---|
| 1398 | end |
---|
[888] | 1399 | % interpolate data using thin plate spline |
---|
[854] | 1400 | [DataOut,VarAttribute,errormsg]=calc_field_tps(Coord,NbCentres,SubRange,FieldVar,CellInfo{icell}.FieldName,cat(3,XI,YI)); |
---|
[888] | 1401 | |
---|
[895] | 1402 | % set to NaN interpolation points which are too far from any initial data (more than 2 CoordMesh) |
---|
[905] | 1403 | Coord=permute(Coord,[1 3 2]); |
---|
| 1404 | Coord=reshape(Coord,size(Coord,1)*size(Coord,2),2); |
---|
| 1405 | if exist('scatteredInterpolant','file')%recent Matlab versions |
---|
| 1406 | F=scatteredInterpolant(Coord,Coord(:,1),'nearest'); |
---|
| 1407 | G=scatteredInterpolant(Coord,Coord(:,2),'nearest'); |
---|
| 1408 | else |
---|
| 1409 | F=TriScatteredInterp(Coord,Coord(:,1),'nearest'); |
---|
| 1410 | G=TriScatteredInterp(Coord,Coord(:,2),'nearest'); |
---|
| 1411 | end |
---|
| 1412 | Distx=F(XI,YI)-XI;% diff of x coordinates with the nearest measurement point |
---|
| 1413 | Disty=G(XI,YI)-YI;% diff of y coordinates with the nearest measurement point |
---|
| 1414 | Dist=Distx.*Distx+Disty.*Disty; |
---|
[854] | 1415 | ListVarName=(fieldnames(DataOut))'; |
---|
| 1416 | VarDimName=cell(size(ListVarName)); |
---|
| 1417 | for ilist=1:numel(ListVarName)% reshape data, excluding coordinates (ilist=1-2), TODO: rationalise |
---|
| 1418 | VarName=ListVarName{ilist}; |
---|
[912] | 1419 | VarDimName{ilist}={'coord_y','coord_x'}; |
---|
[854] | 1420 | ProjData.(VarName)=DataOut.(VarName); |
---|
[912] | 1421 | if ~isempty(thresh2) |
---|
| 1422 | ProjData.(VarName)(Dist>thresh2)=NaN;% put to NaN interpolated positions further than RangeInterp from initial data |
---|
| 1423 | end |
---|
[854] | 1424 | end |
---|
| 1425 | end |
---|
| 1426 | |
---|
| 1427 | case 'grid' |
---|
| 1428 | %% case of input fields defined on a structured grid |
---|
| 1429 | VarName=FieldData.ListVarName{VarIndex(1)};%get the first variable of the cell to get the input matrix dimensions |
---|
| 1430 | DimValue=size(FieldData.(VarName));%input matrix dimensions |
---|
| 1431 | DimValue(DimValue==1)=[];%remove singleton dimensions |
---|
| 1432 | NbDim=numel(DimValue);%update number of space dimensions |
---|
| 1433 | nbcolor=1; %default number of 'color' components: third matrix index without corresponding coordinate |
---|
| 1434 | if NbDim>=3 |
---|
| 1435 | if NbDim>3 |
---|
| 1436 | errormsg='matrices with more than 3 dimensions not handled'; |
---|
| 1437 | return |
---|
| 1438 | else |
---|
| 1439 | if numel(CellInfo{icell}.CoordIndex)==2% the third matrix dimension does not correspond to a space coordinate |
---|
| 1440 | nbcolor=DimValue(3); |
---|
| 1441 | DimValue(3)=[]; %number of 'color' components updated |
---|
| 1442 | NbDim=2;% space dimension set to 2 |
---|
| 1443 | end |
---|
| 1444 | end |
---|
| 1445 | end |
---|
| 1446 | Coord_z=[]; |
---|
| 1447 | Coord_y=[]; |
---|
| 1448 | Coord_x=[]; |
---|
| 1449 | |
---|
| 1450 | if testangle |
---|
| 1451 | ProjMode{icell}='interp_lin'; %request linear interpolation for projection on a tilted plane |
---|
| 1452 | end |
---|
| 1453 | |
---|
| 1454 | if isequal(ProjMode{icell},'projection')% && (~testangle || test90y || test90x) |
---|
| 1455 | if NbDim==2 && ~testXMin && ~testXMax && ~testYMin && ~testYMax% no range restriction |
---|
| 1456 | ListVarName=[ListVarName FieldData.ListVarName(VarIndex)]; |
---|
| 1457 | VarDimName=[VarDimName FieldData.VarDimName(VarIndex)]; |
---|
| 1458 | if isfield(FieldData,'VarAttribute') |
---|
| 1459 | VarAttribute=[VarAttribute FieldData.VarAttribute(VarIndex)]; |
---|
| 1460 | end |
---|
| 1461 | ProjData.(AYName)=FieldData.(AYName); |
---|
| 1462 | ProjData.(AXName)=FieldData.(AXName); |
---|
| 1463 | for ivar=VarIndex |
---|
| 1464 | VarName=FieldData.ListVarName{ivar}; |
---|
| 1465 | ProjData.(VarName)=FieldData.(VarName);% no change by projection |
---|
| 1466 | end |
---|
| 1467 | else |
---|
| 1468 | Coord{1}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(1)}); |
---|
| 1469 | Coord{2}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(2)}); |
---|
| 1470 | if NbDim==3 |
---|
| 1471 | Coord{3}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(3)}); |
---|
| 1472 | end |
---|
| 1473 | if numel(Coord{NbDim-1})==2 |
---|
| 1474 | DY=(Coord{NbDim-1}(2)-Coord{NbDim-1}(1))/(DimValue(1)-1); |
---|
| 1475 | end |
---|
| 1476 | if numel(Coord{NbDim})==2 |
---|
| 1477 | DX=(Coord{NbDim}(2)-Coord{NbDim}(1))/(DimValue(2)-1); |
---|
| 1478 | end |
---|
| 1479 | if testYMax |
---|
| 1480 | YIndexMax=(YMax-Coord{NbDim-1}(1))/DY+1;% matrix index corresponding to the max y value for the new field |
---|
| 1481 | if testYMin%test_direct(indY) |
---|
| 1482 | YIndexMin=(YMin-Coord{NbDim-1}(1))/DY+1;% matrix index corresponding to the min x value for the new field |
---|
| 1483 | else |
---|
| 1484 | YIndexMin=1; |
---|
| 1485 | end |
---|
| 1486 | else |
---|
| 1487 | YIndexMax=Coord{NbDim-1}(end)/DY; |
---|
| 1488 | YIndexMin=1; |
---|
| 1489 | end |
---|
| 1490 | if testXMax |
---|
| 1491 | XIndexMax=(XMax-Coord{NbDim}(1))/DY+1;% matrix index corresponding to the max y value for the new field |
---|
| 1492 | if testYMin%test_direct(indY) |
---|
| 1493 | XIndexMin=(XMin-Coord{NbDim}(1))/DX+1;% matrix index corresponding to the min x value for the new field |
---|
| 1494 | else |
---|
| 1495 | XIndexMin=1; |
---|
| 1496 | end |
---|
| 1497 | else |
---|
| 1498 | XIndexMax=Coord{NbDim}(end)/DX; |
---|
| 1499 | XIndexMin=1; |
---|
| 1500 | end |
---|
| 1501 | YIndexRange(1)=ceil(min(YIndexMin,YIndexMax));%first y index to select from the previous field |
---|
| 1502 | YIndexRange(1)=max(YIndexRange(1),1);% avoid bound lower than the first index |
---|
| 1503 | YIndexRange(2)=floor(max(YIndexMin,YIndexMax));%last y index to select from the previous field |
---|
| 1504 | YIndexRange(2)=min(YIndexRange(2),DimValue(NbDim-1));% limit to the last available index |
---|
| 1505 | XIndexRange(1)=ceil(min(XIndexMin,XIndexMax));%first x index to select from the previous field |
---|
| 1506 | XIndexRange(1)=max(XIndexRange(1),1);% avoid bound lower than the first index |
---|
| 1507 | XIndexRange(2)=floor(max(XIndexMin,XIndexMax));%last x index to select from the previous field |
---|
| 1508 | XIndexRange(2)=min(XIndexRange(2),DimValue(NbDim));% limit to the last available index |
---|
| 1509 | if test90y |
---|
| 1510 | ind_new=[3 2 1]; |
---|
| 1511 | DimCell={AYProjName,AXProjName}; |
---|
| 1512 | iz=ceil((ObjectData.Coord(1,1)-Coord{3}(1))/DX)+1; |
---|
| 1513 | for ivar=VarIndex |
---|
| 1514 | VarName=FieldData.ListVarName{ivar}; |
---|
| 1515 | ListVarName=[ListVarName VarName]; |
---|
| 1516 | VarDimName=[VarDimName {DimCell}]; |
---|
| 1517 | VarAttribute{length(ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes |
---|
| 1518 | ProjData.(VarName)=permute(FieldData.(VarName),ind_new);% permute x and z indices for 90 degree rotation |
---|
| 1519 | ProjData.(VarName)=squeeze(ProjData.(VarName)(iz,:,:));% select the z index iz |
---|
| 1520 | end |
---|
| 1521 | ProjData.(AYName)=[Ybound(1) Ybound(2)]; %record the new (projected ) y coordinates |
---|
| 1522 | ProjData.(AXName)=[Coord{1}(end),Coord{1}(1)]; %record the new (projected ) x coordinates |
---|
| 1523 | else |
---|
| 1524 | if NbDim==3 |
---|
| 1525 | DZ=(Coord{1}(end)-Coord{1}(1))/(numel(Coord{1})-1); |
---|
| 1526 | DimCell(1)=[]; %suppress z variable |
---|
| 1527 | DimValue(1)=[]; |
---|
| 1528 | test_direct=1;%TOdo; GENERALIZE, SEE CASE OF points |
---|
| 1529 | if test_direct(1) |
---|
| 1530 | iz=ceil((ObjectData.Coord(1,3)-Coord{1}(1))/DZ)+1; |
---|
| 1531 | else |
---|
| 1532 | iz=ceil((Coord{1}(1)-ObjectData.Coord(1,3))/DZ)+1; |
---|
| 1533 | end |
---|
| 1534 | end |
---|
| 1535 | for ivar=VarIndex% loop on non coordinate variables |
---|
| 1536 | VarName=FieldData.ListVarName{ivar}; |
---|
| 1537 | ListVarName=[ListVarName VarName]; |
---|
| 1538 | VarDimName=[VarDimName {DimCell}]; |
---|
| 1539 | if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute)>=ivar |
---|
| 1540 | VarAttribute{length(ListVarName)}=FieldData.VarAttribute{ivar}; |
---|
| 1541 | end |
---|
| 1542 | if NbDim==3 |
---|
| 1543 | ProjData.(VarName)=squeeze(FieldData.(VarName)(iz,YIndexRange(1):YIndexRange(end),XIndexRange(1):XIndexRange(end))); |
---|
| 1544 | else |
---|
| 1545 | ProjData.(VarName)=FieldData.(VarName)(YIndexRange(1):YIndexRange(end),XIndexRange(1):XIndexRange(end),:); |
---|
| 1546 | end |
---|
| 1547 | end |
---|
| 1548 | if testXMax |
---|
| 1549 | ProjData.(AXName)=Coord{NbDim}(1)+DX*(XIndexRange-1); %record the new (projected ) x coordinates |
---|
| 1550 | else |
---|
| 1551 | ProjData.(AXName)=FieldData.(AXName); |
---|
| 1552 | end |
---|
| 1553 | if testYMax |
---|
| 1554 | ProjData.(AYName)=Coord{NbDim-1}(1)+DY*(YIndexRange-1); %record the new (projected ) x coordinates |
---|
| 1555 | else |
---|
| 1556 | ProjData.(AYName)=FieldData.(AYName); |
---|
| 1557 | end |
---|
| 1558 | end |
---|
| 1559 | end |
---|
| 1560 | else % case with interpolation on a grid |
---|
| 1561 | if NbDim==2 %2D case |
---|
| 1562 | if isequal(ProjMode{icell},'interp_tps') |
---|
| 1563 | npx_interp_tps=ceil(abs(DX/DAX)); |
---|
| 1564 | npy_interp_tps=ceil(abs(DY/DAY)); |
---|
| 1565 | Minterp_tps=ones(npy_interp_tps,npx_interp_tps)/(npx_interp_tps*npy_interp_tps); |
---|
| 1566 | test_interp_tps=1; |
---|
| 1567 | else |
---|
| 1568 | test_interp_tps=0; |
---|
| 1569 | end |
---|
| 1570 | Coord{1}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(1)}); |
---|
| 1571 | Coord{2}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(2)}); |
---|
| 1572 | if ~(testXMin && testYMin)% % if the range of the projected coordinates is not fully defined by the projection object, find the extrema of the projected field |
---|
| 1573 | xcorner=[min(Coord{NbDim}) max(Coord{NbDim}) max(Coord{NbDim}) min(Coord{NbDim})]-ObjectData.Coord(1,1);% corner absissa of the original grid with respect to the new origin |
---|
| 1574 | ycorner=[min(Coord{NbDim-1}) min(Coord{NbDim-1}) max(Coord{NbDim-1}) max(Coord{NbDim-1})]-ObjectData.Coord(1,2);% corner ordinates of the original grid |
---|
| 1575 | xcor_new=xcorner*cos(PlaneAngle(3))+ycorner*sin(PlaneAngle(3));%coordinates of the corners in new frame |
---|
| 1576 | ycor_new=-xcorner*sin(PlaneAngle(3))+ycorner*cos(PlaneAngle(3)); |
---|
| 1577 | if ~testXMin |
---|
| 1578 | XMin=min(xcor_new); |
---|
| 1579 | end |
---|
| 1580 | if ~testXMax |
---|
| 1581 | XMax=max(xcor_new); |
---|
| 1582 | end |
---|
| 1583 | if ~testYMin |
---|
| 1584 | YMin=min(ycor_new); |
---|
| 1585 | end |
---|
| 1586 | if ~testYMax |
---|
| 1587 | YMax=max(ycor_new); |
---|
| 1588 | end |
---|
| 1589 | end |
---|
| 1590 | coord_x_proj=XMin:DX:XMax; |
---|
| 1591 | coord_y_proj=YMin:DY:YMax; |
---|
| 1592 | ProjData.(AYName)=[coord_y_proj(1) coord_y_proj(end)]; %record the new (projected ) y coordinates |
---|
| 1593 | ProjData.(AXName)=[coord_x_proj(1) coord_x_proj(end)]; %record the new (projected ) x coordinates |
---|
| 1594 | [X,YI]=meshgrid(coord_x_proj,coord_y_proj);%grid in the new coordinates |
---|
| 1595 | XI=ObjectData.Coord(1,1)+(X)*cos(PlaneAngle(3))-YI*sin(PlaneAngle(3));%corresponding coordinates in the original system |
---|
| 1596 | YI=ObjectData.Coord(1,2)+(X)*sin(PlaneAngle(3))+YI*cos(PlaneAngle(3)); |
---|
| 1597 | if numel(Coord{1})==2% x coordiante defiend by its bounds, get the whole set |
---|
| 1598 | Coord{1}=linspace(Coord{1}(1),Coord{1}(2),CellInfo{icell}.CoordSize(1)); |
---|
| 1599 | end |
---|
| 1600 | if numel(Coord{2})==2% y coordiante defiend by its bounds, get the whole set |
---|
| 1601 | Coord{2}=linspace(Coord{2}(1),Coord{2}(2),CellInfo{icell}.CoordSize(2)); |
---|
| 1602 | end |
---|
| 1603 | [X,Y]=meshgrid(Coord{2},Coord{1});%initial coordinates |
---|
| 1604 | %name of error flag variable |
---|
| 1605 | FFName='FF';%default name (if not already used) |
---|
| 1606 | if isfield(ProjData,'FF') |
---|
| 1607 | ind=1; |
---|
| 1608 | while isfield(ProjData,['FF_' num2str(ind)]) |
---|
| 1609 | ind=ind+1; |
---|
| 1610 | end |
---|
| 1611 | FFName=['FF_' num2str(ind)];% append an index to the name of error flag, FF_1,FF_2... |
---|
| 1612 | end |
---|
| 1613 | % project all variables in the cell |
---|
| 1614 | for ivar=VarIndex |
---|
| 1615 | VarName=FieldData.ListVarName{ivar}; |
---|
| 1616 | if size(FieldData.(VarName),3)==1 |
---|
| 1617 | ProjData.(VarName)=interp2(X,Y,double(FieldData.(VarName)),XI,YI,'*linear');%interpolation fct |
---|
| 1618 | else |
---|
| 1619 | ProjData.(VarName)=interp2(X,Y,double(FieldData.(VarName)(:,:,1)),XI,YI,'*linear'); |
---|
| 1620 | for icolor=2:size(FieldData.(VarName),3)% project 'color' components |
---|
| 1621 | ProjData.(VarName)=cat(3,ProjData.(VarName),interp2(X,Y,double(FieldData.(VarName)(:,:,icolor)),XI,YI,'*linear')); |
---|
| 1622 | end |
---|
| 1623 | end |
---|
| 1624 | if isa(FieldData.(VarName),'uint8') |
---|
| 1625 | ProjData.(VarName)=uint8(ProjData.(VarName));%put result to integer 8 bits if the initial field is integer (image) |
---|
| 1626 | elseif isa(FieldData.(VarName),'uint16') |
---|
| 1627 | ProjData.(VarName)=uint16(ProjData.(VarName));%put result to integer 16 bits if the initial field is integer (image) |
---|
| 1628 | end |
---|
| 1629 | ListVarName=[ListVarName VarName]; |
---|
| 1630 | DimCell(1:2)={AYName,AXName}; |
---|
| 1631 | VarDimName=[VarDimName {DimCell}]; |
---|
| 1632 | if isfield(FieldData,'VarAttribute')&&length(FieldData.VarAttribute)>=ivar |
---|
| 1633 | VarAttribute{length(ListVarName)+nbcoord}=FieldData.VarAttribute{ivar}; |
---|
| 1634 | end; |
---|
| 1635 | ProjData.(FFName)=isnan(ProjData.(VarName));%detact NaN (points outside the interpolation range) |
---|
| 1636 | ProjData.(VarName)(ProjData.(FFName))=0; %set to 0 the NaN data |
---|
| 1637 | end |
---|
| 1638 | %update list of variables with error flag |
---|
| 1639 | ListVarName=[ListVarName FFName]; |
---|
| 1640 | VarDimName=[VarDimName {DimCell}]; |
---|
| 1641 | VarAttribute{numel(ListVarName)}.Role='errorflag'; |
---|
| 1642 | elseif ~testangle |
---|
| 1643 | % unstructured z coordinate |
---|
| 1644 | test_sup=(Coord{1}>=ObjectData.Coord(1,3)); |
---|
| 1645 | iz_sup=find(test_sup); |
---|
| 1646 | iz=iz_sup(1); |
---|
| 1647 | if iz>=1 & iz<=npz |
---|
| 1648 | %ProjData.ListDimName=[ProjData.ListDimName ListDimName(2:end)]; |
---|
| 1649 | %ProjData.DimValue=[ProjData.DimValue npY npX]; |
---|
| 1650 | for ivar=VarIndex |
---|
| 1651 | VarName=FieldData.ListVarName{ivar}; |
---|
| 1652 | ListVarName=[ListVarName VarName]; |
---|
| 1653 | VarAttribute{length(ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes |
---|
[866] | 1654 | ProjData.(VarName)=squeeze(FieldData.(VarName)(iz,:,:));% select the z index iz |
---|
[854] | 1655 | %TODO : do a vertical average for a thick plane |
---|
| 1656 | if test_interp(2) || test_interp(3) |
---|
[866] | 1657 | ProjData.(VarName)=interp2(Coord{3},Coord{2},ProjData.(VarName),Coord_x,Coord_y); |
---|
[854] | 1658 | end |
---|
| 1659 | end |
---|
| 1660 | end |
---|
| 1661 | else |
---|
| 1662 | errormsg='projection of structured coordinates on oblique plane not yet implemented'; |
---|
| 1663 | %TODO: use interp_lin3 |
---|
| 1664 | return |
---|
| 1665 | end |
---|
| 1666 | end |
---|
| 1667 | end |
---|
| 1668 | % update the global list of projected variables: |
---|
| 1669 | ProjData.ListVarName=[ProjData.ListVarName ListVarName]; |
---|
| 1670 | ProjData.VarDimName=[ProjData.VarDimName VarDimName]; |
---|
| 1671 | ProjData.VarAttribute=[ProjData.VarAttribute VarAttribute]; |
---|
| 1672 | |
---|
| 1673 | %% projection of velocity components in the rotated coordinates |
---|
| 1674 | if testangle |
---|
| 1675 | ivar_U=[];ivar_V=[];ivar_W=[]; |
---|
| 1676 | for ivar=1:numel(VarAttribute) |
---|
| 1677 | if isfield(VarAttribute{ivar},'Role') |
---|
| 1678 | if strcmp(VarAttribute{ivar}.Role,'vector_x') |
---|
| 1679 | ivar_U=ivar; |
---|
| 1680 | elseif strcmp(VarAttribute{ivar}.Role,'vector_y') |
---|
| 1681 | ivar_V=ivar; |
---|
| 1682 | elseif strcmp(VarAttribute{ivar}.Role,'vector_z') |
---|
| 1683 | ivar_W=ivar; |
---|
| 1684 | end |
---|
| 1685 | end |
---|
| 1686 | end |
---|
| 1687 | if ~isempty(ivar_U) |
---|
| 1688 | if isempty(ivar_V) |
---|
| 1689 | msgbox_uvmat('ERROR','v velocity component missing in proj_field.m') |
---|
| 1690 | return |
---|
| 1691 | else |
---|
| 1692 | UName=ListVarName{ivar_U}; |
---|
| 1693 | VName=ListVarName{ivar_V}; |
---|
[933] | 1694 | UValue=cos(PlaneAngle(3))*ProjData.(UName)+ sin(PlaneAngle(3))*ProjData.(VName); |
---|
[854] | 1695 | ProjData.(VName)=(-sin(PlaneAngle(3))*ProjData.(UName)+ cos(PlaneAngle(3))*ProjData.(VName)); |
---|
[933] | 1696 | ProjData.(UName)=UValue; |
---|
[854] | 1697 | if ~isempty(ivar_W) |
---|
| 1698 | WName=FieldData.ListVarName{ivar_W}; |
---|
[933] | 1699 | VValue=ProjData.(VName)+ ProjData.(WName)*sin(Theta);% |
---|
[866] | 1700 | ProjData.(WName)=NormVec_X*ProjData.(UName)+ NormVec_Y*ProjData.(VName)+ NormVec_Z* ProjData.(WName); |
---|
[933] | 1701 | ProjData.(VName)=VValue; |
---|
[854] | 1702 | end |
---|
| 1703 | end |
---|
| 1704 | end |
---|
| 1705 | end |
---|
| 1706 | end |
---|
| 1707 | % %prepare substraction in case of two input fields |
---|
| 1708 | % SubData.ListVarName={}; |
---|
| 1709 | % SubData.VarDimName={}; |
---|
| 1710 | % SubData.VarAttribute={}; |
---|
| 1711 | % check_remove=zeros(size(ProjData.ListVarName)); |
---|
| 1712 | % for iproj=1:numel(ProjData.VarAttribute) |
---|
| 1713 | % if isfield(ProjData.VarAttribute{iproj},'CheckSub')&&isequal(ProjData.VarAttribute{iproj}.CheckSub,1) |
---|
| 1714 | % VarName=ProjData.ListVarName{iproj}; |
---|
| 1715 | % SubData.ListVarName=[SubData.ListVarName {VarName}]; |
---|
| 1716 | % SubData.VarDimName=[SubData.VarDimName ProjData.VarDimName{iproj}]; |
---|
| 1717 | % SubData.VarAttribute=[SubData.VarAttribute ProjData.VarAttribute{iproj}]; |
---|
| 1718 | % SubData.(VarName)=ProjData.(VarName); |
---|
| 1719 | % check_remove(iproj)=1; |
---|
| 1720 | % end |
---|
| 1721 | % end |
---|
| 1722 | % if ~isempty(find(check_remove)) |
---|
| 1723 | % ind_remove=find(check_remove); |
---|
| 1724 | % ProjData.ListVarName(ind_remove)=[]; |
---|
| 1725 | % ProjData.VarDimName(ind_remove)=[]; |
---|
| 1726 | % ProjData.VarAttribute(ind_remove)=[]; |
---|
| 1727 | % ProjData=sub_field(ProjData,[],SubData); |
---|
| 1728 | % end |
---|
| 1729 | |
---|
| 1730 | %----------------------------------------------------------------- |
---|
| 1731 | %projection in a volume |
---|
| 1732 | function [ProjData,errormsg] = proj_volume(FieldData, ObjectData) |
---|
| 1733 | |
---|
| 1734 | %----------------------------------------------------------------- |
---|
| 1735 | ProjData=FieldData;%default output |
---|
| 1736 | |
---|
| 1737 | %% axis origin |
---|
| 1738 | if isempty(ObjectData.Coord) |
---|
| 1739 | ObjectData.Coord(1,1)=0;%origin of the plane set to [0 0] by default |
---|
| 1740 | ObjectData.Coord(1,2)=0; |
---|
| 1741 | ObjectData.Coord(1,3)=0; |
---|
| 1742 | end |
---|
| 1743 | |
---|
| 1744 | %% rotation angles |
---|
| 1745 | VolumeAngle=[0 0 0]; |
---|
| 1746 | norm_plane=[0 0 1]; |
---|
| 1747 | if isfield(ObjectData,'Angle')&& isequal(size(ObjectData.Angle),[1 3])&& ~isequal(ObjectData.Angle,[0 0 0]) |
---|
| 1748 | PlaneAngle=ObjectData.Angle; |
---|
| 1749 | VolumeAngle=ObjectData.Angle; |
---|
| 1750 | om=norm(VolumeAngle);%norm of rotation angle in radians |
---|
| 1751 | OmAxis=VolumeAngle/om; %unit vector marking the rotation axis |
---|
| 1752 | cos_om=cos(pi*om/180); |
---|
| 1753 | sin_om=sin(pi*om/180); |
---|
| 1754 | coeff=OmAxis(3)*(1-cos_om); |
---|
| 1755 | %components of the unity vector norm_plane normal to the projection plane |
---|
| 1756 | norm_plane(1)=OmAxis(1)*coeff+OmAxis(2)*sin_om; |
---|
| 1757 | norm_plane(2)=OmAxis(2)*coeff-OmAxis(1)*sin_om; |
---|
| 1758 | norm_plane(3)=OmAxis(3)*coeff+cos_om; |
---|
| 1759 | end |
---|
| 1760 | testangle=~isequal(VolumeAngle,[0 0 0]); |
---|
| 1761 | |
---|
| 1762 | %% mesh sizes DX, DY, DZ |
---|
| 1763 | DX=0; |
---|
| 1764 | DY=0; %default |
---|
| 1765 | DZ=0; |
---|
| 1766 | if isfield(ObjectData,'DX')&~isempty(ObjectData.DX) |
---|
| 1767 | DX=abs(ObjectData.DX);%mesh of interpolation points |
---|
| 1768 | end |
---|
| 1769 | if isfield(ObjectData,'DY')&~isempty(ObjectData.DY) |
---|
| 1770 | DY=abs(ObjectData.DY);%mesh of interpolation points |
---|
| 1771 | end |
---|
| 1772 | if isfield(ObjectData,'DZ')&~isempty(ObjectData.DZ) |
---|
| 1773 | DZ=abs(ObjectData.DZ);%mesh of interpolation points |
---|
| 1774 | end |
---|
| 1775 | if ~strcmp(ProjMode,'projection') && (DX==0||DY==0||DZ==0) |
---|
| 1776 | errormsg='grid mesh DX , DY or DZ is missing'; |
---|
| 1777 | return |
---|
| 1778 | end |
---|
| 1779 | |
---|
| 1780 | %% extrema along each axis |
---|
| 1781 | testXMin=0; |
---|
| 1782 | testXMax=0; |
---|
| 1783 | testYMin=0; |
---|
| 1784 | testYMax=0; |
---|
| 1785 | if isfield(ObjectData,'RangeX') |
---|
| 1786 | XMin=min(ObjectData.RangeX); |
---|
| 1787 | XMax=max(ObjectData.RangeX); |
---|
| 1788 | testXMin=XMax>XMin; |
---|
| 1789 | testXMax=1; |
---|
| 1790 | end |
---|
| 1791 | if isfield(ObjectData,'RangeY') |
---|
| 1792 | YMin=min(ObjectData.RangeY); |
---|
| 1793 | YMax=max(ObjectData.RangeY); |
---|
| 1794 | testYMin=YMax>YMin; |
---|
| 1795 | testYMax=1; |
---|
| 1796 | end |
---|
| 1797 | width=0;%default width of the projection band |
---|
| 1798 | if isfield(ObjectData,'RangeZ') |
---|
| 1799 | ZMin=min(ObjectData.RangeZ); |
---|
| 1800 | ZMax=max(ObjectData.RangeZ); |
---|
| 1801 | testZMin=ZMax>ZMin; |
---|
| 1802 | testZMax=1; |
---|
| 1803 | end |
---|
| 1804 | |
---|
| 1805 | %% initiate Matlab structure for physical field |
---|
| 1806 | [ProjData,errormsg]=proj_heading(FieldData,ObjectData); |
---|
| 1807 | if ~isempty(errormsg) |
---|
| 1808 | return |
---|
| 1809 | end |
---|
| 1810 | |
---|
| 1811 | ProjData.NbDim=3; |
---|
| 1812 | ProjData.ListVarName={}; |
---|
| 1813 | ProjData.VarDimName={}; |
---|
| 1814 | if ~isequal(DX,0)&& ~isequal(DY,0) |
---|
| 1815 | ProjData.CoordMesh=sqrt(DX*DY);%define typical data mesh, useful for mouse selection in plots |
---|
| 1816 | elseif isfield(FieldData,'CoordMesh') |
---|
| 1817 | ProjData.CoordMesh=FieldData.CoordMesh; |
---|
| 1818 | end |
---|
| 1819 | |
---|
| 1820 | error=0;%default |
---|
| 1821 | flux=0; |
---|
| 1822 | testfalse=0; |
---|
| 1823 | ListIndex={}; |
---|
| 1824 | |
---|
| 1825 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 1826 | %% group the variables (fields of 'FieldData') in cells of variables with the same dimensions |
---|
| 1827 | %----------------------------------------------------------------- |
---|
| 1828 | idimvar=0; |
---|
| 1829 | % LOOP ON GROUPS OF VARIABLES SHARING THE SAME DIMENSIONS |
---|
| 1830 | % CellVarIndex=cells of variable index arrays |
---|
| 1831 | ivar_new=0; % index of the current variable in the projected field |
---|
| 1832 | icoord=0; |
---|
| 1833 | nbcoord=0;%number of added coordinate variables brought by projection |
---|
| 1834 | nbvar=0; |
---|
| 1835 | for icell=1:length(CellVarIndex) |
---|
| 1836 | NbDim=NbDimVec(icell); |
---|
| 1837 | if NbDim<3 |
---|
| 1838 | continue |
---|
| 1839 | end |
---|
| 1840 | VarIndex=CellVarIndex{icell};% indices of the selected variables in the list FieldData.ListVarName |
---|
| 1841 | VarType=VarTypeCell{icell}; |
---|
| 1842 | ivar_X=VarType.coord_x; |
---|
| 1843 | ivar_Y=VarType.coord_y; |
---|
| 1844 | ivar_Z=VarType.coord_z; |
---|
| 1845 | ivar_U=VarType.vector_x; |
---|
| 1846 | ivar_V=VarType.vector_y; |
---|
| 1847 | ivar_W=VarType.vector_z; |
---|
| 1848 | ivar_C=VarType.scalar ; |
---|
| 1849 | ivar_Anc=VarType.ancillary; |
---|
| 1850 | test_anc=zeros(size(VarIndex)); |
---|
| 1851 | test_anc(ivar_Anc)=ones(size(ivar_Anc)); |
---|
| 1852 | ivar_F=VarType.warnflag; |
---|
| 1853 | ivar_FF=VarType.errorflag; |
---|
| 1854 | check_unstructured_coord=~isempty(ivar_X) && ~isempty(ivar_Y); |
---|
| 1855 | DimCell=FieldData.VarDimName{VarIndex(1)}; |
---|
| 1856 | if ischar(DimCell) |
---|
| 1857 | DimCell={DimCell};%name of dimensions |
---|
| 1858 | end |
---|
| 1859 | |
---|
| 1860 | %% case of input fields with unstructured coordinates |
---|
| 1861 | if check_unstructured_coord |
---|
| 1862 | XName=FieldData.ListVarName{ivar_X}; |
---|
| 1863 | YName=FieldData.ListVarName{ivar_Y}; |
---|
| 1864 | eval(['coord_x=FieldData.' XName ';']) |
---|
| 1865 | eval(['coord_y=FieldData.' YName ';']) |
---|
| 1866 | if length(ivar_Z)==1 |
---|
| 1867 | ZName=FieldData.ListVarName{ivar_Z}; |
---|
| 1868 | eval(['coord_z=FieldData.' ZName ';']) |
---|
| 1869 | end |
---|
| 1870 | |
---|
| 1871 | % translate initial coordinates |
---|
| 1872 | coord_x=coord_x-ObjectData.Coord(1,1); |
---|
| 1873 | coord_y=coord_y-ObjectData.Coord(1,2); |
---|
| 1874 | if ~isempty(ivar_Z) |
---|
| 1875 | coord_z=coord_z-ObjectData.Coord(1,3); |
---|
| 1876 | end |
---|
| 1877 | |
---|
| 1878 | % selection of the vectors in the projection range |
---|
| 1879 | % if length(ivar_Z)==1 && width > 0 |
---|
| 1880 | % %components of the unitiy vector normal to the projection plane |
---|
| 1881 | % fieldZ=NormVec_X*coord_x + NormVec_Y*coord_y+ NormVec_Z*coord_z;% distance to the plane |
---|
| 1882 | % indcut=find(abs(fieldZ) <= width); |
---|
| 1883 | % for ivar=VarIndex |
---|
| 1884 | % VarName=FieldData.ListVarName{ivar}; |
---|
| 1885 | % eval(['FieldData.' VarName '=FieldData.' VarName '(indcut);']) |
---|
| 1886 | % % A VOIR : CAS DE VAR STRUCTUREE MAIS PAS GRILLE REGULIERE : INTERPOLER SUR GRILLE REGULIERE |
---|
| 1887 | % end |
---|
| 1888 | % coord_x=coord_x(indcut); |
---|
| 1889 | % coord_y=coord_y(indcut); |
---|
| 1890 | % coord_z=coord_z(indcut); |
---|
| 1891 | % end |
---|
| 1892 | |
---|
| 1893 | %rotate coordinates if needed: TODO modify |
---|
| 1894 | if testangle |
---|
| 1895 | coord_X=(coord_x *cos(Phi) + coord_y* sin(Phi)); |
---|
| 1896 | coord_Y=(-coord_x *sin(Phi) + coord_y *cos(Phi))*cos(Theta); |
---|
| 1897 | if ~isempty(ivar_Z) |
---|
| 1898 | coord_Y=coord_Y+coord_z *sin(Theta); |
---|
| 1899 | end |
---|
| 1900 | |
---|
| 1901 | coord_X=(coord_X *cos(Psi) - coord_Y* sin(Psi));%A VERIFIER |
---|
| 1902 | coord_Y=(coord_X *sin(Psi) + coord_Y* cos(Psi)); |
---|
| 1903 | |
---|
| 1904 | else |
---|
| 1905 | coord_X=coord_x; |
---|
| 1906 | coord_Y=coord_y; |
---|
| 1907 | coord_Z=coord_z; |
---|
| 1908 | end |
---|
| 1909 | %restriction to the range of x and y if imposed |
---|
| 1910 | testin=ones(size(coord_X)); %default |
---|
| 1911 | testbound=0; |
---|
| 1912 | if testXMin |
---|
| 1913 | testin=testin & (coord_X >= XMin); |
---|
| 1914 | testbound=1; |
---|
| 1915 | end |
---|
| 1916 | if testXMax |
---|
| 1917 | testin=testin & (coord_X <= XMax); |
---|
| 1918 | testbound=1; |
---|
| 1919 | end |
---|
| 1920 | if testYMin |
---|
| 1921 | testin=testin & (coord_Y >= YMin); |
---|
| 1922 | testbound=1; |
---|
| 1923 | end |
---|
| 1924 | if testYMax |
---|
| 1925 | testin=testin & (coord_Y <= YMax); |
---|
| 1926 | testbound=1; |
---|
| 1927 | end |
---|
| 1928 | if testbound |
---|
| 1929 | indcut=find(testin); |
---|
| 1930 | for ivar=VarIndex |
---|
| 1931 | VarName=FieldData.ListVarName{ivar}; |
---|
| 1932 | eval(['FieldData.' VarName '=FieldData.' VarName '(indcut);']) |
---|
| 1933 | end |
---|
| 1934 | coord_X=coord_X(indcut); |
---|
| 1935 | coord_Y=coord_Y(indcut); |
---|
| 1936 | if length(ivar_Z)==1 |
---|
| 1937 | coord_Z=coord_Z(indcut); |
---|
| 1938 | end |
---|
| 1939 | end |
---|
| 1940 | % different cases of projection |
---|
| 1941 | if isequal(ObjectData.ProjMode,'projection')%%%%%%% NOT USED %%%%%%%%%% |
---|
| 1942 | for ivar=VarIndex %transfer variables to the projection plane |
---|
| 1943 | VarName=FieldData.ListVarName{ivar}; |
---|
| 1944 | if ivar==ivar_X %x coordinate |
---|
| 1945 | eval(['ProjData.' VarName '=coord_X;']) |
---|
| 1946 | elseif ivar==ivar_Y % y coordinate |
---|
| 1947 | eval(['ProjData.' VarName '=coord_Y;']) |
---|
| 1948 | elseif isempty(ivar_Z) || ivar~=ivar_Z % other variables (except Z coordinate wyhich is not reproduced) |
---|
| 1949 | eval(['ProjData.' VarName '=FieldData.' VarName ';']) |
---|
| 1950 | end |
---|
| 1951 | if isempty(ivar_Z) || ivar~=ivar_Z |
---|
| 1952 | ProjData.ListVarName=[ProjData.ListVarName VarName]; |
---|
| 1953 | ProjData.VarDimName=[ProjData.VarDimName DimCell]; |
---|
| 1954 | nbvar=nbvar+1; |
---|
| 1955 | if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute) >=ivar |
---|
| 1956 | ProjData.VarAttribute{nbvar}=FieldData.VarAttribute{ivar}; |
---|
| 1957 | end |
---|
| 1958 | end |
---|
| 1959 | end |
---|
| 1960 | elseif isequal(ObjectData.ProjMode,'interp_lin')||isequal(ObjectData.ProjMode,'interp_tps')%interpolate data on a regular grid |
---|
| 1961 | coord_x_proj=XMin:DX:XMax; |
---|
| 1962 | coord_y_proj=YMin:DY:YMax; |
---|
| 1963 | coord_z_proj=ZMin:DZ:ZMax; |
---|
| 1964 | DimCell={'coord_z','coord_y','coord_x'}; |
---|
| 1965 | ProjData.ListVarName={'coord_z','coord_y','coord_x'}; |
---|
| 1966 | ProjData.VarDimName={'coord_z','coord_y','coord_x'}; |
---|
| 1967 | nbcoord=2; |
---|
| 1968 | ProjData.coord_z=[ZMin ZMax]; |
---|
| 1969 | ProjData.coord_y=[YMin YMax]; |
---|
| 1970 | ProjData.coord_x=[XMin XMax]; |
---|
| 1971 | if isempty(ivar_X), ivar_X=0; end; |
---|
| 1972 | if isempty(ivar_Y), ivar_Y=0; end; |
---|
| 1973 | if isempty(ivar_Z), ivar_Z=0; end; |
---|
| 1974 | if isempty(ivar_U), ivar_U=0; end; |
---|
| 1975 | if isempty(ivar_V), ivar_V=0; end; |
---|
| 1976 | if isempty(ivar_W), ivar_W=0; end; |
---|
| 1977 | if isempty(ivar_F), ivar_F=0; end; |
---|
| 1978 | if isempty(ivar_FF), ivar_FF=0; end; |
---|
| 1979 | if ~isequal(ivar_FF,0) |
---|
| 1980 | VarName_FF=FieldData.ListVarName{ivar_FF}; |
---|
| 1981 | eval(['indsel=find(FieldData.' VarName_FF '==0);']) |
---|
| 1982 | coord_X=coord_X(indsel); |
---|
| 1983 | coord_Y=coord_Y(indsel); |
---|
| 1984 | end |
---|
| 1985 | FF=zeros(1,length(coord_y_proj)*length(coord_x_proj)); |
---|
| 1986 | testFF=0; |
---|
| 1987 | [X,Y,Z]=meshgrid(coord_y_proj,coord_z_proj,coord_x_proj);%grid in the new coordinates |
---|
| 1988 | for ivar=VarIndex |
---|
| 1989 | VarName=FieldData.ListVarName{ivar}; |
---|
| 1990 | if ~( ivar==ivar_X || ivar==ivar_Y || ivar==ivar_Z || ivar==ivar_F || ivar==ivar_FF || test_anc(ivar)==1) |
---|
| 1991 | ivar_new=ivar_new+1; |
---|
| 1992 | ProjData.ListVarName=[ProjData.ListVarName {VarName}]; |
---|
| 1993 | ProjData.VarDimName=[ProjData.VarDimName {DimCell}]; |
---|
| 1994 | if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute) >=ivar |
---|
| 1995 | ProjData.VarAttribute{ivar_new+nbcoord}=FieldData.VarAttribute{ivar}; |
---|
| 1996 | end |
---|
| 1997 | if ~isequal(ivar_FF,0) |
---|
| 1998 | eval(['FieldData.' VarName '=FieldData.' VarName '(indsel);']) |
---|
| 1999 | end |
---|
| 2000 | % linear interpolation |
---|
| 2001 | InterpFct=TriScatteredInterp(double(coord_X),double(coord_Y),double(coord_Z),double(FieldData.(VarName))); |
---|
| 2002 | ProjData.(VarName)=InterpFct(X,Y,Z); |
---|
| 2003 | % eval(['varline=reshape(ProjData.' VarName ',1,length(coord_y_proj)*length(coord_x_proj));']) |
---|
| 2004 | % FFlag= isnan(varline); %detect undefined values NaN |
---|
| 2005 | % indnan=find(FFlag); |
---|
| 2006 | % if ~isempty(indnan) |
---|
| 2007 | % varline(indnan)=zeros(size(indnan)); |
---|
| 2008 | % eval(['ProjData.' VarName '=reshape(varline,length(coord_y_proj),length(coord_x_proj));']) |
---|
| 2009 | % FF(indnan)=ones(size(indnan)); |
---|
| 2010 | % testFF=1; |
---|
| 2011 | % end |
---|
| 2012 | if ivar==ivar_U |
---|
| 2013 | ivar_U=ivar_new; |
---|
| 2014 | end |
---|
| 2015 | if ivar==ivar_V |
---|
| 2016 | ivar_V=ivar_new; |
---|
| 2017 | end |
---|
| 2018 | if ivar==ivar_W |
---|
| 2019 | ivar_W=ivar_new; |
---|
| 2020 | end |
---|
| 2021 | end |
---|
| 2022 | end |
---|
| 2023 | if testFF |
---|
| 2024 | ProjData.FF=reshape(FF,length(coord_y_proj),length(coord_x_proj)); |
---|
| 2025 | ProjData.ListVarName=[ProjData.ListVarName {'FF'}]; |
---|
| 2026 | ProjData.VarDimName=[ProjData.VarDimName {DimCell}]; |
---|
| 2027 | ProjData.VarAttribute{ivar_new+1+nbcoord}.Role='errorflag'; |
---|
| 2028 | end |
---|
| 2029 | end |
---|
| 2030 | |
---|
| 2031 | %% case of input fields defined on a structured grid |
---|
| 2032 | else |
---|
| 2033 | VarName=FieldData.ListVarName{VarIndex(1)};%get the first variable of the cell to get the input matrix dimensions |
---|
| 2034 | eval(['DimValue=size(FieldData.' VarName ');'])%input matrix dimensions |
---|
| 2035 | DimValue(DimValue==1)=[];%remove singleton dimensions |
---|
| 2036 | NbDim=numel(DimValue);%update number of space dimensions |
---|
| 2037 | nbcolor=1; %default number of 'color' components: third matrix index without corresponding coordinate |
---|
| 2038 | if NbDim>=3 |
---|
| 2039 | if NbDim>3 |
---|
| 2040 | errormsg='matrices with more than 3 dimensions not handled'; |
---|
| 2041 | return |
---|
| 2042 | else |
---|
| 2043 | if numel(find(VarType.coord))==2% the third matrix dimension does not correspond to a space coordinate |
---|
| 2044 | nbcolor=DimValue(3); |
---|
| 2045 | DimValue(3)=[]; %number of 'color' components updated |
---|
| 2046 | NbDim=2;% space dimension set to 2 |
---|
| 2047 | end |
---|
| 2048 | end |
---|
| 2049 | end |
---|
| 2050 | AYName=FieldData.ListVarName{VarType.coord(NbDim-1)};%name of input x coordinate (name preserved on projection) |
---|
| 2051 | AXName=FieldData.ListVarName{VarType.coord(NbDim)};%name of input y coordinate (name preserved on projection) |
---|
| 2052 | eval(['AX=FieldData.' AXName ';']) |
---|
| 2053 | eval(['AY=FieldData.' AYName ';']) |
---|
| 2054 | ListDimName=FieldData.VarDimName{VarIndex(1)}; |
---|
| 2055 | ProjData.ListVarName=[ProjData.ListVarName {AYName} {AXName}]; %TODO: check if it already exists in Projdata (several cells) |
---|
| 2056 | ProjData.VarDimName=[ProjData.VarDimName {AYName} {AXName}]; |
---|
| 2057 | |
---|
| 2058 | % for idim=1:length(ListDimName) |
---|
| 2059 | % DimName=ListDimName{idim}; |
---|
| 2060 | % if strcmp(DimName,'rgb')||strcmp(DimName,'nb_coord')||strcmp(DimName,'nb_coord_i') |
---|
| 2061 | % nbcolor=DimValue(idim); |
---|
| 2062 | % DimValue(idim)=[]; |
---|
| 2063 | % end |
---|
| 2064 | % if isequal(DimName,'nb_coord_j')% NOTE: CASE OF TENSOR NOT TREATED |
---|
| 2065 | % DimValue(idim)=[]; |
---|
| 2066 | % end |
---|
| 2067 | % end |
---|
| 2068 | Coord_z=[]; |
---|
| 2069 | Coord_y=[]; |
---|
| 2070 | Coord_x=[]; |
---|
| 2071 | |
---|
| 2072 | for idim=1:NbDim %loop on space dimensions |
---|
| 2073 | test_interp(idim)=0;%test for coordiate interpolation (non regular grid), =0 by default |
---|
| 2074 | ivar=VarType.coord(idim);% index of the variable corresponding to the current dimension |
---|
| 2075 | if ~isequal(ivar,0)% a variable corresponds to the dimension #idim |
---|
| 2076 | eval(['Coord{idim}=FieldData.' FieldData.ListVarName{ivar} ';']) ;% coord values for the input field |
---|
| 2077 | if numel(Coord{idim})==2 %input array defined on a regular grid |
---|
| 2078 | DCoord_min(idim)=(Coord{idim}(2)-Coord{idim}(1))/DimValue(idim); |
---|
| 2079 | else |
---|
| 2080 | DCoord=diff(Coord{idim});%array of coordinate derivatives for the input field |
---|
| 2081 | DCoord_min(idim)=min(DCoord); |
---|
| 2082 | DCoord_max=max(DCoord); |
---|
| 2083 | % test_direct(idim)=DCoord_max>0;% =1 for increasing values, 0 otherwise |
---|
| 2084 | if abs(DCoord_max-DCoord_min(idim))>abs(DCoord_max/1000) |
---|
| 2085 | msgbox_uvmat('ERROR',['non monotonic dimension variable # ' num2str(idim) ' in proj_field.m']) |
---|
| 2086 | return |
---|
| 2087 | end |
---|
| 2088 | test_interp(idim)=(DCoord_max-DCoord_min(idim))> 0.0001*abs(DCoord_max);% test grid regularity |
---|
| 2089 | end |
---|
| 2090 | test_direct(idim)=(DCoord_min(idim)>0); |
---|
| 2091 | else % no variable associated with the dimension #idim, the coordinate value is set equal to the matrix index by default |
---|
| 2092 | Coord_i_str=['Coord_' num2str(idim)]; |
---|
| 2093 | DCoord_min(idim)=1;%default |
---|
| 2094 | Coord{idim}=[0.5 DimValue(idim)-0.5]; |
---|
| 2095 | test_direct(idim)=1; |
---|
| 2096 | end |
---|
| 2097 | end |
---|
| 2098 | if DY==0 |
---|
| 2099 | DY=abs(DCoord_min(NbDim-1)); |
---|
| 2100 | end |
---|
| 2101 | npY=1+round(abs(Coord{NbDim-1}(end)-Coord{NbDim-1}(1))/DY);%nbre of points after interpol |
---|
| 2102 | if DX==0 |
---|
| 2103 | DX=abs(DCoord_min(NbDim)); |
---|
| 2104 | end |
---|
| 2105 | npX=1+round(abs(Coord{NbDim}(end)-Coord{NbDim}(1))/DX);%nbre of points after interpol |
---|
| 2106 | for idim=1:NbDim |
---|
| 2107 | if test_interp(idim) |
---|
| 2108 | DimValue(idim)=1+round(abs(Coord{idim}(end)-Coord{idim}(1))/abs(DCoord_min(idim)));%nbre of points after possible interpolation on a regular gri |
---|
| 2109 | end |
---|
| 2110 | end |
---|
| 2111 | Coord_y=linspace(Coord{NbDim-1}(1),Coord{NbDim-1}(end),npY); |
---|
| 2112 | test_direct_y=test_direct(NbDim-1); |
---|
| 2113 | Coord_x=linspace(Coord{NbDim}(1),Coord{NbDim}(end),npX); |
---|
| 2114 | test_direct_x=test_direct(NbDim); |
---|
| 2115 | DAX=DCoord_min(NbDim); |
---|
| 2116 | DAY=DCoord_min(NbDim-1); |
---|
| 2117 | minAX=min(Coord_x); |
---|
| 2118 | maxAX=max(Coord_x); |
---|
| 2119 | minAY=min(Coord_y); |
---|
| 2120 | maxAY=max(Coord_y); |
---|
| 2121 | xcorner=[minAX maxAX minAX maxAX]-ObjectData.Coord(1,1); |
---|
| 2122 | ycorner=[maxAY maxAY minAY minAY]-ObjectData.Coord(1,2); |
---|
| 2123 | xcor_new=xcorner*cos(Phi)+ycorner*sin(Phi);%coord new frame |
---|
| 2124 | ycor_new=-xcorner*sin(Phi)+ycorner*cos(Phi); |
---|
| 2125 | if ~testXMax |
---|
| 2126 | XMax=max(xcor_new); |
---|
| 2127 | end |
---|
| 2128 | if ~testXMin |
---|
| 2129 | XMin=min(xcor_new); |
---|
| 2130 | end |
---|
| 2131 | if ~testYMax |
---|
| 2132 | YMax=max(ycor_new); |
---|
| 2133 | end |
---|
| 2134 | if ~testYMin |
---|
| 2135 | YMin=min(ycor_new); |
---|
| 2136 | end |
---|
| 2137 | DXinit=(maxAX-minAX)/(DimValue(NbDim)-1); |
---|
| 2138 | DYinit=(maxAY-minAY)/(DimValue(NbDim-1)-1); |
---|
| 2139 | if DX==0 |
---|
| 2140 | DX=DXinit; |
---|
| 2141 | end |
---|
| 2142 | if DY==0 |
---|
| 2143 | DY=DYinit; |
---|
| 2144 | end |
---|
| 2145 | if NbDim==3 |
---|
| 2146 | DZ=(Coord{1}(end)-Coord{1}(1))/(DimValue(1)-1); |
---|
| 2147 | if ~test_direct(1) |
---|
| 2148 | DZ=-DZ; |
---|
| 2149 | end |
---|
| 2150 | Coord_z=linspace(Coord{1}(1),Coord{1}(end),DimValue(1)); |
---|
| 2151 | test_direct_z=test_direct(1); |
---|
| 2152 | end |
---|
| 2153 | npX=floor((XMax-XMin)/DX+1); |
---|
| 2154 | npY=floor((YMax-YMin)/DY+1); |
---|
| 2155 | if test_direct_y |
---|
| 2156 | coord_y_proj=linspace(YMin,YMax,npY);%abscissa of the new pixels along the line |
---|
| 2157 | else |
---|
| 2158 | coord_y_proj=linspace(YMax,YMin,npY);%abscissa of the new pixels along the line |
---|
| 2159 | end |
---|
| 2160 | if test_direct_x |
---|
| 2161 | coord_x_proj=linspace(XMin,XMax,npX);%abscissa of the new pixels along the line |
---|
| 2162 | else |
---|
| 2163 | coord_x_proj=linspace(XMax,XMin,npX);%abscissa of the new pixels along the line |
---|
| 2164 | end |
---|
| 2165 | |
---|
| 2166 | % case with no rotation and interpolation |
---|
| 2167 | if isequal(ProjMode,'projection') && isequal(Phi,0) && isequal(Theta,0) && isequal(Psi,0) |
---|
| 2168 | if ~testXMin && ~testXMax && ~testYMin && ~testYMax && NbDim==2 |
---|
| 2169 | ProjData=FieldData; |
---|
| 2170 | else |
---|
| 2171 | indY=NbDim-1; |
---|
| 2172 | if test_direct(indY) |
---|
| 2173 | min_indy=ceil((YMin-Coord{indY}(1))/DYinit)+1; |
---|
| 2174 | YIndexFirst=floor((YMax-Coord{indY}(1))/DYinit)+1; |
---|
| 2175 | Ybound(1)=Coord{indY}(1)+DYinit*(min_indy-1); |
---|
| 2176 | Ybound(2)=Coord{indY}(1)+DYinit*(YIndexFirst-1); |
---|
| 2177 | else |
---|
| 2178 | min_indy=ceil((Coord{indY}(1)-YMax)/DYinit)+1; |
---|
| 2179 | max_indy=floor((Coord{indY}(1)-YMin)/DYinit)+1; |
---|
| 2180 | Ybound(2)=Coord{indY}(1)-DYinit*(max_indy-1); |
---|
| 2181 | Ybound(1)=Coord{indY}(1)-DYinit*(min_indy-1); |
---|
| 2182 | end |
---|
| 2183 | if test_direct(NbDim)==1 |
---|
| 2184 | min_indx=ceil((XMin-Coord{NbDim}(1))/DXinit)+1; |
---|
| 2185 | max_indx=floor((XMax-Coord{NbDim}(1))/DXinit)+1; |
---|
| 2186 | Xbound(1)=Coord{NbDim}(1)+DXinit*(min_indx-1); |
---|
| 2187 | Xbound(2)=Coord{NbDim}(1)+DXinit*(max_indx-1); |
---|
| 2188 | else |
---|
| 2189 | min_indx=ceil((Coord{NbDim}(1)-XMax)/DXinit)+1; |
---|
| 2190 | max_indx=floor((Coord{NbDim}(1)-XMin)/DXinit)+1; |
---|
| 2191 | Xbound(2)=Coord{NbDim}(1)+DXinit*(max_indx-1); |
---|
| 2192 | Xbound(1)=Coord{NbDim}(1)+DXinit*(min_indx-1); |
---|
| 2193 | end |
---|
| 2194 | if NbDim==3 |
---|
| 2195 | DimCell(1)=[]; %suppress z variable |
---|
| 2196 | DimValue(1)=[]; |
---|
| 2197 | %structured coordinates |
---|
| 2198 | if test_direct(1) |
---|
| 2199 | iz=ceil((ObjectData.Coord(1,3)-Coord{1}(1))/DZ)+1; |
---|
| 2200 | else |
---|
| 2201 | iz=ceil((Coord{1}(1)-ObjectData.Coord(1,3))/DZ)+1; |
---|
| 2202 | end |
---|
| 2203 | end |
---|
| 2204 | min_indy=max(min_indy,1);% deals with margin (bound lower than the first index) |
---|
| 2205 | min_indx=max(min_indx,1); |
---|
| 2206 | max_indy=min(max_indy,DimValue(1)); |
---|
| 2207 | max_indx=min(max_indx,DimValue(2)); |
---|
| 2208 | for ivar=VarIndex% loop on non coordinate variables |
---|
| 2209 | VarName=FieldData.ListVarName{ivar}; |
---|
| 2210 | ProjData.ListVarName=[ProjData.ListVarName VarName]; |
---|
| 2211 | ProjData.VarDimName=[ProjData.VarDimName {DimCell}]; |
---|
| 2212 | if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute)>=ivar |
---|
| 2213 | ProjData.VarAttribute{length(ProjData.ListVarName)}=FieldData.VarAttribute{ivar}; |
---|
| 2214 | end |
---|
| 2215 | if NbDim==3 |
---|
| 2216 | eval(['ProjData.' VarName '=squeeze(FieldData.' VarName '(iz,min_indy:max_indy,min_indx:max_indx));']); |
---|
| 2217 | else |
---|
| 2218 | eval(['ProjData.' VarName '=FieldData.' VarName '(min_indy:max_indy,min_indx:max_indx,:);']); |
---|
| 2219 | end |
---|
| 2220 | end |
---|
| 2221 | eval(['ProjData.' AYName '=[Ybound(1) Ybound(2)];']) %record the new (projected ) y coordinates |
---|
| 2222 | eval(['ProjData.' AXName '=[Xbound(1) Xbound(2)];']) %record the new (projected ) x coordinates |
---|
| 2223 | end |
---|
| 2224 | elseif isfield(FieldData,'A') %TO GENERALISE % case with rotation and/or interpolation |
---|
| 2225 | if NbDim==2 %2D case |
---|
| 2226 | [X,Y]=meshgrid(coord_x_proj,coord_y_proj);%grid in the new coordinates |
---|
| 2227 | XIMA=ObjectData.Coord(1,1)+(X)*cos(Phi)-Y*sin(Phi);%corresponding coordinates in the original image |
---|
| 2228 | YIMA=ObjectData.Coord(1,2)+(X)*sin(Phi)+Y*cos(Phi); |
---|
| 2229 | XIMA=(XIMA-minAX)/DXinit+1;% image index along x |
---|
| 2230 | YIMA=(-YIMA+maxAY)/DYinit+1;% image index along y |
---|
| 2231 | XIMA=reshape(round(XIMA),1,npX*npY);%indices reorganized in 'line' |
---|
| 2232 | YIMA=reshape(round(YIMA),1,npX*npY); |
---|
| 2233 | flagin=XIMA>=1 & XIMA<=DimValue(2) & YIMA >=1 & YIMA<=DimValue(1);%flagin=1 inside the original image |
---|
| 2234 | if isequal(ObjectData.ProjMode,'interp_tps') |
---|
| 2235 | npx_interp_tps=ceil(abs(DX/DAX)); |
---|
| 2236 | npy_interp_tps=ceil(abs(DY/DAY)); |
---|
| 2237 | Minterp_tps=ones(npy_interp_tps,npx_interp_tps)/(npx_interp_tps*npy_interp_tps); |
---|
| 2238 | test_interp_tps=1; |
---|
| 2239 | else |
---|
| 2240 | test_interp_tps=0; |
---|
| 2241 | end |
---|
| 2242 | eval(['ProjData.' AYName '=[coord_y_proj(1) coord_y_proj(end)];']) %record the new (projected ) y coordinates |
---|
| 2243 | eval(['ProjData.' AXName '=[coord_x_proj(1) coord_x_proj(end)];']) %record the new (projected ) x coordinates |
---|
| 2244 | for ivar=VarIndex |
---|
| 2245 | VarName=FieldData.ListVarName{ivar}; |
---|
| 2246 | if test_interp(1) || test_interp(2)%interpolate on a regular grid |
---|
| 2247 | eval(['ProjData.' VarName '=interp2(Coord{2},Coord{1},FieldData.' VarName ',Coord_x,Coord_y'');']) %TO TEST |
---|
| 2248 | end |
---|
| 2249 | %filter the field (image) if option 'interp_tps' is used |
---|
| 2250 | if test_interp_tps |
---|
| 2251 | Aclass=class(FieldData.A); |
---|
| 2252 | ProjData.(VarName)=interp_tps2(Minterp_tps,FieldData.(VarName),'valid'); |
---|
| 2253 | if ~isequal(Aclass,'double') |
---|
| 2254 | ProjData.(VarName)=Aclass(FieldData.(VarName));%revert to integer values |
---|
| 2255 | end |
---|
| 2256 | end |
---|
| 2257 | eval(['vec_A=reshape(FieldData.' VarName ',[],nbcolor);'])%put the original image in line |
---|
| 2258 | %ind_in=find(flagin); |
---|
| 2259 | ind_out=find(~flagin); |
---|
| 2260 | ICOMB=(XIMA-1)*DimValue(1)+YIMA; |
---|
| 2261 | ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A |
---|
| 2262 | vec_B(flagin,1:nbcolor)=vec_A(ICOMB,:); |
---|
| 2263 | for icolor=1:nbcolor |
---|
| 2264 | vec_B(ind_out,icolor)=zeros(size(ind_out)); |
---|
| 2265 | end |
---|
| 2266 | ProjData.ListVarName=[ProjData.ListVarName VarName]; |
---|
| 2267 | ProjData.VarDimName=[ProjData.VarDimName {DimCell}]; |
---|
| 2268 | if isfield(FieldData,'VarAttribute')&&length(FieldData.VarAttribute)>=ivar |
---|
| 2269 | ProjData.VarAttribute{length(ProjData.ListVarName)+nbcoord}=FieldData.VarAttribute{ivar}; |
---|
| 2270 | end |
---|
| 2271 | eval(['ProjData.' VarName '=reshape(vec_B,npY,npX,nbcolor);']); |
---|
| 2272 | end |
---|
| 2273 | ProjData.FF=reshape(~flagin,npY,npX);%false flag A FAIRE: tenir compte d'un flga antï¿œrieur |
---|
| 2274 | ProjData.ListVarName=[ProjData.ListVarName 'FF']; |
---|
| 2275 | ProjData.VarDimName=[ProjData.VarDimName {DimCell}]; |
---|
| 2276 | ProjData.VarAttribute{length(ProjData.ListVarName)}.Role='errorflag'; |
---|
| 2277 | else %3D case |
---|
| 2278 | if ~testangle |
---|
| 2279 | % unstructured z coordinate |
---|
| 2280 | test_sup=(Coord{1}>=ObjectData.Coord(1,3)); |
---|
| 2281 | iz_sup=find(test_sup); |
---|
| 2282 | iz=iz_sup(1); |
---|
| 2283 | if iz>=1 & iz<=npz |
---|
| 2284 | %ProjData.ListDimName=[ProjData.ListDimName ListDimName(2:end)]; |
---|
| 2285 | %ProjData.DimValue=[ProjData.DimValue npY npX]; |
---|
| 2286 | for ivar=VarIndex |
---|
| 2287 | VarName=FieldData.ListVarName{ivar}; |
---|
| 2288 | ProjData.ListVarName=[ProjData.ListVarName VarName]; |
---|
| 2289 | ProjData.VarAttribute{length(ProjData.ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes |
---|
| 2290 | eval(['ProjData.' VarName '=squeeze(FieldData.' VarName '(iz,:,:));'])% select the z index iz |
---|
| 2291 | %TODO : do a vertical average for a thick plane |
---|
| 2292 | if test_interp(2) || test_interp(3) |
---|
| 2293 | eval(['ProjData.' VarName '=interp2(Coord{3},Coord{2},ProjData.' VarName ',Coord_x,Coord_y'');']) |
---|
| 2294 | end |
---|
| 2295 | end |
---|
| 2296 | end |
---|
| 2297 | else |
---|
| 2298 | errormsg='projection of structured coordinates on oblique plane not yet implemented'; |
---|
| 2299 | %TODO: use interp3 |
---|
| 2300 | return |
---|
| 2301 | end |
---|
| 2302 | end |
---|
| 2303 | end |
---|
| 2304 | end |
---|
| 2305 | |
---|
| 2306 | %% projection of velocity components in the rotated coordinates |
---|
| 2307 | if testangle |
---|
| 2308 | if isempty(ivar_V) |
---|
| 2309 | msgbox_uvmat('ERROR','v velocity component missing in proj_field.m') |
---|
| 2310 | return |
---|
| 2311 | end |
---|
| 2312 | UName=FieldData.ListVarName{ivar_U}; |
---|
| 2313 | VName=FieldData.ListVarName{ivar_V}; |
---|
| 2314 | eval(['ProjData.' UName '=cos(Phi)*ProjData.' UName '+ sin(Phi)*ProjData.' VName ';']) |
---|
| 2315 | eval(['ProjData.' VName '=cos(Theta)*(-sin(Phi)*ProjData.' UName '+ cos(Phi)*ProjData.' VName ');']) |
---|
| 2316 | if ~isempty(ivar_W) |
---|
| 2317 | WName=FieldData.ListVarName{ivar_W}; |
---|
| 2318 | eval(['ProjData.' VName '=ProjData.' VName '+ ProjData.' WName '*sin(Theta);'])% |
---|
| 2319 | eval(['ProjData.' WName '=NormVec_X*ProjData.' UName '+ NormVec_Y*ProjData.' VName '+ NormVec_Z* ProjData.' WName ';']); |
---|
| 2320 | end |
---|
| 2321 | if ~isequal(Psi,0) |
---|
| 2322 | eval(['ProjData.' UName '=cos(Psi)* ProjData.' UName '- sin(Psi)*ProjData.' VName ';']); |
---|
| 2323 | eval(['ProjData.' VName '=sin(Psi)* ProjData.' UName '+ cos(Psi)*ProjData.' VName ';']); |
---|
| 2324 | end |
---|
| 2325 | end |
---|
| 2326 | end |
---|
| 2327 | |
---|
| 2328 | %------------------------------------------------------------------------ |
---|
| 2329 | %--- transfer the global attributes |
---|
| 2330 | function [ProjData,errormsg]=proj_heading(FieldData,ObjectData) |
---|
| 2331 | %------------------------------------------------------------------------ |
---|
| 2332 | ProjData=[];%default |
---|
| 2333 | errormsg='';%default |
---|
| 2334 | |
---|
| 2335 | %% transfer error |
---|
| 2336 | if isfield(FieldData,'Txt') |
---|
| 2337 | errormsg=FieldData.Txt; %transmit erreur message |
---|
| 2338 | return; |
---|
| 2339 | end |
---|
| 2340 | |
---|
| 2341 | %% transfer global attributes |
---|
| 2342 | if ~isfield(FieldData,'ListGlobalAttribute') |
---|
| 2343 | ProjData.ListGlobalAttribute={}; |
---|
| 2344 | else |
---|
| 2345 | ProjData.ListGlobalAttribute=FieldData.ListGlobalAttribute; |
---|
| 2346 | end |
---|
| 2347 | for iattr=1:length(ProjData.ListGlobalAttribute) |
---|
| 2348 | AttrName=ProjData.ListGlobalAttribute{iattr}; |
---|
| 2349 | if isfield(FieldData,AttrName) |
---|
| 2350 | ProjData.(AttrName)=FieldData.(AttrName); |
---|
| 2351 | end |
---|
| 2352 | end |
---|
| 2353 | |
---|
| 2354 | %% transfer coordinate unit |
---|
| 2355 | if isfield(ProjData,'CoordUnit') |
---|
| 2356 | ProjData=rmfield(ProjData,'CoordUnit');% do not transfer by default (to avoid x/y=1 for profiles) |
---|
| 2357 | end |
---|
| 2358 | if isfield(FieldData,'CoordUnit') |
---|
| 2359 | if isfield(ObjectData,'CoordUnit') && ~strcmp(FieldData.CoordUnit,ObjectData.CoordUnit) |
---|
| 2360 | errormsg=[ObjectData.Type ' in ' ObjectData.CoordUnit ' coordinates, while field in ' FieldData.CoordUnit ]; |
---|
| 2361 | return |
---|
| 2362 | elseif strcmp(ObjectData.Type,'plane')|| strcmp(ObjectData.Type,'volume') |
---|
| 2363 | ProjData.CoordUnit=FieldData.CoordUnit; |
---|
| 2364 | end |
---|
| 2365 | end |
---|
| 2366 | |
---|
| 2367 | %% store the properties of the projection object |
---|
| 2368 | ListObject={'Name','Type','ProjMode','angle','RangeX','RangeY','RangeZ','DX','DY','DZ','Coord'}; |
---|
| 2369 | for ilist=1:length(ListObject) |
---|
| 2370 | if isfield(ObjectData,ListObject{ilist}) |
---|
| 2371 | val=ObjectData.(ListObject{ilist}); |
---|
| 2372 | if ~isempty(val) |
---|
| 2373 | ProjData.(['ProjObject' ListObject{ilist}])=val; |
---|
| 2374 | ProjData.ListGlobalAttribute=[ProjData.ListGlobalAttribute {['ProjObject' ListObject{ilist}]}]; |
---|
| 2375 | end |
---|
| 2376 | end |
---|
| 2377 | end |
---|
| 2378 | |
---|