source: trunk/src/proj_field.m @ 533

Last change on this file since 533 was 533, checked in by sommeria, 9 years ago

bugs repaired in proj_field

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