[37] | 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 | } |
---|