source: trunk/src/proj_field.m @ 698

Last change on this file since 698 was 690, checked in by sommeria, 11 years ago

representation of ancillary data as table improved.

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