source: trunk/SRC/Interpolation/square2quadrilateral.pro

;+
;
; @file_comments
; warm (or map) a unit square onto an arbitrary quadrilateral
; according to the 4-point correspondences:
;  - (0,0) -> (x0,y0)
;  - (1,0) -> (x1,y1)
;  - (1,1) -> (x2,y2)
;  - (0,1) -> (x3,y3)
;
; The mapping is done using perspective transformation which preserve
; lines in all orientations and permit quadrilateral to quadrilateral
; mappings. see ref. below.
;
; @categories
; Picture, Grid
;
; @param x0in {in}{required}
; @param y0in {in}{required}
; @param x1in {in}{required}
; @param y1in {in}{required}
; @param x2in {in}{required}
; @param y2in {in}{required}
; @param x3in {in}{required}
; @param y3in {in}{required}
; the coordinates of the quadrilateral (see above for correspondence with the
; unit square).
; Can be scalar or array.
; (x0,y0), (x1,y1), (x2,y2) and (x3,y3) are given in the anticlockwise order.
;
; @param xxin {in}{optional}
; first coordinates of the point(s) for which we want to do the mapping.
; @param yyin {in}{optional}
; second coordinates of the point(s) for which we want to do the mapping.
;
; @keyword DOUBLE {type=salar 0 or 1}{default=0}
; activate to perform double precision computation
;
; @returns
; (2,n) array: the new coordinates (xout,yout) of the (xin,yin)
; point(s) after mapping.
; If xin is a scalar, then n is equal to the number of elements of
; x0. If xin is an array , then n is equal to the number of
; elements of xin.
; If xin and yin are omitted, <pro>square2quadrilateral</pro> returns the
; matrix A which is used for the inverse transformation.
;
; @restrictions
; I think degenerated quadrilateral (e.g. flat of twisted) is not work.
; This has to be tested.
;
; @examples
;
;   IDL> splot,[0,5],[0,3],/nodata,xstyle=1,ystyle=1
;   IDL> tracegrille, findgen(11)*.1, findgen(11)*.1,color=indgen(12)*20
;   IDL> xin = (findgen(11)*.1)#replicate(1, 11)
;   IDL> yin = replicate(1, 11)#(findgen(11)*.1)
;   IDL> out = square2quadrilateral(2,1,3,0,5,1,2,3, xin, yin)
;   IDL> tracegrille, reform(out[0,*],11,11), reform(out[1,*],11,11),color=indgen(12)*20
;
; @history
; Sebastien Masson (smasson\@lodyc.jussieu.fr)
;  - August 2003
;    Based on "Digital Image Warping" by G. Wolberg
;    IEEE Computer Society Press, Los Alamitos, California
;    Chapter 3, see p 52-56
;
; @version
; $Id$
;
;-
FUNCTION square2quadrilateral, x0in, y0in, x1in, y1in, x2in, y2in $
                             , x3in, y3in, xxin, yyin $
                             , DOUBLE=double
;
; Warning, wrong definition of (x2,y2) and (x3,y3) at the bottom of
; page 54 of Wolberg's book, see figure 3.7 page 56 for the good
; definition.
;
  compile_opt idl2, strictarrsubs
;
  IF keyword_set(double) THEN BEGIN
    x0 = double(x0in)
    x1 = double(x1in)
    x2 = double(x2in)
    x3 = double(x3in)
    y0 = double(y0in)
    y1 = double(y1in)
    y2 = double(y2in)
    y3 = double(y3in)
    IF arg_present(xxin) THEN BEGIN
      xin = double(xxin)
      yin = double(yyin)
    ENDIF
  ENDIF ELSE BEGIN
    x0 = float(x0in)
    x1 = float(x1in)
    x2 = float(x2in)
    x3 = float(x3in)
    y0 = float(y0in)
    y1 = float(y1in)
    y2 = float(y2in)
    y3 = float(y3in)
    IF arg_present(xxin) THEN BEGIN
      xin = float(xxin)
      yin = float(yyin)
    ENDIF
  ENDELSE
;
  IF keyword_set(double) THEN a = dblarr(8, n_elements(x0)) $
  ELSE a = fltarr(8, n_elements(x0))
;
  delx3 = x0-x1+x2-x3
  dely3 = y0-y1+y2-y3
;
  affinemap = where(delx3 EQ 0 AND dely3 EQ 0)
  IF affinemap[0] NE -1 THEN BEGIN
    xx0 = x0[affinemap]
    xx1 = x1[affinemap]
    xx2 = x2[affinemap]
    yy0 = y0[affinemap]
    yy1 = y1[affinemap]
    yy2 = y2[affinemap]
;
    a[0, affinemap] = xx1-xx0
    a[1, affinemap] = xx2-xx1
    a[2, affinemap] = xx0
    a[3, affinemap] = yy1-yy0
    a[4, affinemap] = yy2-yy1
    a[5, affinemap] = yy0
    a[6, affinemap] = 0
    a[7, affinemap] = 0
  ENDIF
;
  projectivemap = where(delx3 NE 0 OR dely3 NE 0)
  IF projectivemap[0] NE -1 THEN BEGIN
    xx0 = x0[projectivemap]
    xx1 = x1[projectivemap]
    xx2 = x2[projectivemap]
    xx3 = x3[projectivemap]
    yy0 = y0[projectivemap]
    yy1 = y1[projectivemap]
    yy2 = y2[projectivemap]
    yy3 = y3[projectivemap]
;
    delx1 = xx1-xx2
    dely1 = yy1-yy2
    delx2 = xx3-xx2
    dely2 = yy3-yy2
    delx3 = delx3[projectivemap]
    dely3 = dely3[projectivemap]
;
    div = delx1*dely2-dely1*delx2
    zero = where(div EQ 0)
    IF zero[0] NE -1 THEN BEGIN
      stop
    ENDIF
    a13 = (delx3*dely2-dely3*delx2)/div
    a23 = (delx1*dely3-dely1*delx3)/div
;
    a[0, projectivemap] = xx1-xx0+a13*xx1
    a[1, projectivemap] = xx3-xx0+a23*xx3
    a[2, projectivemap] = xx0
    a[3, projectivemap] = yy1-yy0+a13*yy1
    a[4, projectivemap] = yy3-yy0+a23*yy3
    a[5, projectivemap] = yy0
    a[6, projectivemap] = a13
    a[7, projectivemap] = a23
  ENDIF
;
  IF NOT arg_present(xxin) THEN return, a
;
  IF n_elements(xin) EQ 1 THEN BEGIN
    xin = replicate(xin, n_elements(x0))
    yin = replicate(yin, n_elements(x0))
  ENDIF
;
  IF keyword_set(double) THEN res = dblarr(2, n_elements(xin)) $
  ELSE res = fltarr(2, n_elements(xin))
  IF n_elements(x0) EQ 1 THEN BEGIN
    div = a[6]*xin[*] + a[7]*yin[*] + 1
    zero = where(div EQ 0)
    IF zero[0] NE -1 THEN BEGIN
      stop
    ENDIF
    res[0, *] = (a[0]*xin[*] + a[1]*yin[*] + a[2])/div
    res[1, *] = (a[3]*xin[*] + a[4]*yin[*] + a[5])/div
  ENDIF ELSE BEGIN
    div = a[6, *]*xin +a[7, *]*yin + 1
    zero = where(div EQ 0)
    IF zero[0] NE -1 THEN BEGIN
      stop
    ENDIF
    res[0, *] = (a[0, *]*xin[*] + a[1, *]*yin[*] + a[2, *])/div
    res[1, *] = (a[3, *]*xin[*] + a[4, *]*yin[*] + a[5, *])/div
  ENDELSE
;
  RETURN, res
END