Code covered by the BSD License

# MinimaMaxima3D

### Samuel Pichardo (view profile)

13 Dec 2007 (Updated )

Find the minima and maxima in a 3D Cartesian data space

[Maxima,MaxPos,Minima,MinPos]=MinimaMaxima3D(Input,Robust,LookInBoundaries,numbermax,numbermin)
```function [Maxima,MaxPos,Minima,MinPos]=MinimaMaxima3D(Input,Robust,LookInBoundaries,numbermax,numbermin)
% 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)

if nargin <1
pf=test.uresTot(test.EvalLims(2,1):test.EvalLims(2,2));
pf=reshape(pf,length(test.EvalCoord{2}.Ry),length(test.EvalCoord{2}.Rx),length(test.EvalCoord{2}.Rz));
Input = abs(pf)*1.5e6;
clear test;
clear pf;
Robust =1;
end

Asize=size(Input);

if length(Asize)<3
error('MinimaMaxima3D can only works with 3D matrices ');
end

if (Asize(1)<3 || Asize(2)<3 || Asize(3)<3)
error('MinimaMaxima3D can only works with matrices with dimensions equal or larger to [3x3x3]');
end

if ~isreal(Input)
warning('ATTENTION, complex values detected!!, using abs(Input)');
Input=abs(Input);
end

if ~exist('Robust','var')
Robust=1;
end

if ~exist('LookInBoundaries','var')
LookInBoundaries=0;
end

if ~exist('numbermax','var')
numbermax=0;
end

if ~exist('numbermin','var')
numbermin=0;
end

[xx_base,yy_base,zz_base]=ndgrid(1:Asize(1),1:Asize(2),1:Asize(3));

IndBase=sub2ind(Asize,xx_base(:),yy_base(:),zz_base(:));

if Robust ~= 0
Numbder_dd=26;
else
Numbder_dd=6;
end

if LookInBoundaries==0
lx=1:Asize(1);
lx_p1=[2:Asize(1),Asize(1)];
lx_m1=[1,1:Asize(1)-1];
ly=1:Asize(2);
ly_p1=[2:Asize(2),Asize(2)];
ly_m1=[1,1:Asize(2)-1];
lz=1:Asize(3);
lz_p1=[2:Asize(3),Asize(3)];
lz_m1=[1,1:Asize(3)-1];
else
lx=1:Asize(1);
lx_p1=[2:Asize(1),Asize(1)-1]; %We replicate the layer N-1 as the layer N+1
lx_m1=[2,1:Asize(1)-1]; %We replicate the layer 2 as the layer -1
ly=1:Asize(2);
ly_p1=[2:Asize(2),Asize(2)-1]; %We replicate the layer N-1 as the layer N+1
ly_m1=[2,1:Asize(2)-1]; %We replicate the layer 2 as the layer -1
lz=1:Asize(3);
lz_p1=[2:Asize(3),Asize(3)-1]; %We replicate the layer N-1 as the layer N+1
lz_m1=[2,1:Asize(3)-1];%We replicate the layer 2 as the layer -1
end

for n_dd=1:Numbder_dd
switch n_dd
case 1
%%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x+1)
[xx,yy,zz]=ndgrid(lx_p1,ly,lz);

case 2
%%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x-1)
[xx,yy,zz]=ndgrid(lx_m1,ly,lz);

case 3
%%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(y)-elem(y+1)
[xx,yy,zz]=ndgrid(lx,ly_p1,lz);

case 4
%%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(y)-elem(y-1)
[xx,yy,zz]=ndgrid(lx,ly_m1,lz);

case 5
%%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(z)-elem(z+1)
[xx,yy,zz]=ndgrid(lx,ly,lz_p1);

case 6
%%%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(z)-elem(z-1)
[xx,yy,zz]=ndgrid(lx,ly,lz_m1);
case 7
%%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x+1,y+1)
[xx,yy,zz]=ndgrid(lx_p1,ly_p1,lz);
case 8
%%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x+1,y-1)
[xx,yy,zz]=ndgrid(lx_p1,ly_m1,lz);
case 9
%%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x-1,y-1)
[xx,yy,zz]=ndgrid(lx_m1,ly_m1,lz);
case 10
%%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x-1,y+1)
[xx,yy,zz]=ndgrid(lx_m1,ly_p1,lz);
case 11
%%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x+1,z+1)
[xx,yy,zz]=ndgrid(lx_p1,ly,lz_p1);
case 12
%%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x+1,z-1)
[xx,yy,zz]=ndgrid(lx_p1,ly,lz_m1);
case 13
%%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x-1,z-1)
[xx,yy,zz]=ndgrid(lx_m1,ly,lz_m1);
case 14
%%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x-1,z+1)
[xx,yy,zz]=ndgrid(lx_m1,ly,lz_p1);
case 15
%%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(y+1,z+1)
[xx,yy,zz]=ndgrid(lx,ly_p1,lz_p1);
case 16
%%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(y+1,z-1)
[xx,yy,zz]=ndgrid(lx,ly_p1,lz_m1);
case 17
%%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(y-1,z-1)
[xx,yy,zz]=ndgrid(lx,ly_m1,lz_m1);
case 18
%%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(y-1,z+1)
[xx,yy,zz]=ndgrid(lx,ly_m1,lz_p1);
case 19
%%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x+1,y+1,z+1)
[xx,yy,zz]=ndgrid(lx_p1,ly_p1,lz_p1);
case 20
%%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x+1,y+1,z-1)
[xx,yy,zz]=ndgrid(lx_p1,ly_p1,lz_m1);
case 21
%%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x+1,y-1,z+1)
[xx,yy,zz]=ndgrid(lx_p1,ly_m1,lz_p1);
case 22
%%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x+1,y-1,z-1)
[xx,yy,zz]=ndgrid(lx_p1,ly_m1,lz_m1);
case 23
%%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x-1,y+1,z+1)
[xx,yy,zz]=ndgrid(lx_m1,ly_p1,lz_p1);
case 24
%%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x-1,y+1,z-1)
[xx,yy,zz]=ndgrid(lx_m1,ly_p1,lz_m1);
case 25
%%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x-1,y-1,z+1)
[xx,yy,zz]=ndgrid(lx_m1,ly_m1,lz_p1);
case 26
%%%%%%%%%%%%%%%%%% %% This index is used to calculated elem(x)-elem(x-1,y-1,z-1)
[xx,yy,zz]=ndgrid(lx_m1,ly_m1,lz_m1);

end

Ind_dd=sub2ind(Asize,xx(:),yy(:),zz(:));

part_deriv = Input(IndBase)-Input(Ind_dd);

if n_dd >1
MatMinMax= (sign_Prev_deriv==sign(part_deriv)).*MatMinMax;
else
MatMinMax=sign(part_deriv);
end

sign_Prev_deriv=sign(part_deriv);
end

%Well , now the easy part, all values MatMinMax ==1 are local maximum and
%the values MatMinMax ==-1 are minimun

AllMaxima=find(MatMinMax==1);
AllMinima=find(MatMinMax==-1);

if numbermax ==0
nmax=length(AllMaxima);
else
nmax=numbermax;
end
nmax=min([nmax,length(AllMaxima)]);
smax=1:nmax;

if numbermin ==0
nmin=length(AllMinima);
else
nmin=numbermin;
end

nmin=min([nmin,length(AllMinima)]);

smin=1:nmin;

[Maxima,IndMax]=sort(Input(AllMaxima),'descend');
Maxima=Maxima(smax);
IndMax=AllMaxima(IndMax(smax));

MaxPos=zeros(nmax,3);
[MaxPos(:,1),MaxPos(:,2),MaxPos(:,3)]=ind2sub(Asize,IndMax);

[Minima,IndMin]=sort(Input(AllMinima));
Minima=Minima(smin);
IndMin=AllMinima(IndMin(smin));

MinPos=zeros(nmin,3);
[MinPos(:,1),MinPos(:,2),MinPos(:,3)]=ind2sub(Asize,IndMin);

```