From addbb92ff06b21d2188eeefdb17692beb85e4b1c Mon Sep 17 00:00:00 2001 From: Eugeniy Mikhailov Date: Tue, 22 Apr 2014 11:48:18 -0400 Subject: Hunter's initial release --- extrema2.m | 174 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 174 insertions(+) create mode 100644 extrema2.m (limited to 'extrema2.m') diff --git a/extrema2.m b/extrema2.m new file mode 100644 index 0000000..53a5e76 --- /dev/null +++ b/extrema2.m @@ -0,0 +1,174 @@ +function [xymax,smax,xymin,smin] = extrema2(xy,varargin) +%EXTREMA2 Gets the extrema points from a surface. +% [XMAX,IMAX,XMIN,IMIN] = EXTREMA2(X) returns the maxima and minima +% elements of the matriz X ignoring NaN's, where +% XMAX - maxima points in descending order (the bigger first and so on) +% IMAX - linear indexes of the XMAX +% XMIN - minima points in descending order +% IMIN - linear indexes of the XMIN. +% The program uses EXTREMA. +% +% The extrema points are searched only through the column, the row and +% the diagonals crossing each matrix element, so it is not a perfect +% mathematical program and for this reason it has an optional argument. +% The user should be aware of these limitations. +% +% [XMAX,IMAX,XMIN,IMIN] = EXTREMA2(X,1) does the same but without +% searching through the diagonals (less strict and perhaps the user gets +% more output points). +% +% DEFINITION (from http://en.wikipedia.org/wiki/Maxima_and_minima): +% In mathematics, maxima and minima, also known as extrema, are points in +% the domain of a function at which the function takes a largest value +% (maximum) or smallest value (minimum), either within a given +% neighbourhood (local extrema) or on the function domain in its entirety +% (global extrema). +% +% Note: To change the linear index to (i,j) use IND2SUB. +% +% Example: +% [x,y] = meshgrid(-2:.2:2,3:-.2:-2); +% z = x.*exp(-x.^2-y.^2); z(10,7)= NaN; z(16:19,13:17) = NaN; +% surf(x,y,z), shading interp +% [zmax,imax,zmin,imin] = extrema2(z); +% hold on +% plot3(x(imax),y(imax),zmax,'bo',x(imin),y(imin),zmin,'ro') +% for i = 1:length(zmax) +% text(x(imax(i)),y(imax(i)),zmax(i),[' ' num2str(zmax(i))]) +% end +% for i = 1:length(zmin) +% text(x(imin(i)),y(imin(i)),zmin(i),[' ' num2str(zmin(i))]) +% end +% hold off +% +% See also EXTREMA, MAX, MIN + +% Written by +% Lic. on Physics Carlos Adrián Vargas Aguilera +% Physical Oceanography MS candidate +% UNIVERSIDAD DE GUADALAJARA +% Mexico, 2005 +% +% nubeobscura@hotmail.com + +% From : http://www.mathworks.com/matlabcentral/fileexchange +% File ID : 12275 +% Submited at: 2006-09-14 +% 2006-11-11 : English translation from spanish. +% 2006-11-17 : Accept NaN's. +% 2006-11-22 : Fixed bug in INDX (by JaeKyu Suhr) +% 2007-04-09 : Change name to MAXIMA2, and definition added. + +M = size(xy); +if length(M) ~= 2 + error('Entry must be a matrix.') +end +N = M(2); +M = M(1); + +% Search peaks through columns: +[smaxcol,smincol] = extremos(xy); + +% Search peaks through rows, on columns with extrema points: +im = unique([smaxcol(:,1);smincol(:,1)]); % Rows with column extrema +[smaxfil,sminfil] = extremos(xy(im,:).'); + +% Convertion from 2 to 1 index: +smaxcol = sub2ind([M,N],smaxcol(:,1),smaxcol(:,2)); +smincol = sub2ind([M,N],smincol(:,1),smincol(:,2)); +smaxfil = sub2ind([M,N],im(smaxfil(:,2)),smaxfil(:,1)); +sminfil = sub2ind([M,N],im(sminfil(:,2)),sminfil(:,1)); + +% Peaks in rows and in columns: +smax = intersect(smaxcol,smaxfil); +smin = intersect(smincol,sminfil); + +% Search peaks through diagonals? +if nargin==1 + % Check peaks on down-up diagonal: + [iext,jext] = ind2sub([M,N],unique([smax;smin])); + [sextmax,sextmin] = extremos_diag(iext,jext,xy,1); + + % Check peaks on up-down diagonal: + smax = intersect(smax,[M; (N*M-M); sextmax]); + smin = intersect(smin,[M; (N*M-M); sextmin]); + + % Peaks on up-down diagonals: + [iext,jext] = ind2sub([M,N],unique([smax;smin])); + [sextmax,sextmin] = extremos_diag(iext,jext,xy,-1); + + % Peaks on columns, rows and diagonals: + smax = intersect(smax,[1; N*M; sextmax]); + smin = intersect(smin,[1; N*M; sextmin]); +end + +% Extrema points: +xymax = xy(smax); +xymin = xy(smin); + +% Descending order: +[temp,inmax] = sort(-xymax); clear temp +xymax = xymax(inmax); +smax = smax(inmax); +[xymin,inmin] = sort(xymin); +smin = smin(inmin); + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +function [smax,smin] = extremos(matriz) +% Peaks through columns or rows. + +smax = []; +smin = []; + +for n = 1:length(matriz(1,:)) + [temp,imaxfil,temp,iminfil] = extrema(matriz(:,n)); clear temp + if ~isempty(imaxfil) % Maxima indexes + imaxcol = repmat(n,length(imaxfil),1); + smax = [smax; imaxfil imaxcol]; + end + if ~isempty(iminfil) % Minima indexes + imincol = repmat(n,length(iminfil),1); + smin = [smin; iminfil imincol]; + end +end + + +function [sextmax,sextmin] = extremos_diag(iext,jext,xy,A) +% Peaks through diagonals (down-up A=-1) + +[M,N] = size(xy); +if A==-1 + iext = M-iext+1; +end +[iini,jini] = cruce(iext,jext,1,1); +[iini,jini] = ind2sub([M,N],unique(sub2ind([M,N],iini,jini))); +[ifin,jfin] = cruce(iini,jini,M,N); +sextmax = []; +sextmin = []; +for n = 1:length(iini) + ises = iini(n):ifin(n); + jses = jini(n):jfin(n); + if A==-1 + ises = M-ises+1; + end + s = sub2ind([M,N],ises,jses); + [temp,imax,temp,imin] = extrema(xy(s)); clear temp + sextmax = [sextmax; s(imax)']; + sextmin = [sextmin; s(imin)']; +end + + +function [i,j] = cruce(i0,j0,I,J) +% Indexes where the diagonal of the element io,jo crosses the left/superior +% (I=1,J=1) or right/inferior (I=M,J=N) side of an MxN matrix. + +arriba = 2*(I*J==1)-1; + +si = (arriba*(j0-J) > arriba*(i0-I)); +i = (I - (J+i0-j0)).*si + J+i0-j0; +j = (I+j0-i0-(J)).*si + J; + + +% Carlos Adrián Vargas Aguilera. nubeobscura@hotmail.com \ No newline at end of file -- cgit v1.2.3