Home > FA-MVEMD > 3D > MinimaMaxima3D.m

# MinimaMaxima3D

## PURPOSE

V 1.0 Dec 13, 07

## SYNOPSIS

function [Maxima,MaxPos,Minima,MinPos]=MinimaMaxima3D(Input,Robust,LookInBoundaries,numbermax,numbermin)

## DESCRIPTION

     V 1.0 Dec 13, 07
     Author Sam Pichardo.
     This  function finds the local minima and maxima in a 3D Cartesian data.
     It's assumed that the data is uniformly distributed.
     The minima and maxima are calculated using a multi-directional derivation.

     Use:

      [Maxima,MaxPos,Minima,MinPos]=MinimaMaxima3D(Input,[Robust],[LookInBoundaries],[numbermax],[numbermin])

     where Input is the 3D data and Robust (optional and with a default value
     of 1) indicates if the multi-directional derivation should include the
     diagonal derivations.

     Input has to have a size larger or equal than [3 x 3 x 3]

     If Robust=1, the total number of derivations taken into account are 26: 6
     for all surrounding elements colliding each of the faces of the unit cube;
     10 for all the surrounding elements in diagonal.

     If Robust =0, then only the 6 elements of the colliding faces are considered

     The function returns in Maxima and MaxPos, respectively,
     the values (numbermax) and subindexes (numbermax x 3) of local maxima
     and position in Input. Maxima (and the subindexes) are sorted in
     descending order.
     Similar situation for Minima and MinimaPos witn a numbermin elements but
     with the execption of being sorted in ascending order.

     IMPORTANT: if numbermin or numbermax are not specified, ALL the minima
     or maxima will be returned. This can be a useless for highly
     oscillating data

     LookInBoundaries (default value of 0) specifies if a search of the minima/maxima should be
     done in the boundaries of the matrix. This situation depends on the
     the desire application. When it is not activated, the algorithm WILL NOT
     FIND ANY MINIMA/MAXIMA on the 6 layers of the boundaries.
     When it is activated, the finding minima and maxima on the boundaries is done by
     replicating the extra layer as the layer 2 (or layer N-1, depending of the boundary)
     By example (and using a 2D matrix for simplicity reasons):
     For the matrix
     [ 4 1 3 7
       5 7 8 8
       9 9 9 9
       5 6 7 9]

     the calculation of the partial derivate following the -x direction will be done by substrascting
     [ 5 7 8 8
       4 1 3 7
       5 7 8 8
       9 9 9 9]
     to the input. And so on for the other dimensions.
     Like this, the value "1" at the coordinate (1,2) will be detected as a
     minima. Same situation for the value "5" at the coordinate (4,1)

## CROSS-REFERENCE INFORMATION

This function calls:

This function is called by:

- EMD3D2V_parallel Purpose:
- EMD3D3V Purpose:
- EMD3D3V_parallel Purpose:
- EMD3D3V_parallel_var Purpose:

## SOURCE CODE

    0001 function [Maxima,MaxPos,Minima,MinPos]=MinimaMaxima3D(Input,Robust,LookInBoundaries,numbermax,numbermin)
    0002 % V 1.0 Dec 13, 07
    0003 % Author Sam Pichardo.
    0004 % This  function finds the local minima and maxima in a 3D Cartesian data.
    0005 % It's assumed that the data is uniformly distributed.
    0006 % The minima and maxima are calculated using a multi-directional derivation.
    0007 %
    0008 % Use:
    0009 %
    0010 %  [Maxima,MaxPos,Minima,MinPos]=MinimaMaxima3D(Input,[Robust],[LookInBoundaries],[numbermax],[numbermin])
    0011 %
    0012 % where Input is the 3D data and Robust (optional and with a default value
    0013 % of 1) indicates if the multi-directional derivation should include the
    0014 % diagonal derivations.
    0015 %
    0016 % Input has to have a size larger or equal than [3 x 3 x 3]
    0017 %
    0018 % If Robust=1, the total number of derivations taken into account are 26: 6
    0019 % for all surrounding elements colliding each of the faces of the unit cube;
    0020 % 10 for all the surrounding elements in diagonal.
    0021 %
    0022 % If Robust =0, then only the 6 elements of the colliding faces are considered
    0023 %
    0024 % The function returns in Maxima and MaxPos, respectively,
    0025 % the values (numbermax) and subindexes (numbermax x 3) of local maxima
    0026 % and position in Input. Maxima (and the subindexes) are sorted in
    0027 % descending order.
    0028 % Similar situation for Minima and MinimaPos witn a numbermin elements but
    0029 % with the execption of being sorted in ascending order.
    0030 %
    0031 % IMPORTANT: if numbermin or numbermax are not specified, ALL the minima
    0032 % or maxima will be returned. This can be a useless for highly
    0033 % oscillating data
    0034 %
    0035 % LookInBoundaries (default value of 0) specifies if a search of the minima/maxima should be
    0036 % done in the boundaries of the matrix. This situation depends on the
    0037 % the desire application. When it is not activated, the algorithm WILL NOT
    0038 % FIND ANY MINIMA/MAXIMA on the 6 layers of the boundaries.
    0039 % When it is activated, the finding minima and maxima on the boundaries is done by
    0040 % replicating the extra layer as the layer 2 (or layer N-1, depending of the boundary)
    0041 % By example (and using a 2D matrix for simplicity reasons):
    0042 % For the matrix
    0043 % [ 4 1 3 7
    0044 %   5 7 8 8
    0045 %   9 9 9 9
    0046 %   5 6 7 9]
    0047 %
    0048 % the calculation of the partial derivate following the -x direction will be done by substrascting
    0049 % [ 5 7 8 8
    0050 %   4 1 3 7
    0051 %   5 7 8 8
    0052 %   9 9 9 9]
    0053 % to the input. And so on for the other dimensions.
    0054 % Like this, the value "1" at the coordinate (1,2) will be detected as a
    0055 % minima. Same situation for the value "5" at the coordinate (4,1)
    0056
    0057
    0058 if nargin <1
    0059     test=load('temp.mat');
    0060     pf=test.uresTot(test.EvalLims(2,1):test.EvalLims(2,2));
    0061     pf=reshape(pf,length(test.EvalCoord{2}.Ry),length(test.EvalCoord{2}.Rx),length(test.EvalCoord{2}.Rz));
    0062     Input = abs(pf)*1.5e6;
    0063     clear test;
    0064     clear pf;
    0065     Robust =1;
    0066 end
    0067
    0068 Asize=size(Input);
    0069
    0070 if length(Asize)<3
    0071     error('MinimaMaxima3D can only works with 3D matrices ');
    0072 end
    0073
    0074
    0075 if (Asize(1)<3 || Asize(2)<3 || Asize(3)<3)
    0076     error('MinimaMaxima3D can only works with matrices with dimensions equal or larger to [3x3x3]');
    0077 end
    0078
    0079 if ~isreal(Input)
    0080     warning('ATTENTION, complex values detected!!, using abs(Input)');
    0081     Input=abs(Input);
    0082 end
    0083
    0084 if ~exist('Robust','var')
    0085     Robust=1;
    0086 end
    0087
    0088 if ~exist('LookInBoundaries','var')
    0089     LookInBoundaries=0;
    0090 end
    0091
    0092 if ~exist('numbermax','var')
    0093     numbermax=0;
    0094 end
    0095
    0096 if ~exist('numbermin','var')
    0097     numbermin=0;
    0098 end
    0099
    0100 [xx_base,yy_base,zz_base]=ndgrid(1:Asize(1),1:Asize(2),1:Asize(3));
    0101
    0102
    0103 IndBase=sub2ind(Asize,xx_base(:),yy_base(:),zz_base(:));
    0104
    0105 if Robust ~= 0
    0106     Numbder_dd=26;
    0107 else
    0108     Numbder_dd=6;
    0109 end
    0110
    0111 if LookInBoundaries==0
    0112     lx=1:Asize(1);
    0113     lx_p1=[2:Asize(1),Asize(1)];
    0114     lx_m1=[1,1:Asize(1)-1];
    0115     ly=1:Asize(2);
    0116     ly_p1=[2:Asize(2),Asize(2)];
    0117     ly_m1=[1,1:Asize(2)-1];
    0118     lz=1:Asize(3);
    0119     lz_p1=[2:Asize(3),Asize(3)];
    0120     lz_m1=[1,1:Asize(3)-1];
    0121 else
    0122     lx=1:Asize(1);
    0123     lx_p1=[2:Asize(1),Asize(1)-1]; %We replicate the layer N-1 as the layer N+1
    0124     lx_m1=[2,1:Asize(1)-1]; %We replicate the layer 2 as the layer -1
    0125     ly=1:Asize(2);
    0126     ly_p1=[2:Asize(2),Asize(2)-1]; %We replicate the layer N-1 as the layer N+1
    0127     ly_m1=[2,1:Asize(2)-1]; %We replicate the layer 2 as the layer -1
    0128     lz=1:Asize(3);
    0129     lz_p1=[2:Asize(3),Asize(3)-1]; %We replicate the layer N-1 as the layer N+1
    0130     lz_m1=[2,1:Asize(3)-1];%We replicate the layer 2 as the layer -1
    0131 end
    0132
    0133 for n_dd=1:Numbder_dd
    0134     switch n_dd
    0135         case 1
    0136             %%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x+1)
    0137             [xx,yy,zz]=ndgrid(lx_p1,ly,lz);
    0138
    0139         case 2
    0140             %%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x-1)
    0141             [xx,yy,zz]=ndgrid(lx_m1,ly,lz);
    0142
    0143         case 3
    0144             %%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(y)-elem(y+1)
    0145             [xx,yy,zz]=ndgrid(lx,ly_p1,lz);
    0146
    0147         case 4
    0148             %%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(y)-elem(y-1)
    0149             [xx,yy,zz]=ndgrid(lx,ly_m1,lz);
    0150
    0151         case 5
    0152             %%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(z)-elem(z+1)
    0153             [xx,yy,zz]=ndgrid(lx,ly,lz_p1);
    0154
    0155          case 6
    0156             %%%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(z)-elem(z-1)
    0157             [xx,yy,zz]=ndgrid(lx,ly,lz_m1);
    0158         case 7
    0159            %%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x+1,y+1)
    0160             [xx,yy,zz]=ndgrid(lx_p1,ly_p1,lz);
    0161         case 8
    0162            %%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x+1,y-1)
    0163             [xx,yy,zz]=ndgrid(lx_p1,ly_m1,lz);
    0164         case 9
    0165            %%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x-1,y-1)
    0166             [xx,yy,zz]=ndgrid(lx_m1,ly_m1,lz);
    0167         case 10
    0168            %%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x-1,y+1)
    0169             [xx,yy,zz]=ndgrid(lx_m1,ly_p1,lz);
    0170         case 11
    0171            %%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x+1,z+1)
    0172             [xx,yy,zz]=ndgrid(lx_p1,ly,lz_p1);
    0173         case 12
    0174            %%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x+1,z-1)
    0175             [xx,yy,zz]=ndgrid(lx_p1,ly,lz_m1);
    0176         case 13
    0177            %%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x-1,z-1)
    0178             [xx,yy,zz]=ndgrid(lx_m1,ly,lz_m1);
    0179         case 14
    0180            %%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x-1,z+1)
    0181             [xx,yy,zz]=ndgrid(lx_m1,ly,lz_p1);
    0182         case 15
    0183            %%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(y+1,z+1)
    0184             [xx,yy,zz]=ndgrid(lx,ly_p1,lz_p1);
    0185         case 16
    0186            %%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(y+1,z-1)
    0187             [xx,yy,zz]=ndgrid(lx,ly_p1,lz_m1);
    0188         case 17
    0189            %%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(y-1,z-1)
    0190             [xx,yy,zz]=ndgrid(lx,ly_m1,lz_m1);
    0191         case 18
    0192            %%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(y-1,z+1)
    0193             [xx,yy,zz]=ndgrid(lx,ly_m1,lz_p1);
    0194          case 19
    0195            %%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x+1,y+1,z+1)
    0196             [xx,yy,zz]=ndgrid(lx_p1,ly_p1,lz_p1);
    0197          case 20
    0198            %%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x+1,y+1,z-1)
    0199             [xx,yy,zz]=ndgrid(lx_p1,ly_p1,lz_m1);
    0200          case 21
    0201            %%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x+1,y-1,z+1)
    0202             [xx,yy,zz]=ndgrid(lx_p1,ly_m1,lz_p1);
    0203          case 22
    0204            %%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x+1,y-1,z-1)
    0205             [xx,yy,zz]=ndgrid(lx_p1,ly_m1,lz_m1);
    0206          case 23
    0207            %%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x-1,y+1,z+1)
    0208             [xx,yy,zz]=ndgrid(lx_m1,ly_p1,lz_p1);
    0209          case 24
    0210            %%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x-1,y+1,z-1)
    0211             [xx,yy,zz]=ndgrid(lx_m1,ly_p1,lz_m1);
    0212          case 25
    0213            %%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x-1,y-1,z+1)
    0214             [xx,yy,zz]=ndgrid(lx_m1,ly_m1,lz_p1);
    0215          case 26
    0216            %%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x-1,y-1,z-1)
    0217             [xx,yy,zz]=ndgrid(lx_m1,ly_m1,lz_m1);
    0218
    0219     end
    0220
    0221     Ind_dd=sub2ind(Asize,xx(:),yy(:),zz(:));
    0222
    0223     part_deriv = Input(IndBase)-Input(Ind_dd);
    0224
    0225     if n_dd >1
    0226         MatMinMax= (sign_Prev_deriv==sign(part_deriv)).*MatMinMax;
    0227     else
    0228         MatMinMax=sign(part_deriv);
    0229     end
    0230
    0231     sign_Prev_deriv=sign(part_deriv);
    0232 end
    0233
    0234 %Well , now the easy part, all values MatMinMax ==1 are local maximum and
    0235 %the values MatMinMax ==-1 are minimun
    0236
    0237 AllMaxima=find(MatMinMax==1);
    0238 AllMinima=find(MatMinMax==-1);
    0239
    0240 if numbermax ==0
    0241     nmax=length(AllMaxima);
    0242 else
    0243     nmax=numbermax;
    0244 end
    0245 nmax=min([nmax,length(AllMaxima)]);
    0246 smax=1:nmax;
    0247
    0248 if numbermin ==0
    0249     nmin=length(AllMinima);
    0250 else
    0251     nmin=numbermin;
    0252 end
    0253
    0254 nmin=min([nmin,length(AllMinima)]);
    0255
    0256 smin=1:nmin;
    0257
    0258 [Maxima,IndMax]=sort(Input(AllMaxima),'descend');
    0259 Maxima=Maxima(smax);
    0260 IndMax=AllMaxima(IndMax(smax));
    0261
    0262 MaxPos=zeros(nmax,3);
    0263 [MaxPos(:,1),MaxPos(:,2),MaxPos(:,3)]=ind2sub(Asize,IndMax);
    0264
    0265 [Minima,IndMin]=sort(Input(AllMinima));
    0266 Minima=Minima(smin);
    0267 IndMin=AllMinima(IndMin(smin));
    0268
    0269 MinPos=zeros(nmin,3);
    0270 [MinPos(:,1),MinPos(:,2),MinPos(:,3)]=ind2sub(Asize,IndMin);
    0271