source: trunk/SRC/Interpolation/quadrilateral2square.pro @ 325

;+
;
; @file_comments
; warm (or map) an arbitrary quadrilateral onto a unit square
; according to the 4-point correspondences:
;       (x0,y0) -> (0,0)
;       (x1,y1) -> (1,0)
;       (x2,y2) -> (1,1)
;       (x3,y3) -> (0,1)
; This is the inverse function of <pro>square2quadrilateral</pro>.
;
; The mapping is done using perspective transformation which preserve
; lines in all orientations and permit quadrilateral to quadrilateral
; mappings. see ref. bellow.
;
; @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}{required}
; the coordinates of the point(s) for which we want to do the mapping.
; Can be scalar or array.
;
; @param yyin {in}{required}
; the coordinates of the point(s) for which we want to do the mapping.
; Can be scalar or array.
;
; @keyword PERF {type=salar 0 or 1}{default=0}
; activate to print the elapsed time spent within quadrilateral2square
;
; @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.
;
; @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
;
; IDL> inorg=quadrilateral2square(2,1,3,0,5,1,2,3,out[0,*],out[1,*])
; IDL> tracegrille, reform(inorg[0,*],11,11), reform(inorg[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 quadrilateral2square, x0in, y0in, x1in, y1in, x2in, y2in, x3in, y3in, xxin, yyin, PERF = perf, DOUBLE = double
;
  compile_opt idl2, strictarrsubs
;
tempsone = systime(1)
;
; 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.
;
  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)
    xin = double(xxin)
    yin = double(yyin)
  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)
    xin = float(xxin)
    yin = float(yyin)
  ENDELSE
;
; get the matrix A
;
  a = square2quadrilateral(x0in, y0in, x1in, y1in, x2in, y2in, x3in, y3in, DOUBLE = double)
;
; compute the adjoint matrix
;
  IF keyword_set(double) THEN adj = dblarr(9, n_elements(x0)) $
  ELSE adj = fltarr(9, n_elements(x0))
;
  adj[0, *] = a[4, *]        -a[7, *]*a[5, *]
  adj[1, *] = a[7, *]*a[2, *]-a[1, *]
  adj[2, *] = a[1, *]*a[5, *]-a[4, *]*a[2, *]
  adj[3, *] = a[6, *]*a[5, *]-a[3, *]
  adj[4, *] = a[0, *]        -a[6, *]*a[2, *]
  adj[5, *] = a[3, *]*a[2, *]-a[0, *]*a[5, *]
  adj[6, *] = a[3, *]*a[7, *]-a[6, *]*a[4, *]
  adj[7, *] = a[6, *]*a[1, *]-a[0, *]*a[7, *]
  adj[8, *] = a[0, *]*a[4, *]-a[3, *]*a[1, *]
;
  IF n_elements(xin) EQ 1 THEN BEGIN
    xin = replicate(xin, n_elements(x0))
    yin = replicate(yin, n_elements(x0))
  ENDIF
;
; compute xprime, yprime and wprime
;
  IF n_elements(x0) EQ 1 THEN BEGIN
    wpr = 1./(adj[6]*xin + adj[7]*yin + adj[8])
  ENDIF ELSE BEGIN
    wpr = 1./(adj[6, *]*xin + adj[7, *]*yin + adj[8, *])
  ENDELSE
  xpr = xin*wpr
  ypr = yin*wpr
;
  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
    res[0, *] = xpr*adj[0] + ypr*adj[1] +wpr*adj[2]
    res[1, *] = xpr*adj[3] + ypr*adj[4] +wpr*adj[5]
  ENDIF ELSE BEGIN
    res[0, *] = xpr*adj[0, *] + ypr*adj[1, *] +wpr*adj[2, *]
    res[1, *] = xpr*adj[3, *] + ypr*adj[4, *] +wpr*adj[5, *]
  ENDELSE
;
  IF keyword_set(perf) THEN print, 'time quadrilateral2square', systime(1)-tempsone

  RETURN, res
END