;+
;
; @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 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 image, grid manipulation
;
; @examples 
;
;     res = square2quadrilateral(x0,y0,x1,y1,x2,y2,x3,y3,xin,yin)
; 
;     @param x0in {in}{required}  the coordinates of the quadrilateral
;     @param y0in {in}{required}  the coordinates of the quadrilateral
;     @param x1in {in}{required}  the coordinates of the quadrilateral
;     @param y1in {in}{required}  the coordinates of the quadrilateral
;     @param x2in {in}{required}  the coordinates of the quadrilateral
;     @param y2in {in}{required}  the coordinates of the quadrilateral
;     @param x3in {in}{required}  the coordinates of the quadrilateral
;     @param y3in  {in}{required}  the coordinates of the quadrilateral
;     (see above for correspondance 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.
;
; @returns
;
;     (2,n) array: the new coodinates (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
;      
;-
;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
FUNCTION quadrilateral2square, x0in, y0in, x1in, y1in, x2in, y2in, x3in, y3in, xxin, yyin, PERF = perf
;
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)
;
; 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