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