dawn@0: /* 
dawn@0: * rsum_inplace.c
dawn@0: * in-place running sum
dawn@0: * 
dawn@0: * Alain de Cheveigné, CNRS/Ircam	Jun 2002
dawn@0: * (c) 2002 CNRS
dawn@0: */
dawn@0: 
dawn@0: /* Replaces each sample of the input matrix by the sum of itself and its N-1 
dawn@0: * neighbors. The operation is done in-place on the input argument.  The last
dawn@0: * N-1 samples of each row are invalid. */
dawn@0: 
dawn@0: #include <math.h>
dawn@0: #include "mex.h"
dawn@0: 
dawn@0: /* #define MWRAP */
dawn@0: #include "mwrap_check.h"
dawn@0: 
dawn@0: /* Input Arguments */
dawn@0: #define X_IN	prhs[0]
dawn@0: #define N_IN	prhs[1]
dawn@0: 
dawn@0: /* Output Arguments */
dawn@0: 
dawn@0: static void rsmth(
dawn@0: 				double *xp,		/* matrix to sum */
dawn@0: 				double N,		/* window size */
dawn@0: 				int m, 			/* rows */ 
dawn@0: 				int n	 		/* columns */
dawn@0: 			) 
dawn@0: {
dawn@0: 	int j,k, Ni;
dawn@0: 	double Nr, tmp, sum, *x, *xx, *bump;
dawn@0: 	
dawn@0: 	Ni = (int) floor(N);
dawn@0: 	Nr = N - (double) Ni;
dawn@0: 	
dawn@0: 	
dawn@0: 	if (Nr == 0.0) {
dawn@0: 		/* N is integer: simple summation */
dawn@0: 		for (j=0;j<n;j++) {								/* for each column */
dawn@0: 			/* prefill running sum */
dawn@0: 			sum=0.0;
dawn@0: 			x=&xp[j*m];
dawn@0: 			bump=x+Ni-1;
dawn@0: 			while (x<bump) {				/* for each row */
dawn@0: 				sum += GET(*x);
dawn@0: 				x++;
dawn@0: 			}
dawn@0: 			/* N */
dawn@0: 			x=&xp[j*m];
dawn@0: 			bump=x+m-(Ni-1);
dawn@0: 			xx=x+Ni-1;
dawn@0: 			while(x<bump) {
dawn@0: 				tmp = GET(*x);
dawn@0: 				sum += GET(*xx);
dawn@0: 				SET(*x) = sum;
dawn@0: 				sum -= tmp;
dawn@0: 				x++; xx++;
dawn@0: 			}
dawn@0: 		}
dawn@0: 	} else {
dawn@0: 		/* N is fractionary: interpolate */
dawn@0: 		for (j=0;j<n;j++) {								/* for each column */
dawn@0: 			/* prefill running sum */
dawn@0: 			sum=0.0;
dawn@0: 			x=&xp[j*m];
dawn@0: 			bump=x+Ni-1;
dawn@0: 			while (x<bump) {				/* for each row */
dawn@0: 				sum += GET(*x);
dawn@0: 				x++;
dawn@0: 			}
dawn@0: 			/* N */
dawn@0: 			x=&xp[j*m];
dawn@0: 			bump=x+m-(Ni-1);
dawn@0: 			xx=x+Ni-1;
dawn@0: 			while(x<bump) {
dawn@0: 				tmp = GET(*x);
dawn@0: 				sum += GET(*xx);
dawn@0: 				SET(*x) = sum + Nr * GET(*xx);
dawn@0: 				sum -= tmp;
dawn@0: 				x++; xx++;
dawn@0: 			}
dawn@0: 		}
dawn@0: 	}	
dawn@0: 	return;						
dawn@0: }
dawn@0: 
dawn@0: void mexFunction(
dawn@0: 				int nlhs, mxArray *plhs[],
dawn@0: 				int nrhs, const mxArray *prhs[]
dawn@0: 				)
dawn@0: {
dawn@0: 	double	*xp, *Np, N; 
dawn@0: 	int n, m;
dawn@0: 	
dawn@0: 	/* Check for proper number of arguments */
dawn@0: 	if (nrhs != 2) {
dawn@0: 		mexErrMsgTxt("RSUM_INPLACE takes 2 input arguments");
dawn@0: 	}
dawn@0: 	
dawn@0: 	/* Check type of input */
dawn@0: 	if (!mxIsNumeric(X_IN) || mxIsComplex(X_IN) || 
dawn@0: 		mxIsSparse(X_IN)	|| !mxIsDouble(X_IN) ) {
dawn@0: 		mexErrMsgTxt("RSUM_INPLACE: X should be doubles");
dawn@0: 	}
dawn@0: 	if (!mxIsNumeric(N_IN) || mxIsComplex(N_IN) || 
dawn@0: 		mxIsSparse(N_IN)	|| !mxIsDouble(N_IN) ) {
dawn@0: 		mexErrMsgTxt("RSUM_INPLACE: N should be doubles");
dawn@0: 	}
dawn@0: 	m=mxGetM(X_IN);		/* rows */
dawn@0: 	n=mxGetN(X_IN);		/* columns */
dawn@0: 
dawn@0: 	if (mxGetM(N_IN)*mxGetN(N_IN) != 1) {
dawn@0: 		mexErrMsgTxt("RSUM_INPLACE: N should be a scalar");
dawn@0: 	}
dawn@0: 			
dawn@0: 	xp = mxGetPr(X_IN);
dawn@0: 	Np = mxGetPr(N_IN);
dawn@0: 	N = *Np;
dawn@0: 
dawn@0: 	if (m < ceil(N)) {
dawn@0: 		mexErrMsgTxt("RSUM_INPLACE: X should have at least N rows");
dawn@0: 	}
dawn@0: 		
dawn@0: 	checkin_matrix((mxArray *) X_IN);
dawn@0: 	
dawn@0: 	/* Do the actual computations in a subroutine */
dawn@0: 	rsmth(xp, N, m, n);
dawn@0: 
dawn@0: 	checkout_matrix((mxArray *) X_IN);
dawn@0: 	return;
dawn@0: }
dawn@0: