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