source: trunk/src/proj_field.m @ 507

Last change on this file since 507 was 494, checked in by sommeria, 12 years ago

various bugs corrected after testing in Windows OS. Introduction
of filter tps

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