[59] | 1 | ;+ |
---|
| 2 | ; |
---|
[125] | 3 | ; @file_comments |
---|
| 4 | ; warm (or map) a unit square onto an arbitrary quadrilateral |
---|
[59] | 5 | ; according to the 4-point correspondences: |
---|
| 6 | ; (0,0) -> (x0,y0) |
---|
| 7 | ; (1,0) -> (x1,y1) |
---|
| 8 | ; (1,1) -> (x2,y2) |
---|
| 9 | ; (0,1) -> (x3,y3) |
---|
| 10 | ; The mapping is done using perspective transformation which preserve |
---|
| 11 | ; lines in all orientations and permit quadrilateral to quadrilateral |
---|
| 12 | ; mappings. see ref. bellow. |
---|
| 13 | ; |
---|
[101] | 14 | ; @categories image, grid manipulation |
---|
[59] | 15 | ; |
---|
[125] | 16 | ; @param x0in {in}{required} the coordinates of the quadrilateral |
---|
[59] | 17 | ; (see above for correspondance with the unit square). Can be |
---|
| 18 | ; scalar or array. (x0,y0), (x1,y1), (x2,y2) and (x3,y3) are |
---|
| 19 | ; given in the anticlockwise order. |
---|
[125] | 20 | ; @param y0in {in}{required} the coordinates of the quadrilateral |
---|
[101] | 21 | ; (see above for correspondance with the unit square). Can be |
---|
| 22 | ; scalar or array. (x0,y0), (x1,y1), (x2,y2) and (x3,y3) are |
---|
| 23 | ; given in the anticlockwise order. |
---|
[125] | 24 | ; @param x1in {in}{required} the coordinates of the quadrilateral |
---|
[101] | 25 | ; (see above for correspondance with the unit square). Can be |
---|
| 26 | ; scalar or array. (x0,y0), (x1,y1), (x2,y2) and (x3,y3) are |
---|
| 27 | ; given in the anticlockwise order. |
---|
[125] | 28 | ; @param y1in {in}{required} the coordinates of the quadrilateral |
---|
[101] | 29 | ; (see above for correspondance with the unit square). Can be |
---|
| 30 | ; scalar or array. (x0,y0), (x1,y1), (x2,y2) and (x3,y3) are |
---|
| 31 | ; given in the anticlockwise order. |
---|
[125] | 32 | ; @param x2in {in}{required} the coordinates of the quadrilateral |
---|
[101] | 33 | ; (see above for correspondance with the unit square). Can be |
---|
| 34 | ; scalar or array. (x0,y0), (x1,y1), (x2,y2) and (x3,y3) are |
---|
| 35 | ; given in the anticlockwise order. |
---|
[125] | 36 | ; @param y2in {in}{required} the coordinates of the quadrilateral |
---|
[101] | 37 | ; (see above for correspondance with the unit square). Can be |
---|
| 38 | ; scalar or array. (x0,y0), (x1,y1), (x2,y2) and (x3,y3) are |
---|
| 39 | ; given in the anticlockwise order. |
---|
[125] | 40 | ; @param x3in {in}{required} the coordinates of the quadrilateral |
---|
[101] | 41 | ; (see above for correspondance with the unit square). Can be |
---|
| 42 | ; scalar or array. (x0,y0), (x1,y1), (x2,y2) and (x3,y3) are |
---|
| 43 | ; given in the anticlockwise order. |
---|
[125] | 44 | ; @param y3in {in}{required} the coordinates of the quadrilateral |
---|
[101] | 45 | ; (see above for correspondance with the unit square). Can be |
---|
| 46 | ; scalar or array. (x0,y0), (x1,y1), (x2,y2) and (x3,y3) are |
---|
| 47 | ; given in the anticlockwise order. |
---|
[59] | 48 | ; |
---|
[125] | 49 | ; @param xxin {in}{optional} the coordinates of the point(s) for which we want to do the |
---|
[59] | 50 | ; mapping. Can be scalar or array. |
---|
[125] | 51 | ; @param yyin {in}{optional} the coordinates of the point(s) for which we want to do the |
---|
[101] | 52 | ; mapping. Can be scalar or array. |
---|
[59] | 53 | ; |
---|
[101] | 54 | ; @returns |
---|
[59] | 55 | ; |
---|
| 56 | ; (2,n) array: the new coodinates (xout, yout) of the (xin,yin) |
---|
[125] | 57 | ; point(s) after mapping. |
---|
[59] | 58 | ; If xin is a scalar, then n is equal to the number of elements of |
---|
| 59 | ; x0. If xin is an array , then n is equal to the number of |
---|
| 60 | ; elements of xin. |
---|
| 61 | ; If xin and yin are omited, square2quadrilateral returns the |
---|
[125] | 62 | ; matrix A which is used for the inverse transformation. |
---|
[59] | 63 | ; |
---|
| 64 | ; |
---|
[101] | 65 | ; @restrictions I think degenerated quadrilateral (e.g. flat of |
---|
[59] | 66 | ; twisted) is not work. This has to be tested. |
---|
| 67 | ; |
---|
[125] | 68 | ; @examples |
---|
[59] | 69 | ; |
---|
| 70 | ; IDL> splot,[0,5],[0,3],/nodata,xstyle=1,ystyle=1 |
---|
| 71 | ; IDL> tracegrille, findgen(11)*.1, findgen(11)*.1,color=indgen(12)*20 |
---|
| 72 | ; IDL> xin = (findgen(11)*.1)#replicate(1, 11) |
---|
| 73 | ; IDL> yin = replicate(1, 11)#(findgen(11)*.1) |
---|
| 74 | ; IDL> out = square2quadrilateral(2,1,3,0,5,1,2,3, xin, yin) |
---|
| 75 | ; IDL> tracegrille, reform(out[0,*],11,11), reform(out[1,*],11,11),color=indgen(12)*20 |
---|
| 76 | ; |
---|
[101] | 77 | ; @history |
---|
| 78 | ; Sebastien Masson (smasson\@lodyc.jussieu.fr) |
---|
[59] | 79 | ; August 2003 |
---|
| 80 | ; Based on "Digital Image Warping" by G. Wolberg |
---|
| 81 | ; IEEE Computer Society Press, Los Alamitos, California |
---|
| 82 | ; Chapter 3, see p 52-56 |
---|
[118] | 83 | ; |
---|
[125] | 84 | ; |
---|
[118] | 85 | ; @version $Id$ |
---|
| 86 | ; |
---|
[59] | 87 | ;- |
---|
| 88 | ;------------------------------------------------------------ |
---|
| 89 | ;------------------------------------------------------------ |
---|
| 90 | ;------------------------------------------------------------ |
---|
| 91 | FUNCTION square2quadrilateral, x0in, y0in, x1in, y1in, x2in, y2in, x3in, y3in, xxin, yyin |
---|
| 92 | ; |
---|
| 93 | ; Warning, wrong definition of (x2,y2) and (x3,y3) at the bottom of |
---|
| 94 | ; page 54 of Wolberg's book, see figure 3.7 page 56 for the good |
---|
| 95 | ; definition. |
---|
| 96 | ; |
---|
[114] | 97 | ; |
---|
| 98 | compile_opt idl2, strictarrsubs |
---|
| 99 | ; |
---|
[59] | 100 | IF keyword_set(double) THEN BEGIN |
---|
| 101 | x0 = double(x0in) |
---|
| 102 | x1 = double(x1in) |
---|
| 103 | x2 = double(x2in) |
---|
| 104 | x3 = double(x3in) |
---|
| 105 | y0 = double(y0in) |
---|
| 106 | y1 = double(y1in) |
---|
| 107 | y2 = double(y2in) |
---|
| 108 | y3 = double(y3in) |
---|
| 109 | IF arg_present(xxin) THEN BEGIN |
---|
| 110 | xin = double(xxin) |
---|
| 111 | yin = double(yyin) |
---|
| 112 | ENDIF |
---|
| 113 | ENDIF ELSE BEGIN |
---|
| 114 | x0 = float(x0in) |
---|
| 115 | x1 = float(x1in) |
---|
| 116 | x2 = float(x2in) |
---|
| 117 | x3 = float(x3in) |
---|
| 118 | y0 = float(y0in) |
---|
| 119 | y1 = float(y1in) |
---|
| 120 | y2 = float(y2in) |
---|
| 121 | y3 = float(y3in) |
---|
| 122 | IF arg_present(xxin) THEN BEGIN |
---|
| 123 | xin = float(xxin) |
---|
| 124 | yin = float(yyin) |
---|
| 125 | ENDIF |
---|
| 126 | ENDELSE |
---|
| 127 | ; |
---|
[114] | 128 | IF keyword_set(double) THEN a = dblarr(8, n_elements(x0)) $ |
---|
[125] | 129 | ELSE a = fltarr(8, n_elements(x0)) |
---|
[59] | 130 | ; |
---|
| 131 | delx3 = x0-x1+x2-x3 |
---|
| 132 | dely3 = y0-y1+y2-y3 |
---|
| 133 | ; |
---|
| 134 | affinemap = where(delx3 EQ 0 AND dely3 EQ 0) |
---|
| 135 | IF affinemap[0] NE -1 THEN BEGIN |
---|
| 136 | xx0 = x0[affinemap] |
---|
| 137 | xx1 = x1[affinemap] |
---|
| 138 | xx2 = x2[affinemap] |
---|
| 139 | yy0 = y0[affinemap] |
---|
| 140 | yy1 = y1[affinemap] |
---|
| 141 | yy2 = y2[affinemap] |
---|
| 142 | ; |
---|
| 143 | a[0, affinemap] = xx1-xx0 |
---|
| 144 | a[1, affinemap] = xx2-xx1 |
---|
| 145 | a[2, affinemap] = xx0 |
---|
| 146 | a[3, affinemap] = yy1-yy0 |
---|
| 147 | a[4, affinemap] = yy2-yy1 |
---|
| 148 | a[5, affinemap] = yy0 |
---|
| 149 | a[6, affinemap] = 0 |
---|
| 150 | a[7, affinemap] = 0 |
---|
| 151 | ENDIF |
---|
| 152 | ; |
---|
| 153 | projectivemap = where(delx3 NE 0 OR dely3 NE 0) |
---|
| 154 | IF projectivemap[0] NE -1 THEN BEGIN |
---|
| 155 | xx0 = x0[projectivemap] |
---|
| 156 | xx1 = x1[projectivemap] |
---|
| 157 | xx2 = x2[projectivemap] |
---|
| 158 | xx3 = x3[projectivemap] |
---|
| 159 | yy0 = y0[projectivemap] |
---|
| 160 | yy1 = y1[projectivemap] |
---|
| 161 | yy2 = y2[projectivemap] |
---|
| 162 | yy3 = y3[projectivemap] |
---|
[125] | 163 | ; |
---|
[59] | 164 | delx1 = xx1-xx2 |
---|
| 165 | dely1 = yy1-yy2 |
---|
| 166 | delx2 = xx3-xx2 |
---|
| 167 | dely2 = yy3-yy2 |
---|
| 168 | delx3 = delx3[projectivemap] |
---|
| 169 | dely3 = dely3[projectivemap] |
---|
| 170 | ; |
---|
| 171 | div = delx1*dely2-dely1*delx2 |
---|
| 172 | zero = where(div EQ 0) |
---|
| 173 | IF zero[0] NE -1 THEN BEGIN |
---|
| 174 | stop |
---|
| 175 | ENDIF |
---|
| 176 | a13 = (delx3*dely2-dely3*delx2)/div |
---|
| 177 | a23 = (delx1*dely3-dely1*delx3)/div |
---|
| 178 | ; |
---|
| 179 | a[0, projectivemap] = xx1-xx0+a13*xx1 |
---|
| 180 | a[1, projectivemap] = xx3-xx0+a23*xx3 |
---|
| 181 | a[2, projectivemap] = xx0 |
---|
| 182 | a[3, projectivemap] = yy1-yy0+a13*yy1 |
---|
| 183 | a[4, projectivemap] = yy3-yy0+a23*yy3 |
---|
| 184 | a[5, projectivemap] = yy0 |
---|
| 185 | a[6, projectivemap] = a13 |
---|
| 186 | a[7, projectivemap] = a23 |
---|
| 187 | ENDIF |
---|
[125] | 188 | ; |
---|
[59] | 189 | IF NOT arg_present(xxin) THEN return, a |
---|
| 190 | ; |
---|
| 191 | IF n_elements(xin) EQ 1 THEN BEGIN |
---|
[125] | 192 | xin = replicate(xin, n_elements(x0)) |
---|
| 193 | yin = replicate(yin, n_elements(x0)) |
---|
[59] | 194 | ENDIF |
---|
| 195 | ; |
---|
| 196 | IF keyword_set(double) THEN res = dblarr(2, n_elements(xin)) $ |
---|
| 197 | ELSE res = fltarr(2, n_elements(xin)) |
---|
| 198 | IF n_elements(x0) EQ 1 THEN BEGIN |
---|
| 199 | div = a[6]*xin[*] + a[7]*yin[*] + 1 |
---|
| 200 | zero = where(div EQ 0) |
---|
| 201 | IF zero[0] NE -1 THEN BEGIN |
---|
| 202 | stop |
---|
| 203 | ENDIF |
---|
| 204 | res[0, *] = (a[0]*xin[*] + a[1]*yin[*] + a[2])/div |
---|
| 205 | res[1, *] = (a[3]*xin[*] + a[4]*yin[*] + a[5])/div |
---|
| 206 | ENDIF ELSE BEGIN |
---|
| 207 | div = a[6, *]*xin +a[7, *]*yin + 1 |
---|
| 208 | zero = where(div EQ 0) |
---|
| 209 | IF zero[0] NE -1 THEN BEGIN |
---|
| 210 | stop |
---|
| 211 | ENDIF |
---|
| 212 | res[0, *] = (a[0, *]*xin[*] + a[1, *]*yin[*] + a[2, *])/div |
---|
| 213 | res[1, *] = (a[3, *]*xin[*] + a[4, *]*yin[*] + a[5, *])/div |
---|
| 214 | ENDELSE |
---|
| 215 | ; |
---|
| 216 | RETURN, res |
---|
| 217 | END |
---|