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 | } |
---|