dawn@0: /* 
dawn@0: * rsmooth.c
dawn@0: * smooth columnwise by convolution by a running square window
dawn@0: * multiple passes allowed (--> triangular, gaussian, etc.)
dawn@0: * 
dawn@0: * Alain de Cheveigné, CNRS Ircam, Dec 2001.
dawn@0: * (c) CNRS 2001.
dawn@0: */
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: 
dawn@0: #define X_IN	prhs[0]
dawn@0: #define SMOOTH_IN prhs[1]
dawn@0: #define NPASSES_IN	prhs[2]
dawn@0: #define CLIP_IN prhs[3]
dawn@0: 
dawn@0: /* Output Arguments */
dawn@0: 
dawn@0: #define Y_OUT plhs[0]
dawn@0: 
dawn@0: static void rsmooth(
dawn@0: 						double *xp,
dawn@0: 						int smooth,
dawn@0: 						int npasses,
dawn@0: 						double *yp,
dawn@0: 						int m,
dawn@0: 						int n,
dawn@0: 						int clip
dawn@0: 						)
dawn@0: {
dawn@0: 	int h,i,j,k,mm;
dawn@0: 	double *work1, *work2, *a, *b, *tmp, sum;
dawn@0: 		
dawn@0: 	mm = m + smooth + npasses*(smooth-1);		/* nrows of work buffer */
dawn@0: 	work1 = (double *) mxCalloc(mm*n,sizeof(double));
dawn@0: 	work2 = (double *) mxCalloc(mm*n,sizeof(double));
dawn@0: 	
dawn@0: 	CHECKIN(work1, (char *) work1 + sizeof(double)*mm*n);
dawn@0: 	CHECKIN(work2, (char *) work2 + sizeof(double)*mm*n);
dawn@0: 		
dawn@0: 	a = work1;
dawn@0: 	b = work2;
dawn@0: 		
dawn@0: 	/* transfer data to buffer, preceded with leading zeros */	 
dawn@0: 	for (k=0; k<n; k++) {
dawn@0: 		for (j=0;j<m;j++) {
dawn@0: 			SET(a[ j+smooth+mm*k ]) = GET(xp[ j+m*k ]);
dawn@0: 		}
dawn@0: 	}
dawn@0: 		
dawn@0: 	/* smooth repeatedly */
dawn@0: 	for (h=0;h<npasses;h++) {
dawn@0: 	  for (k=0; k<n; k++) {	/* for each column */
dawn@0: 			sum=0;
dawn@0: 			for (j=smooth;j<mm;j++) {	/* for each row (inner loop) */
dawn@0: 				sum += - GET(a[j-smooth+mm*k]) + GET(a[j+mm*k]);
dawn@0: 				SET(b[j+mm*k]) = sum/smooth;
dawn@0: 			}
dawn@0: 		}
dawn@0: 		tmp=a; a=b; b=tmp;	/* swap for next round */
dawn@0: 	}
dawn@0: 				
dawn@0: 	/* transfer to output matrix */
dawn@0: 	if (clip) {
dawn@0: 		for (j=0;j<m;j++) {
dawn@0: 			for (k=0; k<n; k++) {
dawn@0: 				 SET(yp[j+m*k]) = GET(a[ j+(int)(npasses*(smooth-1)/2)+smooth + mm*k]);			 
dawn@0: 			}
dawn@0: 		}
dawn@0: 	} else {
dawn@0: 		for (j=0;j<m+npasses*(smooth-1);j++) {
dawn@0: 			for (k=0; k<n; k++) {
dawn@0: 				 SET(yp[j+(m+npasses*(smooth-1))*k]) = GET(a[j + smooth + mm*k]);			 
dawn@0: 			}
dawn@0: 		}
dawn@0: 	}
dawn@0: 	CHECKOUT(work1);
dawn@0: 	CHECKOUT(work2);
dawn@0: 	
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, *yp;
dawn@0: 	int m,n,npasses,smooth,clip;
dawn@0: 	
dawn@0: 	/* Check for proper number of arguments */
dawn@0: 	if (nrhs > 4 || nrhs < 2) {
dawn@0: 		mexErrMsgTxt("RSMOOTH takes at least 2 and most 4 input arguments");
dawn@0: 	} else if (nlhs > 1) {
dawn@0: 		mexErrMsgTxt("RSMOOTH allows at most 1 output argument");
dawn@0: 	}
dawn@0: 	
dawn@0: 	/* Check type of input */
dawn@0: 	m = mxGetM(X_IN); /* rows */
dawn@0: 	n = mxGetN(X_IN); /* columns */ 
dawn@0: 	if (!mxIsNumeric(X_IN) || mxIsComplex(X_IN) || 
dawn@0: 			mxIsSparse(X_IN) || !mxIsDouble(X_IN) || m<=1) {
dawn@0: 		mexErrMsgTxt("RSMOOTH: X should be a column vector or matrix of doubles");
dawn@0: 	}
dawn@0:  
dawn@0: 	if (mxGetM(SMOOTH_IN)*mxGetN(SMOOTH_IN) > 1 || mxGetPr(SMOOTH_IN)[0] <1) {
dawn@0: 		mexErrMsgTxt("RSMOOTH: SMOOTH should be a positive scalar");
dawn@0: 	}	 
dawn@0: 	smooth = mxGetPr(SMOOTH_IN)[0];
dawn@0: 	
dawn@0: 	if (nrhs > 2) {
dawn@0: 		npasses = mxGetPr(NPASSES_IN)[0];
dawn@0: 		if (mxGetM(NPASSES_IN)*mxGetN(NPASSES_IN) > 1 || npasses <1) {
dawn@0: 			mexErrMsgTxt("RSMOOTH: NPASSES should be a positive scalar");
dawn@0: 		} 
dawn@0: 	} else {
dawn@0: 		npasses=1;
dawn@0: 	}
dawn@0: 	
dawn@0: 	if (nrhs > 3) {
dawn@0: 		clip=1; 
dawn@0: 	} else {
dawn@0: 		clip=0;
dawn@0: 	}
dawn@0: 	
dawn@0: 	/* Create matrix to return */
dawn@0: 	if (clip) {
dawn@0: 		Y_OUT = mxCreateDoubleMatrix(m, n, mxREAL); 
dawn@0: 	} else {
dawn@0: 		Y_OUT = mxCreateDoubleMatrix(m+npasses*(smooth-1), n, mxREAL); 
dawn@0: 	}
dawn@0: 		
dawn@0: 	/* Assign pointers to the various parameters */
dawn@0: 	xp = mxGetPr(X_IN);
dawn@0: 	yp = mxGetPr(Y_OUT);
dawn@0: 	
dawn@0: 	checkin_matrix((mxArray *) X_IN);
dawn@0: 	checkin_matrix(Y_OUT);
dawn@0: 	
dawn@0: 	/* Do the actual computations in a subroutine */
dawn@0: 	rsmooth(xp,smooth,npasses,yp, m,n,clip);
dawn@0:  
dawn@0: 	checkout_matrix((mxArray *) X_IN);
dawn@0: 	checkout_matrix(Y_OUT);
dawn@0: 	return;
dawn@0: }
dawn@0: 
dawn@0: