1 | /* file: imwbackward.c |
---|
2 | ** author: Andrea Vedaldi |
---|
3 | ** description: Backward projection of an image. |
---|
4 | **/ |
---|
5 | |
---|
6 | /* TODO. |
---|
7 | * - Make a faster version for the uniform grid case. |
---|
8 | * - Rename variables. |
---|
9 | * - Check potential bug below!!!! |
---|
10 | */ |
---|
11 | #include"mex.h" |
---|
12 | #include"math.h" |
---|
13 | #include<mexutils.c> |
---|
14 | #include<stdlib.h> |
---|
15 | |
---|
16 | /** Matlab driver. |
---|
17 | **/ |
---|
18 | #define greater(a,b) (a) > (b) |
---|
19 | #define max(a,b) (((a)>(b))?(a):(b)) |
---|
20 | #define min(a,b) (((a)<(b))?(a):(b)) |
---|
21 | #define getM(arg) mxGetM(in[arg]) |
---|
22 | #define getN(arg) mxGetN(in[arg]) |
---|
23 | #define getPr(arg) mxGetPr(in[arg]) |
---|
24 | |
---|
25 | int |
---|
26 | findNeighbor(double x, const double* X, int K) { |
---|
27 | int i = 0 ; |
---|
28 | int j = K - 1 ; |
---|
29 | int pivot = 0 ; |
---|
30 | double y = 0 ; |
---|
31 | if(x < X[i]) return i-1 ; |
---|
32 | if(x >= X[j]) return j ; |
---|
33 | |
---|
34 | while(i < j - 1) { |
---|
35 | pivot = (i+j) >> 1 ; |
---|
36 | y = X[pivot] ; |
---|
37 | /*mexPrintf("%d %d %d %f %f\n",i,j,pivot,x,y) ;*/ |
---|
38 | if(x < y) { |
---|
39 | j = pivot ; |
---|
40 | } else { |
---|
41 | i = pivot ; |
---|
42 | } |
---|
43 | } |
---|
44 | return i ; |
---|
45 | } |
---|
46 | |
---|
47 | void |
---|
48 | mexFunction(int nout, mxArray *out[], |
---|
49 | int nin, const mxArray *in[]) |
---|
50 | { |
---|
51 | enum { X=0,Y,I,iwXp,iwYp } ; |
---|
52 | enum { wI=0 } ; |
---|
53 | |
---|
54 | int M, N, Mp, Np, ip, jp ; |
---|
55 | double *X_pt, *Y_pt, *I_pt, *iwXp_pt, *iwYp_pt, *wI_pt ; |
---|
56 | double Xmin, Xmax, Ymin, Ymax ; |
---|
57 | const double NaN = mxGetNaN() ; |
---|
58 | |
---|
59 | /* ----------------------------------------------------------------- |
---|
60 | * Check the arguments |
---|
61 | * -------------------------------------------------------------- */ |
---|
62 | if(nin != 5) { |
---|
63 | mexErrMsgTxt("Five input argumets required.") ; |
---|
64 | } |
---|
65 | |
---|
66 | if(!uIsRealMatrix(in[I], -1, -1)) { |
---|
67 | mexErrMsgTxt("I must be a real matrix of class DOUBLE") ; |
---|
68 | } |
---|
69 | |
---|
70 | if(!uIsRealMatrix(in[iwXp], -1, -1)) { |
---|
71 | mexErrMsgTxt("iwXp must be a real matrix") ; |
---|
72 | } |
---|
73 | |
---|
74 | M = getM(I) ; |
---|
75 | N = getN(I) ; |
---|
76 | Mp = getM(iwXp) ; |
---|
77 | Np = getN(iwXp) ; |
---|
78 | |
---|
79 | if(!uIsRealMatrix(in[iwYp], Mp, Np)) { |
---|
80 | mexErrMsgTxt |
---|
81 | ("iwXp and iwYp must be a real matrices of the same dimension") ; |
---|
82 | } |
---|
83 | |
---|
84 | if(!uIsRealVector(in[X],N) || !uIsRealVector(in[Y],M)) { |
---|
85 | mexErrMsgTxt |
---|
86 | ("X and Y must be vectors of the same dimensions " |
---|
87 | "of the columns/rows of I, respectivelye") ; |
---|
88 | } |
---|
89 | |
---|
90 | X_pt = getPr(X); |
---|
91 | Y_pt = getPr(Y) ; |
---|
92 | I_pt = getPr(I) ; |
---|
93 | iwXp_pt = getPr(iwXp) ; |
---|
94 | iwYp_pt = getPr(iwYp) ; |
---|
95 | |
---|
96 | /* Allocate the result. */ |
---|
97 | out[wI] = mxCreateDoubleMatrix(Mp, Np, mxREAL) ; |
---|
98 | wI_pt = mxGetPr(out[wI]) ; |
---|
99 | |
---|
100 | /* ----------------------------------------------------------------- |
---|
101 | * Do the job |
---|
102 | * -------------------------------------------------------------- */ |
---|
103 | Xmin = X_pt[0] ; |
---|
104 | Xmax = X_pt[N-1] ; |
---|
105 | Ymin = Y_pt[0] ; |
---|
106 | Ymax = Y_pt[M-1] ; |
---|
107 | |
---|
108 | for(jp = 0 ; jp < Np ; ++jp) { |
---|
109 | for(ip = 0 ; ip < Mp ; ++ip) { |
---|
110 | /* Search for the four neighbors of the backprojected point. */ |
---|
111 | double x = *iwXp_pt++ ; |
---|
112 | double y = *iwYp_pt++ ; |
---|
113 | double z = NaN ; |
---|
114 | |
---|
115 | /* This messy code allows the identity transformation |
---|
116 | * to be processed as expected. */ |
---|
117 | if(x >= Xmin && x <= Xmax && |
---|
118 | y >= Ymin && y <= Ymax) { |
---|
119 | int j = findNeighbor(x, X_pt, N) ; |
---|
120 | int i = findNeighbor(y, Y_pt, M) ; |
---|
121 | double* pt = I_pt + j*M + i ; |
---|
122 | |
---|
123 | /* Weights. */ |
---|
124 | double x0 = X_pt[j] ; |
---|
125 | double x1 = X_pt[j+1] ; |
---|
126 | double y0 = Y_pt[i] ; |
---|
127 | double y1 = Y_pt[i+1] ; |
---|
128 | double wx = (x-x0)/(x1-x0) ; |
---|
129 | double wy = (y-y0)/(y1-y0) ; |
---|
130 | |
---|
131 | /* Load all possible neighbors. */ |
---|
132 | double z00 = 0.0 ; |
---|
133 | double z10 = 0.0 ; |
---|
134 | double z01 = 0.0 ; |
---|
135 | double z11 = 0.0 ; |
---|
136 | |
---|
137 | /* 2006: possible bug? should pt++ in any case? */ |
---|
138 | if(j > -1) { |
---|
139 | if(i > -1 ) z00 = *pt ; |
---|
140 | pt++ ; |
---|
141 | if(i < M-1) z10 = *pt ; |
---|
142 | } |
---|
143 | |
---|
144 | pt += M - 1; |
---|
145 | |
---|
146 | if(j < N-1) { |
---|
147 | if(i > -1 ) z01 = *pt ; |
---|
148 | pt++ ; |
---|
149 | if(i < M-1) z11 = *pt ; |
---|
150 | } |
---|
151 | |
---|
152 | /* Bilinear interpolation. */ |
---|
153 | z = |
---|
154 | (1-wy)*( (1-wx) * z00 + wx * z01) + |
---|
155 | wy*( (1-wx) * z10 + wx * z11) ; |
---|
156 | } |
---|
157 | |
---|
158 | *(wI_pt + jp*Mp + ip) = z ; |
---|
159 | } |
---|
160 | } |
---|
161 | } |
---|