[59] | 1 | ;+ |
---|
| 2 | ; |
---|
[125] | 3 | ; @file_comments |
---|
| 4 | ; warm (or map) an arbitrary quadrilateral onto a unit square |
---|
[59] | 5 | ; according to the 4-point correspondences: |
---|
| 6 | ; (x0,y0) -> (0,0) |
---|
| 7 | ; (x1,y1) -> (1,0) |
---|
| 8 | ; (x2,y2) -> (1,1) |
---|
| 9 | ; (x3,y3) -> (0,1) |
---|
[231] | 10 | ; This is the inverse function of <pro>square2quadrilateral</pro>. |
---|
| 11 | ; |
---|
[59] | 12 | ; The mapping is done using perspective transformation which preserve |
---|
| 13 | ; lines in all orientations and permit quadrilateral to quadrilateral |
---|
[358] | 14 | ; mappings. see ref. below. |
---|
[59] | 15 | ; |
---|
[231] | 16 | ; @categories |
---|
[157] | 17 | ; Picture, Grid |
---|
[59] | 18 | ; |
---|
[118] | 19 | ; @param x0in {in}{required} |
---|
| 20 | ; @param y0in {in}{required} |
---|
| 21 | ; @param x1in {in}{required} |
---|
| 22 | ; @param y1in {in}{required} |
---|
| 23 | ; @param x2in {in}{required} |
---|
| 24 | ; @param y2in {in}{required} |
---|
| 25 | ; @param x3in {in}{required} |
---|
| 26 | ; @param y3in {in}{required} |
---|
| 27 | ; the coordinates of the quadrilateral |
---|
[163] | 28 | ; (see above for correspondence with the unit square). Can be |
---|
[118] | 29 | ; scalar or array. (x0,y0), (x1,y1), (x2,y2) and (x3,y3) are |
---|
| 30 | ; given in the anticlockwise order. |
---|
[59] | 31 | ; |
---|
[136] | 32 | ; @param xxin {in}{required} |
---|
| 33 | ; the coordinates of the point(s) for which we want to do the mapping. |
---|
| 34 | ; Can be scalar or array. |
---|
[59] | 35 | ; |
---|
[136] | 36 | ; @param yyin {in}{required} |
---|
| 37 | ; the coordinates of the point(s) for which we want to do the mapping. |
---|
| 38 | ; Can be scalar or array. |
---|
| 39 | ; |
---|
[282] | 40 | ; @keyword PERF {type=salar 0 or 1}{default=0} |
---|
| 41 | ; activate to print the elapsed time spent within quadrilateral2square |
---|
[118] | 42 | ; |
---|
[282] | 43 | ; @keyword DOUBLE {type=salar 0 or 1}{default=0} |
---|
| 44 | ; activate to perform double precision computation |
---|
| 45 | ; |
---|
[101] | 46 | ; @returns |
---|
[242] | 47 | ; (2,n) array: the new coordinates (xout,yout) of the (xin,yin) point(s) after |
---|
[136] | 48 | ; mapping. |
---|
| 49 | ; If xin is a scalar, then n is equal to the number of elements of x0. |
---|
[242] | 50 | ; If xin is an array, then n is equal to the number of elements of xin. |
---|
[59] | 51 | ; |
---|
[125] | 52 | ; @restrictions |
---|
| 53 | ; I think degenerated quadrilateral (e.g. flat of twisted) is not work. |
---|
| 54 | ; This has to be tested. |
---|
[59] | 55 | ; |
---|
[125] | 56 | ; @examples |
---|
[59] | 57 | ; |
---|
| 58 | ; IDL> splot,[0,5],[0,3],/nodata,xstyle=1,ystyle=1 |
---|
| 59 | ; IDL> tracegrille, findgen(11)*.1, findgen(11)*.1,color=indgen(12)*20 |
---|
| 60 | ; IDL> xin = (findgen(11)*.1)#replicate(1, 11) |
---|
| 61 | ; IDL> yin = replicate(1, 11)#(findgen(11)*.1) |
---|
| 62 | ; IDL> out = square2quadrilateral(2,1,3,0,5,1,2,3, xin, yin) |
---|
| 63 | ; IDL> tracegrille, reform(out[0,*],11,11), reform(out[1,*],11,11),color=indgen(12)*20 |
---|
| 64 | ; |
---|
| 65 | ; IDL> inorg=quadrilateral2square(2,1,3,0,5,1,2,3,out[0,*],out[1,*]) |
---|
| 66 | ; IDL> tracegrille, reform(inorg[0,*],11,11), reform(inorg[1,*],11,11),color=indgen(12)*20 |
---|
| 67 | ; |
---|
[101] | 68 | ; @history |
---|
| 69 | ; Sebastien Masson (smasson\@lodyc.jussieu.fr) |
---|
[59] | 70 | ; August 2003 |
---|
| 71 | ; Based on "Digital Image Warping" by G. Wolberg |
---|
| 72 | ; IEEE Computer Society Press, Los Alamitos, California |
---|
| 73 | ; Chapter 3, see p 52-56 |
---|
[118] | 74 | ; |
---|
[125] | 75 | ; |
---|
[231] | 76 | ; @version |
---|
| 77 | ; $Id$ |
---|
[118] | 78 | ; |
---|
[59] | 79 | ;- |
---|
[327] | 80 | FUNCTION quadrilateral2square, x0in, y0in, x1in, y1in, x2in, y2in $ |
---|
| 81 | , x3in, y3in, xxin, yyin $ |
---|
| 82 | , PERF=perf, DOUBLE=double |
---|
[59] | 83 | ; |
---|
[114] | 84 | compile_opt idl2, strictarrsubs |
---|
| 85 | ; |
---|
[59] | 86 | tempsone = systime(1) |
---|
| 87 | ; |
---|
| 88 | ; Warning, wrong definition of (x2,y2) and (x3,y3) at the bottom of |
---|
| 89 | ; page 54 of Wolberg's book, see figure 3.7 page 56 for the good |
---|
| 90 | ; definition. |
---|
| 91 | ; |
---|
| 92 | IF keyword_set(double) THEN BEGIN |
---|
| 93 | x0 = double(x0in) |
---|
| 94 | x1 = double(x1in) |
---|
| 95 | x2 = double(x2in) |
---|
| 96 | x3 = double(x3in) |
---|
| 97 | y0 = double(y0in) |
---|
| 98 | y1 = double(y1in) |
---|
| 99 | y2 = double(y2in) |
---|
| 100 | y3 = double(y3in) |
---|
| 101 | xin = double(xxin) |
---|
| 102 | yin = double(yyin) |
---|
| 103 | ENDIF ELSE BEGIN |
---|
| 104 | x0 = float(x0in) |
---|
| 105 | x1 = float(x1in) |
---|
| 106 | x2 = float(x2in) |
---|
| 107 | x3 = float(x3in) |
---|
| 108 | y0 = float(y0in) |
---|
| 109 | y1 = float(y1in) |
---|
| 110 | y2 = float(y2in) |
---|
| 111 | y3 = float(y3in) |
---|
| 112 | xin = float(xxin) |
---|
| 113 | yin = float(yyin) |
---|
| 114 | ENDELSE |
---|
| 115 | ; |
---|
| 116 | ; get the matrix A |
---|
| 117 | ; |
---|
[282] | 118 | a = square2quadrilateral(x0in, y0in, x1in, y1in, x2in, y2in, x3in, y3in, DOUBLE = double) |
---|
[59] | 119 | ; |
---|
| 120 | ; compute the adjoint matrix |
---|
[125] | 121 | ; |
---|
[59] | 122 | IF keyword_set(double) THEN adj = dblarr(9, n_elements(x0)) $ |
---|
| 123 | ELSE adj = fltarr(9, n_elements(x0)) |
---|
| 124 | ; |
---|
| 125 | adj[0, *] = a[4, *] -a[7, *]*a[5, *] |
---|
| 126 | adj[1, *] = a[7, *]*a[2, *]-a[1, *] |
---|
| 127 | adj[2, *] = a[1, *]*a[5, *]-a[4, *]*a[2, *] |
---|
| 128 | adj[3, *] = a[6, *]*a[5, *]-a[3, *] |
---|
| 129 | adj[4, *] = a[0, *] -a[6, *]*a[2, *] |
---|
| 130 | adj[5, *] = a[3, *]*a[2, *]-a[0, *]*a[5, *] |
---|
| 131 | adj[6, *] = a[3, *]*a[7, *]-a[6, *]*a[4, *] |
---|
| 132 | adj[7, *] = a[6, *]*a[1, *]-a[0, *]*a[7, *] |
---|
| 133 | adj[8, *] = a[0, *]*a[4, *]-a[3, *]*a[1, *] |
---|
[125] | 134 | ; |
---|
[59] | 135 | IF n_elements(xin) EQ 1 THEN BEGIN |
---|
[125] | 136 | xin = replicate(xin, n_elements(x0)) |
---|
| 137 | yin = replicate(yin, n_elements(x0)) |
---|
[59] | 138 | ENDIF |
---|
| 139 | ; |
---|
| 140 | ; compute xprime, yprime and wprime |
---|
| 141 | ; |
---|
| 142 | IF n_elements(x0) EQ 1 THEN BEGIN |
---|
| 143 | wpr = 1./(adj[6]*xin + adj[7]*yin + adj[8]) |
---|
| 144 | ENDIF ELSE BEGIN |
---|
| 145 | wpr = 1./(adj[6, *]*xin + adj[7, *]*yin + adj[8, *]) |
---|
| 146 | ENDELSE |
---|
| 147 | xpr = xin*wpr |
---|
| 148 | ypr = yin*wpr |
---|
| 149 | ; |
---|
| 150 | IF keyword_set(double) THEN res = dblarr(2, n_elements(xin)) $ |
---|
| 151 | ELSE res = fltarr(2, n_elements(xin)) |
---|
| 152 | ; |
---|
| 153 | IF n_elements(x0) EQ 1 THEN BEGIN |
---|
| 154 | res[0, *] = xpr*adj[0] + ypr*adj[1] +wpr*adj[2] |
---|
| 155 | res[1, *] = xpr*adj[3] + ypr*adj[4] +wpr*adj[5] |
---|
| 156 | ENDIF ELSE BEGIN |
---|
| 157 | res[0, *] = xpr*adj[0, *] + ypr*adj[1, *] +wpr*adj[2, *] |
---|
| 158 | res[1, *] = xpr*adj[3, *] + ypr*adj[4, *] +wpr*adj[5, *] |
---|
| 159 | ENDELSE |
---|
| 160 | ; |
---|
| 161 | IF keyword_set(perf) THEN print, 'time quadrilateral2square', systime(1)-tempsone |
---|
| 162 | |
---|
| 163 | RETURN, res |
---|
| 164 | END |
---|