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
% 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: email@example.com
% 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;