[646]  1  %'mask_proj': restrict input fields to a mask region, set to 0 outside


[639]  2  %


[646]  3  % function [ProjData,errormsg]=mask_proj(FieldData,MaskData)


[809]  4 


 5  %=======================================================================


[908]  6  % Copyright 20082015, LEGI UMR 5519 / CNRS UJF GINP, Grenoble, France


[809]  7  % http://www.legi.grenobleinp.fr


 8  % Joel.Sommeria  Joel.Sommeria (A) legi.cnrs.fr


[639]  9  %


 10  % This file is part of the toolbox UVMAT.


[809]  11  %


[639]  12  % UVMAT is free software; you can redistribute it and/or modify


[809]  13  % it under the terms of the GNU General Public License as published


 14  % by the Free Software Foundation; either version 2 of the license,


 15  % or (at your option) any later version.


 16  %


[639]  17  % UVMAT is distributed in the hope that it will be useful,


 18  % but WITHOUT ANY WARRANTY; without even the implied warranty of


 19  % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the


[809]  20  % GNU General Public License (see LICENSE.txt) for more details.


 21  %=======================================================================


[639]  22 


 23  function [ProjData,errormsg]=mask_proj(FieldData,MaskData)


 24  errormsg='';%default


 25  ProjData=FieldData;


 26 


 27 


 28  %% group the variables (fields of 'FieldData') in cells of variables with the same dimensions


 29  [CellInfo,NbDimArray,errormsg]=find_field_cells(FieldData);


 30  if ~isempty(errormsg)


 31  errormsg=['error in proj_field/proj_plane:' errormsg];


 32  return


 33  end


 34  [Npy,Npx]=size(MaskData.A);


[675]  35 


[639]  36  for icell=1:numel(CellInfo)


 37  if NbDimArray(icell)==2


[672]  38  if isfield(CellInfo{icell},'VarIndex_errorflag')


 39  FFName=FieldData.ListVarName{CellInfo{icell}.VarIndex_errorflag};


[675]  40  check_new_error_flag=0;


[672]  41  else


 42  FFName='FF';%default error (mask) flag name (if not already used)


 43  if isfield(FieldData,'FF')


 44  ind=1;


 45  while isfield(FieldData,['FF_' num2str(ind)])


 46  ind=ind+1;


 47  end


 48  FFName=['FF_' num2str(ind)];% append an index to the name of error flag, FF_1,FF_2...


 49  end


 50  ProjData.ListVarName=[FieldData.ListVarName {FFName}];


[675]  51  ProjData.VarAttribute{numel(ProjData.ListVarName)}.Role='errorflag';


 52  check_new_error_flag=1;


[672]  53  end


[639]  54  switch CellInfo{icell}.CoordType;


[675]  55  case 'scattered'


 56  XName=FieldData.ListVarName{CellInfo{icell}.Coord_x};


 57  YName=FieldData.ListVarName{CellInfo{icell}.Coord_y};


[782]  58  DX=(MaskData.Coord_x(2)MaskData.Coord_x(1))/(Npx1);


 59  DY=(MaskData.Coord_y(2)MaskData.Coord_y(1))/(Npy1);


 60  mask_ind_i=round(0.5+(FieldData.(XName)MaskData.Coord_x(1))/DX);%nbpoint terms


 61  mask_ind_j=round(0.5+(FieldData.(YName)MaskData.Coord_y(1))/DY);%nbpoint terms


[639]  62  checkin=mask_ind_j+Npy*(mask_ind_i1);%array of mask indices for the nbpoints


 63  checkin=checkin(mask_ind_i>=1 & mask_ind_i<=Npx & mask_ind_j>=1 & mask_ind_j<=Npy);%reduced array of mask indices (inside the image)


 64  checkfalse=true(size(FieldData.(XName)));


 65  MaskData.A=reshape(MaskData.A,1,[]);


 66  checkfalse(MaskData.A(checkin)>200)=0;


 67  for ivar=1:numel(CellInfo{icell}.VarIndex)


 68  VarName=FieldData.ListVarName{CellInfo{icell}.VarIndex(ivar)};


[672]  69  ProjData.(VarName)(checkfalse)=0;


[639]  70  end


[675]  71  if check_new_error_flag% an error flag needs to be created in the current cell


[672]  72  ProjData.(FFName)=zeros(size(ProjData.(VarName)));


[675]  73  ProjData.VarDimName=[FieldData.VarDimName FieldData.VarDimName(CellInfo{icell}.VarIndex(1))];


[672]  74  end


[675]  75  ProjData.(FFName)(checkfalse)=1;% update the existing error flag


[639]  76  case 'grid'


[675]  77  XName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(2)};


 78  YName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(1)};


[639]  79  Var1Name=FieldData.ListVarName{CellInfo{icell}.VarIndex(1)};


 80  [Npy_field,Npx_field]=size(FieldData.(Var1Name));


[675]  81  XArray=linspace(FieldData.(XName)(1),FieldData.(XName)(end),Npx_field);


 82  YArray=linspace(FieldData.(YName)(1),FieldData.(YName)(end),Npy_field);


[782]  83  XMask=linspace(MaskData.Coord_x(1),MaskData.Coord_x(end),Npx);


 84  YMask=linspace(MaskData.Coord_y(1),MaskData.Coord_y(end),Npy);


[639]  85  [XMask,YMask]=meshgrid(XMask,YMask);


 86  Mask = interp2(XMask,YMask,MaskData.A,XArray,YArray','nearest');


[672]  87  Mask=Mask>200;


[639]  88  for ivar=1:numel(CellInfo{icell}.VarIndex)


 89  VarName=FieldData.ListVarName{CellInfo{icell}.VarIndex(ivar)};


[672]  90  if ~strcmp(VarName,FFName)


 91  ProjData.(VarName)=FieldData.(VarName).*Mask;


 92  end


[639]  93  end


[675]  94  if check_new_error_flag% an error flag needs to be created in the current cell


 95  ProjData.(FFName)= ~Mask;


 96  ProjData.VarDimName=[FieldData.VarDimName FieldData.VarDimName(CellInfo{icell}.VarIndex(1))];


 97  else


[672]  98  ProjData.(FFName)=FieldData.(FFName)  ~Mask;


 99  end


[639]  100  end


 101  end


 102  end


 103 


[672]  104 

