comparison toolboxes/FullBNT-1.0.7/bnt/potentials/Tables/rep_mult.c @ 0:e9a9cd732c1e tip

first hg version after svn

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