source: trunk/src/proj_field.m @ 748

Last change on this file since 748 was 748, checked in by sommeria, 10 years ago

update for 3D plots, panel Coordiantes introduces, while coordiantes now called Axes

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