;+
;
; @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. bellow.
;
; @categories image, grid manipulation
;
; @examples 
;
;     res = square2quadrilateral(x0,y0,x1,y1,x2,y2,x3,y3[,xin,yin])
; 
FUNCTION square2quadrilateral, x0in, y0in, x1in, y1in, x2in, y2in, x3in, y3in, xxin, yyin
;     @param x0in {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 y0in {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 x1in {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 y1in {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 x2in {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 y2in {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 x3in {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 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.
;     If xin and yin are omited, square2quadrilateral 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
;      
;-
;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
FUNCTION square2quadrilateral, x0in, y0in, x1in, y1in, x2in, y2in, x3in, y3in, xxin, yyin
;
; 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)
    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 = dlbarr(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