12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970 |
- function matrixOut = smooth2a(matrixIn,Nr,Nc)
- % Smooths 2D array data. Ignores NaN's.
- %
- %function matrixOut = smooth2a(matrixIn,Nr,Nc)
- %
- % This function smooths the data in matrixIn using a mean filter over a
- % rectangle of size (2*Nr+1)-by-(2*Nc+1). Basically, you end up replacing
- % element "i" by the mean of the rectange centered on "i". Any NaN
- % elements are ignored in the averaging. If element "i" is a NaN, then it
- % will be preserved as NaN in the output. At the edges of the matrix,
- % where you cannot build a full rectangle, as much of the rectangle that
- % fits on your matrix is used (similar to the default on Matlab's builtin
- % function "smooth").
- %
- % "matrixIn": original matrix
- % "Nr": number of points used to smooth rows
- % "Nc": number of points to smooth columns. If not specified, Nc = Nr.
- %
- % "matrixOut": smoothed version of original matrix
- %
- %
- % Written by Greg Reeves, March 2009.
- % Division of Biology
- % Caltech
- %
- % Inspired by "smooth2", written by Kelly Hilands, October 2004
- % Applied Research Laboratory
- % Penn State University
- %
- % Developed from code written by Olof Liungman, 1997
- % Dept. of Oceanography, Earth Sciences Centre
- % G?teborg University, Sweden
- % E-mail: olof.liungman@oce.gu.se
- %
- % Initial error statements and definitions
- %
- if nargin < 2, error('Not enough input arguments!'), end
- N(1) = Nr;
- if nargin < 3, N(2) = N(1); else N(2) = Nc; end
- if length(N(1)) ~= 1, error('Nr must be a scalar!'), end
- if length(N(2)) ~= 1, error('Nc must be a scalar!'), end
- %
- % Building matrices that will compute running sums. The left-matrix, eL,
- % smooths along the rows. The right-matrix, eR, smooths along the
- % columns. You end up replacing element "i" by the mean of a (2*Nr+1)-by-
- % (2*Nc+1) rectangle centered on element "i".
- %
- [row,col] = size(matrixIn);
- eL = spdiags(ones(row,2*N(1)+1),(-N(1):N(1)),row,row);
- eR = spdiags(ones(col,2*N(2)+1),(-N(2):N(2)),col,col);
- %
- % Setting all "NaN" elements of "matrixIn" to zero so that these will not
- % affect the summation. (If this isn't done, any sum that includes a NaN
- % will also become NaN.)
- %
- A = isnan(matrixIn);
- matrixIn(A) = 0;
- %
- % For each element, we have to count how many non-NaN elements went into
- % the sums. This is so we can divide by that number to get a mean. We use
- % the same matrices to do this (ie, "eL" and "eR").
- %
- nrmlize = eL*(~A)*eR;
- nrmlize(A) = NaN;
- %
- % Actually taking the mean.
- %
- matrixOut = eL*matrixIn*eR;
- matrixOut = matrixOut./nrmlize;
|