source: trunk/SRC/Interpolation/square2quadrilateral.pro @ 228

Last change on this file since 228 was 163, checked in by navarro, 18 years ago

header improvements : type of parameters and keywords, default values, spell checking + idldoc assistant (IDL online_help)

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