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

;+
; NAME:inquad
;
; PURPOSE: to find if an (x,y) point is in a quadrilateral (x1,x2,x3,x4)
;
; CATEGORY:grid manipulation
;
; CALLING SEQUENCE:
;
;     res = inquad(x, y, x1, y1, x2, y2, x3, y3, x4, y4)
;
; INPUTS:
;
;     x,y: 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. 
;
;     x1, y1, x2, y2, x3, y3, x4, y4: the coordinates of the
;     quadrilateral given in the CLOCKWISE order. Scalar or array.
;
; KEYWORD PARAMETERS:
;
;    /DOUBLE: use double precision to perform the computation 
;
;    /ONSPHERE: to specify that the quadilateral are on a sphere and
;    that teir coordinates are longitude-latitude coordinates. In this
;    case, est-west periodicity, poles singularity and other pbs
;    related to longitude-latitude coordinates are managed
;    automatically. 
;
;    ZOOMRADIUS:the zoom (circle centred on the (x,y) with a radius of
;    zoomradius degree where we look for the the quadrilateral which;    contains the (x,y) point) used for the satellite projection
;    when /onsphere is activated. Default is 4 and seems to be the
;    minimum which can be used. Can be increase if the cell size is
;    larger than 5 degrees.
;   
;    /NOPRINT: to suppress the print messages.
;
; OUTPUTS:
;
;    res, a n element 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)
;
; COMMON BLOCKS:none
;
; SIDE EFFECTS:
;
; RESTRICTIONS: I think degenerated quadrilateral (e.g. flat of
; twisted) is not work. This has to be tested.
;
; EXAMPLE:
;
;       x = 1.*[1, 2, 6, 7, 3]
;       y = 1.*[1, 3, 3, 4, 7]
;       x1 = 1.*[0,4,2]
;       y1 = 1.*[1,4,8]
;       x2 = 1.*[1,6,4]
;       y2 = 1.*[5,6,8]
;       x3 = 1.*[3,8,4]
;       y3 = 1.*[4,4,6]
;       x4 = 1.*[2,6,2]
;       y4 = 1.*[0,2,6]
;       splot, [0,10], [0,10], xstyle = 1, ystyle = 1,/nodata
;       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]]
;       oplot, x, y, color = 20, psym = 1, thick = 2
;       print, inquad(x, y, x1, y1, x2, y2, x3, y3, x4, y4)
;
;      On a sphere see clickincell.pro...
;
; MODIFICATION HISTORY:
;      Sebastien Masson (smasson@lodyc.jussieu.fr)
;      August 2003
;      Based on Convert_clic_ij.pro written by Gurvan Madec 
;
;-
FUNCTION inquad, x, y, x1, y1, x2, y2, x3, y3, x4, y4, ONSPHERE = onsphere,  DOUBLE = double, ZOOMRADIUS = zoomradius, NOPRINT = noprint, NEWCOORD = newcoord
;
;
  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(zoomradius) THEN zoomradius = 4
    map_set, y[0], x[0], 0, /satellite, sat_p = [1+zoomradius*20/6371.229, 0, 0], /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 print, 'The point is out of the cells'
; 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 print, 'The point is out of the cells'
; 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 quadilateral 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 quadilateral((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
; quadirlateral. 
; row j of test corresponds to all the points localized in cell j 
  test = $
; (x-x1)
  ((x[*]#replicate(1, nquad)-replicate(1, ntofind)#x1[*]) $
; *(y2-y1)
  *(replicate(1, ntofind)#(y2-y1)[*]) $
; GE (x2-x1)
  GE ((replicate(1, ntofind)#(x2-x1)[*]) $
; *(y-y1)
  *(y[*]#replicate(1, nquad)-replicate(1, ntofind)#y1[*])))
;-------
  test = temporary(test) $
; *(x-x2)
  *((x[*]#replicate(1, nquad)-replicate(1, ntofind)#x2[*]) $
; *(y3-y2)
  *(replicate(1, ntofind)#(y3-y2)[*]) $
; GE (x3-x2)
  GE ((replicate(1, ntofind)#(x3-x2)[*]) $
; *(y-y2)
  *(y[*]#replicate(1, nquad)-replicate(1, ntofind)#y2[*])))
;-------
  test = temporary(test) $
; *(x-x3)
  *((x[*]#replicate(1, nquad)-replicate(1, ntofind)#x3[*]) $
; *(y4-y3)
  *(replicate(1, ntofind)#(y4-y3)[*]) $
; GE (x4-x3)
  GE ((replicate(1, ntofind)#(x4-x3)[*]) $
; *(y-y3)
  *(y[*]#replicate(1, nquad)-replicate(1, ntofind)#y3[*])))
;-------
  test = temporary(test) $
; *(x-x4)
  *((x[*]#replicate(1, nquad)-replicate(1, ntofind)#x4[*]) $
; *(y1-y4)
  *(replicate(1, ntofind)#(y1-y4)[*]) $
; GE (x1-x4)
  GE ((replicate(1, ntofind)#(x1-x4)[*]) $
; *(y-y4)
  *(y[*]#replicate(1, nquad)-replicate(1, ntofind)#y4[*])))
;
; 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 print, 'Points number '+strjoin(strtrim(where(chtest EQ 0), 1), ', ')+' are out of the grid'
      stop
    ENDIF
; points in more than one cell
    IF (where(chtest GT 1))[0] NE -1 THEN BEGIN
      IF NOT keyword_set(noprint) THEN print, 'Points number '+strjoin(strtrim(where(chtest GT 1), 1), ', ')+' are in more than one cell'
      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 print, 'The point is out of the cells'
        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 print, 'The point is in more than one cell'
      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 accordind to forsort
    found = found[sort(forsort)]
  ENDELSE
;
  IF keyword_set(onsphere) THEN BEGIN
    IF arg_present(newcoord) THEN BEGIN
      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