Daniel@0: /* rep_mult.c  repmat first two operands to the size provided by */
Daniel@0: /* the third operand, then perform point multiply                */
Daniel@0: /* 3 input, 1 output                                             */
Daniel@0: /* C = rep_mult(A, B, sizes)                                     */
Daniel@0: 
Daniel@0: #include "mex.h"
Daniel@0: 
Daniel@0: void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
Daniel@0: {
Daniel@0: 	double		*xp, *yp, *zp, *pSizes;
Daniel@0: 	int			xnd, ynd, numElements = 1;
Daniel@0: 	const int	*xdim, *ydim;
Daniel@0: 	int         i, j, ndim;
Daniel@0: 	int			*s, *sx, *sy, *cpsx, *cpsy;
Daniel@0: 	int			*subs, *s1, *cpsx2, *cpsy2;
Daniel@0: 
Daniel@0: 	if (nrhs != 3)
Daniel@0: 		mexErrMsgTxt("Incorrect number of inputs.");
Daniel@0: 	
Daniel@0: 	if (nlhs > 1)
Daniel@0: 		mexErrMsgTxt("Too many output arguments.");
Daniel@0: 	
Daniel@0: 	xnd = mxGetNumberOfDimensions(prhs[0]);
Daniel@0: 	ynd = mxGetNumberOfDimensions(prhs[1]);
Daniel@0: 	xdim = mxGetDimensions(prhs[0]);
Daniel@0: 	ydim = mxGetDimensions(prhs[1]);
Daniel@0: 	ndim = mxGetNumberOfElements(prhs[2]);
Daniel@0: 
Daniel@0: 	pSizes = mxGetPr(prhs[2]);
Daniel@0: 
Daniel@0: 	sx = (int *)malloc(sizeof(int)*ndim);
Daniel@0: 	sy = (int *)malloc(sizeof(int)*ndim);
Daniel@0: 	s =  (int *)malloc(sizeof(int)*ndim);
Daniel@0: 	s1 = (int *)malloc(sizeof(int)*ndim);
Daniel@0: 	*(cpsx = (int *)malloc(sizeof(int)*ndim)) = 1;
Daniel@0: 	*(cpsy = (int *)malloc(sizeof(int)*ndim)) = 1;
Daniel@0: 	subs =   (int *)malloc(sizeof(int)*ndim);
Daniel@0: 	cpsx2 =  (int *)malloc(sizeof(int)*ndim);
Daniel@0: 	cpsy2 =  (int *)malloc(sizeof(int)*ndim);
Daniel@0: 	for(i=0; i<ndim; i++){
Daniel@0: 		subs[i] = 0;
Daniel@0: 		sx[i] = (i < xnd) ? xdim[i] : 1;
Daniel@0: 		sy[i] = (i < ynd) ? ydim[i] : 1;
Daniel@0: 		s[i] = (int)pSizes[i];
Daniel@0: 		s1[i] = s[i] - 1;
Daniel@0: 		numElements *= s[i];
Daniel@0: 	}
Daniel@0: 				
Daniel@0: 	for(i=0; i<ndim-1; i++){
Daniel@0: 		cpsx[i+1] = cpsx[i]*sx[i]--;
Daniel@0: 		cpsy[i+1] = cpsy[i]*sy[i]--;
Daniel@0: 		cpsx2[i] = cpsx[i]*sx[i];
Daniel@0: 		cpsy2[i] = cpsy[i]*sy[i];
Daniel@0: 	}
Daniel@0: 	cpsx2[ndim-1] = cpsx[ndim-1]*(--sx[ndim-1]);
Daniel@0: 	cpsy2[ndim-1] = cpsy[ndim-1]*(--sy[ndim-1]);
Daniel@0: 	
Daniel@0: 	plhs[0] = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxREAL);
Daniel@0: 	zp = mxGetPr(plhs[0]);
Daniel@0: 	xp = mxGetPr(prhs[0]);
Daniel@0: 	yp = mxGetPr(prhs[1]);
Daniel@0: 
Daniel@0: 	for(j=0; j<numElements; j++){
Daniel@0: 		*zp++ = *xp * *yp;
Daniel@0: 		for(i=0; i<ndim; i++){
Daniel@0: 			if(subs[i] == s1[i]){
Daniel@0: 				subs[i] = 0;
Daniel@0: 				if(sx[i])
Daniel@0: 					xp -= cpsx2[i];
Daniel@0: 				if(sy[i])
Daniel@0: 					yp -= cpsy2[i];
Daniel@0: 			}
Daniel@0: 			else{
Daniel@0: 				subs[i]++;
Daniel@0: 				if(sx[i])
Daniel@0: 					xp += cpsx[i];
Daniel@0: 				if(sy[i])
Daniel@0: 					yp += cpsy[i];
Daniel@0: 				break;
Daniel@0: 			}
Daniel@0: 		}
Daniel@0: 	}
Daniel@0: 	free(sx);
Daniel@0: 	free(sy);
Daniel@0: 	free(s);
Daniel@0: 	free(s1);
Daniel@0: 	free(cpsx);
Daniel@0: 	free(cpsy);
Daniel@0: 	free(subs);
Daniel@0: 	free(cpsx2);
Daniel@0: 	free(cpsy2);
Daniel@0: }