;;;; geopoint-utils.scm
;;;; Kon Lovett, May '17

(module geopoint-utils

(;export
  geopoint-in-closed-polygon? )

(import scheme)

(import chicken)

(use numbers)

(use geopoint geopolygon)

(use type-checks)

;;;

;;

;; Is the point inside a closed polygon?
;; Ray casting algorithm (http://en.wikipedia.org/wiki/Point_in_polygon)
;; (http://alienryderflex.com/polygon/)
(define (geopoint-in-closed-polygon? gp gpoly)
  ;
  (define (intersects? pi pj lat lon)
    (and
      (or
        (and (< (geopoint-latitude pi) lat) (<= lat (geopoint-latitude pj) ) )
        (and (< (geopoint-latitude pj) lat) (<= lat (geopoint-latitude pi)) ) )
      (or
        (<= (geopoint-longitude pi) lon)
        (<= (geopoint-longitude pj) lon) )
      (>=
        lon
        (+
          (geopoint-longitude pi)
          (*
            (/
              (- lat (geopoint-latitude pi))
              (- (geopoint-latitude pj) (geopoint-latitude pi)))
            (- (geopoint-longitude pj) (geopoint-longitude pi))))) ) )
  ;
	; test for intersection of ray with every segment of the polygon
	; start with the "closing" segment, then every [i-i],[i] segment
  ;
	(let ((gpoly (if (list? gpoly) (list->vector gpoly) gpoly)))
	  (check-geopolygon 'geopoint-in-closed-polygon? gpoly)
	  (check-geopoint 'geopoint-in-closed-polygon? gp)
    (let ((len (vector-length gpoly) )
          (lat (geopoint-latitude gp) )
          (lon (geopoint-longitude gp) ) )
      ; assumes an open-poly is "closed" so a closed-poly must be treated as "open"
      (let ((len (if (geopolygon-closed? gpoly) (fx- len 1) len)))
        (let loop ((i 0) (j (fx- len 1)) (poly? #f))
          (if (fx= i len)
            poly?
            (let ((pi (vector-ref gpoly i))
                  (pj (vector-ref gpoly j)) )
              (loop (fx+ i 1) i (if (intersects? pi pj lat lon) (not poly?) poly?)) ) ) ) ) ) ) )

) ;module geopoint-utils