1 | ;+ |
---|
2 | ; |
---|
3 | ; @file_comments |
---|
4 | ; warm (or map) a unit square onto an arbitrary quadrilateral |
---|
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 | ; |
---|
14 | ; @categories |
---|
15 | ; Picture, Grid |
---|
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 (see above for correspondence with the |
---|
26 | ; unit square). |
---|
27 | ; Can be scalar or array. |
---|
28 | ; (x0,y0), (x1,y1), (x2,y2) and (x3,y3) are given in the anticlockwise order. |
---|
29 | ; |
---|
30 | ; @param xxin {in}{optional} |
---|
31 | ; first coordinates of the point(s) for which we want to do the mapping. |
---|
32 | ; @param yyin {in}{optional} |
---|
33 | ; second coordinates of the point(s) for which we want to do the mapping. |
---|
34 | ; |
---|
35 | ; @keyword DOUBLE {type=salar 0 or 1}{default=0} |
---|
36 | ; activate to perform double precision computation |
---|
37 | ; |
---|
38 | ; @returns |
---|
39 | ; (2,n) array: the new coordinates (xout,yout) of the (xin,yin) |
---|
40 | ; point(s) after mapping. |
---|
41 | ; If xin is a scalar, then n is equal to the number of elements of |
---|
42 | ; x0. If xin is an array , then n is equal to the number of |
---|
43 | ; elements of xin. |
---|
44 | ; If xin and yin are omitted, <pro>square2quadrilateral</pro> returns the |
---|
45 | ; matrix A which is used for the inverse transformation. |
---|
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 | ; @history |
---|
61 | ; Sebastien Masson (smasson\@lodyc.jussieu.fr) |
---|
62 | ; August 2003 |
---|
63 | ; Based on "Digital Image Warping" by G. Wolberg |
---|
64 | ; IEEE Computer Society Press, Los Alamitos, California |
---|
65 | ; Chapter 3, see p 52-56 |
---|
66 | ; |
---|
67 | ; |
---|
68 | ; @version |
---|
69 | ; $Id$ |
---|
70 | ; |
---|
71 | ;- |
---|
72 | ; |
---|
73 | FUNCTION square2quadrilateral, x0in, y0in, x1in, y1in, x2in, y2in, x3in, y3in, xxin, yyin, DOUBLE = double |
---|
74 | ; |
---|
75 | ; Warning, wrong definition of (x2,y2) and (x3,y3) at the bottom of |
---|
76 | ; page 54 of Wolberg's book, see figure 3.7 page 56 for the good |
---|
77 | ; definition. |
---|
78 | ; |
---|
79 | compile_opt idl2, strictarrsubs |
---|
80 | ; |
---|
81 | IF keyword_set(double) THEN BEGIN |
---|
82 | x0 = double(x0in) |
---|
83 | x1 = double(x1in) |
---|
84 | x2 = double(x2in) |
---|
85 | x3 = double(x3in) |
---|
86 | y0 = double(y0in) |
---|
87 | y1 = double(y1in) |
---|
88 | y2 = double(y2in) |
---|
89 | y3 = double(y3in) |
---|
90 | IF arg_present(xxin) THEN BEGIN |
---|
91 | xin = double(xxin) |
---|
92 | yin = double(yyin) |
---|
93 | ENDIF |
---|
94 | ENDIF ELSE BEGIN |
---|
95 | x0 = float(x0in) |
---|
96 | x1 = float(x1in) |
---|
97 | x2 = float(x2in) |
---|
98 | x3 = float(x3in) |
---|
99 | y0 = float(y0in) |
---|
100 | y1 = float(y1in) |
---|
101 | y2 = float(y2in) |
---|
102 | y3 = float(y3in) |
---|
103 | IF arg_present(xxin) THEN BEGIN |
---|
104 | xin = float(xxin) |
---|
105 | yin = float(yyin) |
---|
106 | ENDIF |
---|
107 | ENDELSE |
---|
108 | ; |
---|
109 | IF keyword_set(double) THEN a = dblarr(8, n_elements(x0)) $ |
---|
110 | ELSE a = fltarr(8, n_elements(x0)) |
---|
111 | ; |
---|
112 | delx3 = x0-x1+x2-x3 |
---|
113 | dely3 = y0-y1+y2-y3 |
---|
114 | ; |
---|
115 | affinemap = where(delx3 EQ 0 AND dely3 EQ 0) |
---|
116 | IF affinemap[0] NE -1 THEN BEGIN |
---|
117 | xx0 = x0[affinemap] |
---|
118 | xx1 = x1[affinemap] |
---|
119 | xx2 = x2[affinemap] |
---|
120 | yy0 = y0[affinemap] |
---|
121 | yy1 = y1[affinemap] |
---|
122 | yy2 = y2[affinemap] |
---|
123 | ; |
---|
124 | a[0, affinemap] = xx1-xx0 |
---|
125 | a[1, affinemap] = xx2-xx1 |
---|
126 | a[2, affinemap] = xx0 |
---|
127 | a[3, affinemap] = yy1-yy0 |
---|
128 | a[4, affinemap] = yy2-yy1 |
---|
129 | a[5, affinemap] = yy0 |
---|
130 | a[6, affinemap] = 0 |
---|
131 | a[7, affinemap] = 0 |
---|
132 | ENDIF |
---|
133 | ; |
---|
134 | projectivemap = where(delx3 NE 0 OR dely3 NE 0) |
---|
135 | IF projectivemap[0] NE -1 THEN BEGIN |
---|
136 | xx0 = x0[projectivemap] |
---|
137 | xx1 = x1[projectivemap] |
---|
138 | xx2 = x2[projectivemap] |
---|
139 | xx3 = x3[projectivemap] |
---|
140 | yy0 = y0[projectivemap] |
---|
141 | yy1 = y1[projectivemap] |
---|
142 | yy2 = y2[projectivemap] |
---|
143 | yy3 = y3[projectivemap] |
---|
144 | ; |
---|
145 | delx1 = xx1-xx2 |
---|
146 | dely1 = yy1-yy2 |
---|
147 | delx2 = xx3-xx2 |
---|
148 | dely2 = yy3-yy2 |
---|
149 | delx3 = delx3[projectivemap] |
---|
150 | dely3 = dely3[projectivemap] |
---|
151 | ; |
---|
152 | div = delx1*dely2-dely1*delx2 |
---|
153 | zero = where(div EQ 0) |
---|
154 | IF zero[0] NE -1 THEN BEGIN |
---|
155 | stop |
---|
156 | ENDIF |
---|
157 | a13 = (delx3*dely2-dely3*delx2)/div |
---|
158 | a23 = (delx1*dely3-dely1*delx3)/div |
---|
159 | ; |
---|
160 | a[0, projectivemap] = xx1-xx0+a13*xx1 |
---|
161 | a[1, projectivemap] = xx3-xx0+a23*xx3 |
---|
162 | a[2, projectivemap] = xx0 |
---|
163 | a[3, projectivemap] = yy1-yy0+a13*yy1 |
---|
164 | a[4, projectivemap] = yy3-yy0+a23*yy3 |
---|
165 | a[5, projectivemap] = yy0 |
---|
166 | a[6, projectivemap] = a13 |
---|
167 | a[7, projectivemap] = a23 |
---|
168 | ENDIF |
---|
169 | ; |
---|
170 | IF NOT arg_present(xxin) THEN return, a |
---|
171 | ; |
---|
172 | IF n_elements(xin) EQ 1 THEN BEGIN |
---|
173 | xin = replicate(xin, n_elements(x0)) |
---|
174 | yin = replicate(yin, n_elements(x0)) |
---|
175 | ENDIF |
---|
176 | ; |
---|
177 | IF keyword_set(double) THEN res = dblarr(2, n_elements(xin)) $ |
---|
178 | ELSE res = fltarr(2, n_elements(xin)) |
---|
179 | IF n_elements(x0) EQ 1 THEN BEGIN |
---|
180 | div = a[6]*xin[*] + a[7]*yin[*] + 1 |
---|
181 | zero = where(div EQ 0) |
---|
182 | IF zero[0] NE -1 THEN BEGIN |
---|
183 | stop |
---|
184 | ENDIF |
---|
185 | res[0, *] = (a[0]*xin[*] + a[1]*yin[*] + a[2])/div |
---|
186 | res[1, *] = (a[3]*xin[*] + a[4]*yin[*] + a[5])/div |
---|
187 | ENDIF ELSE BEGIN |
---|
188 | div = a[6, *]*xin +a[7, *]*yin + 1 |
---|
189 | zero = where(div EQ 0) |
---|
190 | IF zero[0] NE -1 THEN BEGIN |
---|
191 | stop |
---|
192 | ENDIF |
---|
193 | res[0, *] = (a[0, *]*xin[*] + a[1, *]*yin[*] + a[2, *])/div |
---|
194 | res[1, *] = (a[3, *]*xin[*] + a[4, *]*yin[*] + a[5, *])/div |
---|
195 | ENDELSE |
---|
196 | ; |
---|
197 | RETURN, res |
---|
198 | END |
---|