source: trunk/src/proj_field.m @ 399

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

implementation of thin plate interpolation (proj on planes with mode 'filter'), rationalisation of variable formats in civ_matlab

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