source: trunk/SRC/Interpolation/inquad.pro @ 327

;+
;
; @file_comments
; to find if an (x,y) point is in a quadrilateral (x1,x2,x3,x4)
;
; @categories
; Grid
;
; @param x {in}{required}
; @param y {in}{required}
; the coordinates of the point we want to know where it is.
; Must be a scalar if /ONSPHERE activated else can be scalar or array.
;
; @param x1 {in}{required}
; @param y1 {in}{required}
; @param x2 {in}{required}
; @param y2 {in}{required}
; @param x3 {in}{required}
; @param y3 {in}{required}
; @param x4 {in}{required}
; @param y4 {in}{required}
; the coordinates of the quadrilateral given in the CLOCKWISE order.
; Scalar or array.
;
; @keyword DOUBLE
; use double precision to perform the computation
;
; @keyword ONSPHERE
; to specify that the quadrilateral are on a sphere and
; that their coordinates are longitude-latitude coordinates. In this
; case, east-west periodicity, poles singularity and other pbs
; related to longitude-latitude coordinates are managed
; automatically.
;
; @keyword DELTA {default=4}
; to speed up the program, we reduce the aera where we look for potential
; quadrilaterals containing (x,y). Delta defines the limit of the box
; centred on (x,y) with a zonal and meridional extent of delta degrees.
;
; @keyword NOPRINT
; to suppress the print messages.
;
; @keyword NEWCOORD
;
; @returns
; a n elements vector where n is the number of elements of
; x. res[i]=j means that the point number i is located in the
; quadrilateral number j with (0 <= j <= n_elements(x0)-1)
;
; @restrictions
; I think degenerated quadrilateral (e.g. flat of twisted) is not work.
; This has to be tested.
;
; @examples
;
; IDL> x = 1.*[1, 2, 6, 7, 3]
; IDL> y = 1.*[1, 3, 3, 4, 7]
; IDL> x1 = 1.*[0,4,2]
; IDL> y1 = 1.*[1,4,8]
; IDL> x2 = 1.*[1,6,4]
; IDL> y2 = 1.*[5,6,8]
; IDL> x3 = 1.*[3,8,4]
; IDL> y3 = 1.*[4,4,6]
; IDL> x4 = 1.*[2,6,2]
; IDL> y4 = 1.*[0,2,6]
; IDL> splot, [0,10], [0,10], xstyle = 1, ystyle = 1,/nodata
; IDL> for i=0,2 do oplot, [x4[i],x1[i],x2[i],x3[i],x4[i]],[y4[i],y1[i],y2[i],y3[i],y4[i]]
; IDL> oplot, x, y, color = 20, psym = 1, thick = 2
; IDL> print, inquad(x, y, x1, y1, x2, y2, x3, y3, x4, y4)
;
; On a sphere see
; <pro>clickincell</pro> ...
;
; @history
;      Sebastien Masson (smasson\@lodyc.jussieu.fr)
;      August 2003
;      Based on Convert_clic_ij.pro written by Gurvan Madec
;
; @version
; $Id$
;
;-
FUNCTION inquad, x, y, x1, y1, x2, y2, x3, y3, x4, y4 $
               , ONSPHERE=onsphere,  DOUBLE=double, DELTA=delta $
               , NOPRINT=noprint, NEWCOORD=newcoord
;
  compile_opt idl2, strictarrsubs
;
  ntofind = n_elements(x)
  nquad = n_elements(x2)
;
  IF keyword_set(onsphere) THEN BEGIN
; save the inputs parameters
    xin = x
    yin = y
    x1in = x1
    y1in = y1
    x2in = x2
    y2in = y2
    x3in = x3
    y3in = y3
    x4in = x4
    y4in = y4
; for map_set
    x = x MOD 360
    x1 = x1 MOD 360
    x2 = x2 MOD 360
    x3 = x3 MOD 360
    x4 = x4 MOD 360
; save !map
    save = {map:!map, x:!x, y:!y, z:!z, p:!p}
; do a satellite projection
    IF NOT keyword_set(delta) THEN delta = 4
    map_set, y[0], x[0], 0, /stereo, limit = [y[0]-delta/2, x[0]-delta/2, y[0]+delta/2, x[0]+delta/2], /noerase, /iso, /noborder
; use normal coordinates to reject cells which are out of the projection.
    tmp  = convert_coord(x, y, /DATA, /TO_NORMAL, DOUBLE = double)
    tmp1 = convert_coord(x1, y1, /DATA, /TO_NORMAL, DOUBLE = double)
    tmp2 = convert_coord(x2, y2, /DATA, /TO_NORMAL, DOUBLE = double)
    tmp3 = convert_coord(x3, y3, /DATA, /TO_NORMAL, DOUBLE = double)
    tmp4 = convert_coord(x4, y4, /DATA, /TO_NORMAL, DOUBLE = double)
; remove cell which have one corner with coordinates equal to NaN
    test = finite(tmp1[0, *]+tmp1[1, *]+tmp2[0, *]+tmp2[1, *] $
                  +tmp3[0, *]+tmp3[1, *]+tmp4[0, *]+tmp4[1, *])
    good = where(temporary(test) EQ 1)
;
    IF good[0] EQ -1 THEN BEGIN
      IF NOT keyword_set(noprint) THEN BEGIN
         ras = report('The point is out of the cells')
      ENDIF
; restore the input parameters
      x = temporary(xin)
      y = temporary(yin)
      x1 = temporary(x1in)
      y1 = temporary(y1in)
      x2 = temporary(x2in)
      y2 = temporary(y2in)
      x3 = temporary(x3in)
      y3 = temporary(y3in)
      x4 = temporary(x4in)
      y4 = temporary(y4in)
; restore old !map...
      !map = save.map
      !x = save.x
      !y = save.y
      !z = save.z
      !p = save.p
      RETURN,  -1
    ENDIF
;
    x  = tmp[0]
    y  = tmp[1]
    x1 = tmp1[0, good]
    y1 = tmp1[1, good]
    x2 = tmp2[0, good]
    y2 = tmp2[1, good]
    x3 = tmp3[0, good]
    y3 = tmp3[1, good]
    x4 = tmp4[0, good]
    y4 = tmp4[1, good]
;
    tmp1 = -1 & tmp2 = -1 & tmp3 = -1 & tmp4 = -1
; remove cells which are obviously bad
    test = (x1 GT x AND x2 GT x AND x3 GT x AND x4 GT x) $
      OR (x1 LT x AND x2 LT x AND x3 LT x AND x4 LT x) $
      OR (y1 GT y AND y2 GT y AND y3 GT y AND y4 GT y) $
      OR (y1 LT y AND y2 LT y AND y3 LT y AND y4 LT y)
    good2 = where(temporary(test) EQ 0)
;
    IF good2[0] EQ -1 THEN BEGIN
      IF NOT keyword_set(noprint) THEN BEGIN
         ras = report('The point is out of the cells')
      ENDIF
; restore the input parameters
      x = temporary(xin)
      y = temporary(yin)
      x1 = temporary(x1in)
      y1 = temporary(y1in)
      x2 = temporary(x2in)
      y2 = temporary(y2in)
      x3 = temporary(x3in)
      y3 = temporary(y3in)
      x4 = temporary(x4in)
      y4 = temporary(y4in)
; restore old !map...
      !map = save.map
      !x = save.x
      !y = save.y
      !z = save.z
      !p = save.p
      RETURN,  -1
    ENDIF
;
    nquad = n_elements(good2)
    x1 = x1[good2]
    y1 = y1[good2]
    x2 = x2[good2]
    y2 = y2[good2]
    x3 = x3[good2]
    y3 = y3[good2]
    x4 = x4[good2]
    y4 = y4[good2]
  ENDIF
;
;
; the point is inside the quadrilateral if test eq 1
; with test equal to:
;     test = ((x-x1)*(y2-y1) GE (x2-x1)*(y-y1)) $
;       *((x-x2)*(y3-y2) GT (x3-x2)*(y-y2)) $
;       *((x-x3)*(y4-y3) GT (x4-x3)*(y-y3)) $
;       *((x-x4)*(y1-y4) GE (x1-x4)*(y-y4))
;
; computation of test without any do loop for ntofind points (x,y) and
; nquad quadrilateral((x1,x2,x3,x4),(y1,y2,y3,y4))
; test dimensions are (ntofind, nquad)
; column i of test corresponds to the intersection of point i with all
; quadrilateral.
; row j of test corresponds to all the points localized in cell j

  IF keyword_set(double) THEN one = 1.d ELSE one = 1. 
  nquad_1 = replicate(one, nquad)
  ntofind_1 = replicate(one, ntofind)
  x_nquad = x[*]#replicate(one, nquad)
  y_nquad = y[*]#replicate(one, nquad)

  test = $
;  (x-x1)*(y2-y1) GE (x2-x1)*(y-y1)
  ( (x_nquad - ntofind_1#x1[*]) * (ntofind_1#(y2-y1)[*]) ) GE ( (ntofind_1#(x2-x1)[*]) * (y_nquad - ntofind_1#y1[*]) ) AND $
; *(x-x2)*(y3-y2) GE (x3-x2)*(y-y2)
  ( (x_nquad - ntofind_1#x2[*]) * (ntofind_1#(y3-y2)[*]) ) GE ( (ntofind_1#(x3-x2)[*]) * (y_nquad - ntofind_1#y2[*]) ) AND $
; *(x-x3)*(y4-y3) GE (x4-x3)*(y-y3)
  ( (x_nquad - ntofind_1#x3[*]) * (ntofind_1#(y4-y3)[*]) ) GE ( (ntofind_1#(x4-x3)[*]) * (y_nquad - ntofind_1#y3[*]) ) AND $
; *(x-x4)*(y1-y4) GE (x1-x4)*(y-y4)
  ( (x_nquad - ntofind_1#x4[*]) * (ntofind_1#(y1-y4)[*]) ) GE ( (ntofind_1#(x1-x4)[*]) * (y_nquad - ntofind_1#y4[*]) )
;
  nquad_1 = 1 & ntofind_1 = 1 & x_nquad = 1 & y_nquad = 1   ; free memory
;
; check test if ntofind gt 1
; if ntofind gt 1, each point must be localised in one uniq cell.
  IF ntofind GT 1 THEN BEGIN
; each column of test must have only 1 position equal to one
    chtest = total(test, 2)
; points out of the cells
    IF (where(chtest EQ 0))[0] NE -1 THEN BEGIN
      IF NOT keyword_set(noprint) THEN BEGIN
         ras = report('Points number '+strjoin(strtrim(where(chtest EQ 0), 1), ', ')+' are out of the grid')
      ENDIF
      stop
    ENDIF
; points in more than one cell
    IF (where(chtest GT 1))[0] NE -1 THEN BEGIN
      IF NOT keyword_set(noprint) THEN BEGIN
        ras = report('Points number '+strjoin(strtrim(where(chtest GT 1), 1), ', ')+' are in more than one cell')
      ENDIF
      stop
    ENDIF
  ENDIF
; find the points for which test eq 1
  found = where(temporary(test) EQ 1)
; if ntofind eq 1, the point may be localised in more than one grid
; cell ou may also be out of the cells
  IF ntofind EQ 1 THEN BEGIN
    CASE 1 OF
      found[0] EQ -1:BEGIN
        IF NOT keyword_set(noprint) THEN BEGIN
           ras = report('The point is out of the cells')
        ENDIF
        IF keyword_set(onsphere) THEN BEGIN
; restore old !map...
          !map = save.map
          !x = save.x
          !y = save.y
          !z = save.z
          !p = save.p
        ENDIF
        return,  -1
      END
      n_elements(found) GT ntofind:BEGIN
        IF NOT keyword_set(noprint) THEN BEGIN
          ras = report('The point is in more than one cell')
        ENDIF
      END
      ELSE:
    ENDCASE
  ENDIF ELSE BEGIN
; if ntofind GT 1, found must be sorted
; i position of found. this corresponds to one x,y point
    forsort = found MOD ntofind
; j position of found. this corresponds to cell in which is one x,y
; point
    found = temporary(found)/ntofind
; found must be sorted according to forsort
    found = found[sort(forsort)]
  ENDELSE
;
  IF keyword_set(onsphere) THEN BEGIN
    IF arg_present(newcoord) THEN BEGIN
      found = found[0]
      newcoord = [[x1[found], y1[found]] $
                  , [x2[found], y2[found]] $
                  , [x3[found], y3[found]] $
                  , [x4[found], y4[found]] $
                  , [x, y]]
    ENDIF
;
    found = good[good2[found]]
; restore the input parameters
    x = temporary(xin)
    y = temporary(yin)
    x1 = temporary(x1in)
    y1 = temporary(y1in)
    x2 = temporary(x2in)
    y2 = temporary(y2in)
    x3 = temporary(x3in)
    y3 = temporary(y3in)
    x4 = temporary(x4in)
    y4 = temporary(y4in)
; restore old !map...
    !map = save.map
    !x = save.x
    !y = save.y
    !z = save.z
    !p = save.p
  ENDIF
;
  RETURN, found
END