[37] | 1 | /* C implementation of repmat. |
---|
| 2 | * by Tom Minka |
---|
| 3 | */ |
---|
| 4 | /* |
---|
| 5 | mex -c mexutil.c |
---|
| 6 | mex repmat.c mexutil.obj |
---|
| 7 | to check for warnings: |
---|
| 8 | gcc -Wall -I/cygdrive/c/MATLAB6p1/extern/include -c repmat.c |
---|
| 9 | */ |
---|
| 10 | #include "mexutil.h" |
---|
| 11 | #include <string.h> |
---|
| 12 | |
---|
| 13 | /* ALWAYS_2D == 1 gives the usual matlab behavior where repmat(0,3) is |
---|
| 14 | * a 3x3 matrix. |
---|
| 15 | * ALWAYS_2D == 0 gives a 3x1 matrix instead. |
---|
| 16 | * repmat(x,sizes,1) is another way to get this behavior. |
---|
| 17 | */ |
---|
| 18 | #define ALWAYS_2D 1 |
---|
| 19 | |
---|
| 20 | /* repeat a block of memory rep times */ |
---|
| 21 | void memrep(char *dest, size_t chunk, int rep) |
---|
| 22 | { |
---|
| 23 | #if 0 |
---|
| 24 | /* slow way */ |
---|
| 25 | int i; |
---|
| 26 | char *p = dest; |
---|
| 27 | for(i=1;i<rep;i++) { |
---|
| 28 | p += chunk; |
---|
| 29 | memcpy(p, dest, chunk); |
---|
| 30 | } |
---|
| 31 | #else |
---|
| 32 | /* fast way */ |
---|
| 33 | if(rep == 1) return; |
---|
| 34 | memcpy(dest + chunk, dest, chunk); |
---|
| 35 | if(rep & 1) { |
---|
| 36 | dest += chunk; |
---|
| 37 | memcpy(dest + chunk, dest, chunk); |
---|
| 38 | } |
---|
| 39 | /* now repeat using a block twice as big */ |
---|
| 40 | memrep(dest, chunk<<1, rep>>1); |
---|
| 41 | #endif |
---|
| 42 | } |
---|
| 43 | |
---|
| 44 | void repmat(char *dest, const char *src, int ndim, int *destdimsize, |
---|
| 45 | int *dimsize, const int *dims, int *rep) |
---|
| 46 | { |
---|
| 47 | int d = ndim-1; |
---|
| 48 | int i, chunk; |
---|
| 49 | /* copy the first repetition into dest */ |
---|
| 50 | if(d == 0) { |
---|
| 51 | chunk = dimsize[0]; |
---|
| 52 | memcpy(dest,src,chunk); |
---|
| 53 | } |
---|
| 54 | else { |
---|
| 55 | /* recursively repeat each slice of src */ |
---|
| 56 | for(i=0;i<dims[d];i++) { |
---|
| 57 | repmat(dest + i*destdimsize[d-1], src + i*dimsize[d-1], |
---|
| 58 | ndim-1, destdimsize, dimsize, dims, rep); |
---|
| 59 | } |
---|
| 60 | chunk = destdimsize[d-1]*dims[d]; |
---|
| 61 | } |
---|
| 62 | /* copy the result rep-1 times */ |
---|
| 63 | memrep(dest,chunk,rep[d]); |
---|
| 64 | } |
---|
| 65 | |
---|
| 66 | void mexFunction(int nlhs, mxArray *plhs[], |
---|
| 67 | int nrhs, const mxArray *prhs[]) |
---|
| 68 | { |
---|
| 69 | const mxArray *srcmat; |
---|
| 70 | int ndim, *dimsize, eltsize; |
---|
| 71 | const int *dims; |
---|
| 72 | int ndimdest, *destdims, *destdimsize; |
---|
| 73 | char *src, *dest; |
---|
| 74 | int *rep; |
---|
| 75 | int i,nrep; |
---|
| 76 | int extra_rep = 1; |
---|
| 77 | int empty; |
---|
| 78 | |
---|
| 79 | if(nrhs < 2) mexErrMsgTxt("Usage: repmat(A, [M N ...])"); |
---|
| 80 | srcmat = prhs[0]; |
---|
| 81 | if(mxIsSparse(srcmat) || mxIsCell(srcmat) || mxIsStruct(srcmat)) { |
---|
| 82 | /* call Matlab's repmat */ |
---|
| 83 | mexCallMATLAB(nlhs,plhs,nrhs,prhs,"xrepmat");return; |
---|
| 84 | /* mexErrMsgTxt("Sorry, can't handle sparse matrices yet."); */ |
---|
| 85 | } |
---|
| 86 | ndim = mxGetNumberOfDimensions(srcmat); |
---|
| 87 | dims = mxGetDimensions(srcmat); |
---|
| 88 | eltsize = mxGetElementSize(srcmat); |
---|
| 89 | |
---|
| 90 | /* compute dimension sizes */ |
---|
| 91 | dimsize = mxCalloc(ndim, sizeof(int)); |
---|
| 92 | dimsize[0] = eltsize*dims[0]; |
---|
| 93 | for(i=1;i<ndim;i++) dimsize[i] = dimsize[i-1]*dims[i]; |
---|
| 94 | |
---|
| 95 | /* determine repetition vector */ |
---|
| 96 | ndimdest = ndim; |
---|
| 97 | if(nrhs == 2) { |
---|
| 98 | /* prhs[1] is a vector of reps */ |
---|
| 99 | nrep = mxGetN(prhs[1]); |
---|
| 100 | if(nrep > ndimdest) ndimdest = nrep; |
---|
| 101 | rep = mxCalloc(ndimdest, sizeof(int)); |
---|
| 102 | for(i=0;i<nrep;i++) { |
---|
| 103 | double repv = mxGetPr(prhs[1])[i]; |
---|
| 104 | rep[i] = (int)repv; |
---|
| 105 | } |
---|
| 106 | #if ALWAYS_2D |
---|
| 107 | if(nrep == 1) { |
---|
| 108 | /* special behavior */ |
---|
| 109 | nrep = 2; |
---|
| 110 | rep[1] = rep[0]; |
---|
| 111 | } |
---|
| 112 | #endif |
---|
| 113 | } |
---|
| 114 | else { |
---|
| 115 | /* concatenate all prhs's */ |
---|
| 116 | int ri=0; |
---|
| 117 | nrep = 0; |
---|
| 118 | for(i=0;i<nrhs-1;i++) { |
---|
| 119 | nrep += mxGetNumberOfElements(prhs[i+1]); |
---|
| 120 | } |
---|
| 121 | if(nrep > ndimdest) ndimdest = nrep; |
---|
| 122 | rep = mxCalloc(ndimdest, sizeof(int)); |
---|
| 123 | for(i=0;i<nrhs-1;i++) { |
---|
| 124 | double *p = mxGetPr(prhs[i+1]); |
---|
| 125 | int j, sz = mxGetNumberOfElements(prhs[i+1]); |
---|
| 126 | for(j=0;j<sz;j++) rep[ri++] = (int)p[j]; |
---|
| 127 | } |
---|
| 128 | } |
---|
| 129 | for(i=nrep;i<ndimdest;i++) rep[i] = 1; |
---|
| 130 | |
---|
| 131 | /* compute output size */ |
---|
| 132 | destdims = mxCalloc(ndimdest, sizeof(int)); |
---|
| 133 | for(i=0;i<ndim;i++) destdims[i] = dims[i]*rep[i]; |
---|
| 134 | for(;i<ndimdest;i++) { |
---|
| 135 | destdims[i] = rep[i]; |
---|
| 136 | extra_rep *= rep[i]; |
---|
| 137 | } |
---|
| 138 | destdimsize = mxCalloc(ndim, sizeof(int)); |
---|
| 139 | destdimsize[0] = eltsize*destdims[0]; |
---|
| 140 | for(i=1;i<ndim;i++) destdimsize[i] = destdimsize[i-1]*destdims[i]; |
---|
| 141 | |
---|
| 142 | /* for speed, array should be uninitialized */ |
---|
| 143 | plhs[0] = mxCreateNumericArrayE(ndimdest, destdims, mxGetClassID(srcmat), |
---|
| 144 | mxIsComplex(srcmat)?mxCOMPLEX:mxREAL); |
---|
| 145 | |
---|
| 146 | /* if any rep[i] == 0, output should be empty array. |
---|
| 147 | Added by KPM 11/13/02. |
---|
| 148 | */ |
---|
| 149 | empty = 0; |
---|
| 150 | for (i=0; i < nrep; i++) { |
---|
| 151 | if (rep[i]==0) |
---|
| 152 | empty = 1; |
---|
| 153 | } |
---|
| 154 | if (empty) |
---|
| 155 | return; |
---|
| 156 | |
---|
| 157 | src = (char*)mxGetData(srcmat); |
---|
| 158 | dest = (char*)mxGetData(plhs[0]); |
---|
| 159 | repmat(dest,src,ndim,destdimsize,dimsize,dims,rep); |
---|
| 160 | if(ndimdest > ndim) memrep(dest,destdimsize[ndim-1],extra_rep); |
---|
| 161 | if(mxIsComplex(srcmat)) { |
---|
| 162 | src = (char*)mxGetPi(srcmat); |
---|
| 163 | dest = (char*)mxGetPi(plhs[0]); |
---|
| 164 | repmat(dest,src,ndim,destdimsize,dimsize,dims,rep); |
---|
| 165 | if(ndimdest > ndim) memrep(dest,destdimsize[ndim-1],extra_rep); |
---|
| 166 | } |
---|
| 167 | } |
---|