[37] | 1 | /* file: tpsumx.c |
---|
| 2 | ** author: Andrea Vedaldi |
---|
| 3 | ** description: Compute the matrix U for a TPS. |
---|
| 4 | **/ |
---|
| 5 | |
---|
| 6 | #include"mex.h" |
---|
| 7 | #include"math.h" |
---|
| 8 | #include<mexutils.c> |
---|
| 9 | #include<stdlib.h> |
---|
| 10 | |
---|
| 11 | /** Matlab driver. |
---|
| 12 | **/ |
---|
| 13 | #define greater(a,b) (a) > (b) |
---|
| 14 | #define max(a,b) (((a)>(b))?(a):(b)) |
---|
| 15 | #define min(a,b) (((a)<(b))?(a):(b)) |
---|
| 16 | #define getM(arg) mxGetM(in[arg]) |
---|
| 17 | #define getN(arg) mxGetN(in[arg]) |
---|
| 18 | #define getPr(arg) mxGetPr(in[arg]) |
---|
| 19 | |
---|
| 20 | void |
---|
| 21 | mexFunction(int nout, mxArray *out[], |
---|
| 22 | int nin, const mxArray *in[]) |
---|
| 23 | { |
---|
| 24 | enum { X=0,Y } ; |
---|
| 25 | enum { U } ; |
---|
| 26 | |
---|
| 27 | int MP, NP, MCP, NCP ; |
---|
| 28 | int i,j ; |
---|
| 29 | double *X_pt, *Y_pt, *U_pt ; |
---|
| 30 | const double small = 2.220446049250313e-16 ; |
---|
| 31 | |
---|
| 32 | /* ----------------------------------------------------------------- |
---|
| 33 | * Check the arguments |
---|
| 34 | * -------------------------------------------------------------- */ |
---|
| 35 | if(nin != 2) { |
---|
| 36 | mexErrMsgTxt("Two input arguments required") ; |
---|
| 37 | } |
---|
| 38 | |
---|
| 39 | if(!uIsRealMatrix(in[X], 2, -1)) { |
---|
| 40 | mexErrMsgTxt("X must be a 2xNP real matrix") ; |
---|
| 41 | } |
---|
| 42 | |
---|
| 43 | if(!uIsRealMatrix(in[Y], 2, -1)) { |
---|
| 44 | mexErrMsgTxt("Y must be a 2xNCP real matrix") ; |
---|
| 45 | } |
---|
| 46 | |
---|
| 47 | NP = getN(X) ; |
---|
| 48 | NCP = getN(Y) ; |
---|
| 49 | |
---|
| 50 | X_pt = getPr(X); |
---|
| 51 | Y_pt = getPr(Y) ; |
---|
| 52 | |
---|
| 53 | /* Allocate the result. */ |
---|
| 54 | out[U] = mxCreateDoubleMatrix(NP, NCP, mxREAL) ; |
---|
| 55 | U_pt = mxGetPr(out[U]) ; |
---|
| 56 | |
---|
| 57 | /* ----------------------------------------------------------------- |
---|
| 58 | * Check the arguments |
---|
| 59 | * -------------------------------------------------------------- */ |
---|
| 60 | for(j = 0 ; j < NCP ; ++j) { |
---|
| 61 | double xcp = *Y_pt++ ; |
---|
| 62 | double ycp = *Y_pt++ ; |
---|
| 63 | for(i = 0 ; i < NP ; ++i) { |
---|
| 64 | double dx = *X_pt++ - xcp ; |
---|
| 65 | double dy = *X_pt++ - ycp ; |
---|
| 66 | double r2 = dx*dx + dy*dy ; |
---|
| 67 | *U_pt++ = (r2 <= small) ? 0 : r2 * log( r2 ) ; |
---|
| 68 | } |
---|
| 69 | X_pt -= 2*NP ; |
---|
| 70 | } |
---|
| 71 | } |
---|