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)
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