|
wolffd@0
|
1 /* rep_mult.c repmat first two operands to the size provided by */
|
|
wolffd@0
|
2 /* the third operand, then perform point multiply */
|
|
wolffd@0
|
3 /* 3 input, 1 output */
|
|
wolffd@0
|
4 /* C = rep_mult(A, B, sizes) */
|
|
wolffd@0
|
5
|
|
wolffd@0
|
6 #include "mex.h"
|
|
wolffd@0
|
7
|
|
wolffd@0
|
8 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
|
|
wolffd@0
|
9 {
|
|
wolffd@0
|
10 double *xp, *yp, *zp, *pSizes;
|
|
wolffd@0
|
11 int xnd, ynd, numElements = 1;
|
|
wolffd@0
|
12 const int *xdim, *ydim;
|
|
wolffd@0
|
13 int i, j, ndim;
|
|
wolffd@0
|
14 int *s, *sx, *sy, *cpsx, *cpsy;
|
|
wolffd@0
|
15 int *subs, *s1, *cpsx2, *cpsy2;
|
|
wolffd@0
|
16
|
|
wolffd@0
|
17 if (nrhs != 3)
|
|
wolffd@0
|
18 mexErrMsgTxt("Incorrect number of inputs.");
|
|
wolffd@0
|
19
|
|
wolffd@0
|
20 if (nlhs > 1)
|
|
wolffd@0
|
21 mexErrMsgTxt("Too many output arguments.");
|
|
wolffd@0
|
22
|
|
wolffd@0
|
23 xnd = mxGetNumberOfDimensions(prhs[0]);
|
|
wolffd@0
|
24 ynd = mxGetNumberOfDimensions(prhs[1]);
|
|
wolffd@0
|
25 xdim = mxGetDimensions(prhs[0]);
|
|
wolffd@0
|
26 ydim = mxGetDimensions(prhs[1]);
|
|
wolffd@0
|
27 ndim = mxGetNumberOfElements(prhs[2]);
|
|
wolffd@0
|
28
|
|
wolffd@0
|
29 pSizes = mxGetPr(prhs[2]);
|
|
wolffd@0
|
30
|
|
wolffd@0
|
31 sx = (int *)malloc(sizeof(int)*ndim);
|
|
wolffd@0
|
32 sy = (int *)malloc(sizeof(int)*ndim);
|
|
wolffd@0
|
33 s = (int *)malloc(sizeof(int)*ndim);
|
|
wolffd@0
|
34 s1 = (int *)malloc(sizeof(int)*ndim);
|
|
wolffd@0
|
35 *(cpsx = (int *)malloc(sizeof(int)*ndim)) = 1;
|
|
wolffd@0
|
36 *(cpsy = (int *)malloc(sizeof(int)*ndim)) = 1;
|
|
wolffd@0
|
37 subs = (int *)malloc(sizeof(int)*ndim);
|
|
wolffd@0
|
38 cpsx2 = (int *)malloc(sizeof(int)*ndim);
|
|
wolffd@0
|
39 cpsy2 = (int *)malloc(sizeof(int)*ndim);
|
|
wolffd@0
|
40 for(i=0; i<ndim; i++){
|
|
wolffd@0
|
41 subs[i] = 0;
|
|
wolffd@0
|
42 sx[i] = (i < xnd) ? xdim[i] : 1;
|
|
wolffd@0
|
43 sy[i] = (i < ynd) ? ydim[i] : 1;
|
|
wolffd@0
|
44 s[i] = (int)pSizes[i];
|
|
wolffd@0
|
45 s1[i] = s[i] - 1;
|
|
wolffd@0
|
46 numElements *= s[i];
|
|
wolffd@0
|
47 }
|
|
wolffd@0
|
48
|
|
wolffd@0
|
49 for(i=0; i<ndim-1; i++){
|
|
wolffd@0
|
50 cpsx[i+1] = cpsx[i]*sx[i]--;
|
|
wolffd@0
|
51 cpsy[i+1] = cpsy[i]*sy[i]--;
|
|
wolffd@0
|
52 cpsx2[i] = cpsx[i]*sx[i];
|
|
wolffd@0
|
53 cpsy2[i] = cpsy[i]*sy[i];
|
|
wolffd@0
|
54 }
|
|
wolffd@0
|
55 cpsx2[ndim-1] = cpsx[ndim-1]*(--sx[ndim-1]);
|
|
wolffd@0
|
56 cpsy2[ndim-1] = cpsy[ndim-1]*(--sy[ndim-1]);
|
|
wolffd@0
|
57
|
|
wolffd@0
|
58 plhs[0] = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxREAL);
|
|
wolffd@0
|
59 zp = mxGetPr(plhs[0]);
|
|
wolffd@0
|
60 xp = mxGetPr(prhs[0]);
|
|
wolffd@0
|
61 yp = mxGetPr(prhs[1]);
|
|
wolffd@0
|
62
|
|
wolffd@0
|
63 for(j=0; j<numElements; j++){
|
|
wolffd@0
|
64 *zp++ = *xp * *yp;
|
|
wolffd@0
|
65 for(i=0; i<ndim; i++){
|
|
wolffd@0
|
66 if(subs[i] == s1[i]){
|
|
wolffd@0
|
67 subs[i] = 0;
|
|
wolffd@0
|
68 if(sx[i])
|
|
wolffd@0
|
69 xp -= cpsx2[i];
|
|
wolffd@0
|
70 if(sy[i])
|
|
wolffd@0
|
71 yp -= cpsy2[i];
|
|
wolffd@0
|
72 }
|
|
wolffd@0
|
73 else{
|
|
wolffd@0
|
74 subs[i]++;
|
|
wolffd@0
|
75 if(sx[i])
|
|
wolffd@0
|
76 xp += cpsx[i];
|
|
wolffd@0
|
77 if(sy[i])
|
|
wolffd@0
|
78 yp += cpsy[i];
|
|
wolffd@0
|
79 break;
|
|
wolffd@0
|
80 }
|
|
wolffd@0
|
81 }
|
|
wolffd@0
|
82 }
|
|
wolffd@0
|
83 free(sx);
|
|
wolffd@0
|
84 free(sy);
|
|
wolffd@0
|
85 free(s);
|
|
wolffd@0
|
86 free(s1);
|
|
wolffd@0
|
87 free(cpsx);
|
|
wolffd@0
|
88 free(cpsy);
|
|
wolffd@0
|
89 free(subs);
|
|
wolffd@0
|
90 free(cpsx2);
|
|
wolffd@0
|
91 free(cpsy2);
|
|
wolffd@0
|
92 }
|
