;; An implementation of parametric curves.
;; 
;; This code is inspired by the Haskell rsagl library.
;;
;; Copyright 2012 Ivan Raikov and the Okinawa Institute of
;; Science and Technology.
;;
;; This program is free software: you can redistribute it and/or
;; modify it under the terms of the GNU General Public License as
;; published by the Free Software Foundation, either version 3 of the
;; License, or (at your option) any later version.
;;
;; This program is distributed in the hope that it will be useful, but
;; WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
;; General Public License for more details.
;;
;; A full copy of the GPL license can be found at
;; <http://www.gnu.org/licenses/>.
;;

(module parametric-curve

	(
	 parametric-curve?
	 simple-curve linear-curve 
	 translate-curve scale-curve
         sample-curve inverse-sample-curve
	 sample-curve* inverse-sample-curve*
	 iterate-curve range-curve bbox-curve
	 )
	
	(import scheme chicken data-structures)
  
	(require-library srfi-1 srfi-4 extras bvsp-spline)

	(import (only srfi-1 fold every list-tabulate zip concatenate)
		(only srfi-4 f64vector make-f64vector list->f64vector f64vector->list f64vector-length f64vector-ref f64vector-set!)
		(only extras fprintf pp)
		(prefix bvsp-spline bvsp-spline:)
		)

;; A boundary-value preserving spline 
(define-record-type spline (make-spline n k x y d d2)
  spline? 
  (n spline-n)
  (k spline-k)
  (x spline-x)
  (y spline-y)
  (d spline-d)
  (d2 spline-d2)
  )


;; A parametric equation with an associated interval, interpolating
;; spline, inverse interpolating spline.
(define-record-type peq (make-peq fn xmin xmax ymin ymax spl invspl)
  peq? 
  (fn      peq-fn )
  (xmin    peq-xmin)
  (xmax    peq-xmax)
  (ymin    peq-ymin)
  (ymax    peq-ymax)
  (spl     peq-spl)
  (invspl  peq-invspl)
  )


;; A N-dimensional parametric curve is a set of N parametric equations.
(define (parametric-curve? x) (every peq? x))


;; Defines a simple parametric equation.
(define (simple-peq n k fn xmin xmax)
  (if (> xmin xmax)
      (simple-peq n k fn xmax xmin)
      (let ((dx (/ (- xmax xmin) n)))
	(let* ((x (list-tabulate (+ 1 n) (lambda (i) (+ xmin (* i dx)))))
	       (y (map fn x))
	       (ymin (fold min +inf.0 y))
	       (ymax (fold max -inf.0 y))
	       (xv (list->f64vector x))
	       (yv (list->f64vector y)))

	  (let ((spl
		 (let-values (((d d2 constr errc diagn)
			       (bvsp-spline:compute n k xv yv)))
		   (if (zero? errc)
		       (make-spline n k xv yv d d2) 
		       (error 'simple-peq "unable to compute interpolating spline"))
		   ))

		(invspl
		 (let-values (((inv-d inv-d2 constr errc diagn)
			       (bvsp-spline:compute n k yv xv)))
		   (if (zero? errc)
		       (make-spline n k yv xv inv-d inv-d2)
		       (let ((y1 (car y)))
			 (if (every (lambda (y) (< (abs (- y y1)) 1e-12)) y) (car x)
			     (error 'simple-peq "unable to compute inverse interpolating spline" errc))
			 ))
		   ))
		)

	      (make-peq fn xmin xmax ymin ymax spl invspl)
	      ))
	  ))
      )


;; Samples a parametric equation at the given point of interest.
(define (sample-peq c)
  (let ((s (peq-spl c))
	(min (peq-xmin c))
	(max (peq-xmax c)))
    (lambda (xp)
      (and (>= xp min) (<= xp max)
	   (let ((v (bvsp-spline:evaluate (spline-n s) (spline-k s) 
					  (spline-x s) (spline-y s) 
					  (spline-d s) (spline-d2 s) 
					  (f64vector xp) 0)))
	     (and v (f64vector-ref v 0)))
	   ))
    ))


;; Inverse sampling: y -> x
(define (inverse-sample-peq c)
  (let (
	(s (peq-invspl c))
	(min (peq-ymin c))
	(max (peq-ymax c))
	)

    (if (number? s)
	(lambda (yp) (and (>= yp min) (<= yp max) s))
	(lambda (yp)
	  (and (>= yp min) (<= yp max)
	       (let ((v (bvsp-spline:evaluate (spline-n s) (spline-k s) 
					      (spline-x s) (spline-y s) 
					      (spline-d s) (spline-d2 s) 
					      (f64vector yp) 0)))
		 (and v (f64vector-ref v 0)))
	       ))
	)
    ))


;; Samples a parametric equation at the given points of interest.
(define (sample-peq* c)
  (let ((s (peq-spl c)))
    (lambda (xps)
      (let ((xpsv (if (list? xps) (list->f64vector xps) xps)))
	(bvsp-spline:evaluate (spline-n s) (spline-k s) 
			      (spline-x s) (spline-y s) 
			      (spline-d s) (spline-d2 s) 
			      xpsv 0))
      )))


;; Inverse sampling: [y] -> [x]
(define (inverse-sample-peq* c)
  (let ((s (peq-invspl c))
	(min (peq-ymin c))
	(max (peq-ymax c)))

    (if (number? s) 

	(lambda (yps) 
	  (let ((ypsv (if (list? yps) (list->f64vector yps) yps)))
	    (let* ((n (f64vector-length ypsv))
		   (r (make-f64vector n)))
	      (let recur ((i 0))
		(if (< i n) 
		    (let ((yv (f64vector-ref ypsv i)))
		      (if (and (>= yv min) (<= yv max)) 
			  (f64vector-set! r i s)
			  (f64vector-set! r i +inf.0))
		      (recur (+ 1 i)))
		    r))
	      )))
		  
	  (lambda (yps)
	    (let ((ypsv (if (list? yps) (list->f64vector yps) yps)))
	      (bvsp-spline:evaluate (spline-n s) (spline-k s) 
				    (spline-x s) (spline-y s) 
				    (spline-d s) (spline-d2 s) 
				    ypsv 0)))
	  )))


;; Samples a parametric equation at regular intervals in the range xmin..xmax inclusive.
(define (iterate-peq c n)
  (let* ((f   (sample-peq* c))
	 (xmin (peq-xmin c))
	 (xmax (peq-xmax c))
	 (dx   (if (zero? (- xmax xmin)) 0 (/ (- xmax xmin) (- n 1)))))
    (f (list-tabulate n (lambda (i) (+ xmin (* dx i)))))
    ))
  
;; Sampling of a parametric equation at regular intervals in the given range.
(define (range-peq c n xmin xmax)
  (if (> xmin xmax)
      (range-peq c n xmax xmin)
      (let (
	    (f1   (inverse-sample-peq* c))
	    (f2   (sample-peq* c))
	    (dx   (if (zero? (- xmax xmin)) 0 (/ (- xmax xmin) (- n 1))))
	    )
	(let* (
	       (v (f1 (list-tabulate n (lambda (i) (+ xmin (* dx i))))))
	       (u (f2 (f64vector->list v)))
	       )
	  (cons v u) 
	  ))
      ))
  

;; Transforms an equation using the given function.
(define (map-peq f x)
  (let ((fx (sample-peq x)) (splx (peq-spl x)))
    (let ((f1 (lambda (u) (f (fx u)))))
      (simple-peq (spline-n splx) (spline-k splx) f1 (peq-xmin x) (peq-xmax x))
      )))


;; Combines two equations using the given function.
(define (zip-peq f x y)
  (let ((fx   (sample-peq x)) 
	(splx (peq-spl x)) 
	(fy   (sample-peq y)))
    (let ((f1 (lambda (u) (f (fx u) (fy u)))))
      (simple-peq (spline-n splx) (spline-k splx) f1 (peq-xmin x) (peq-xmax x))
      )))


;; Translates a parametric equation.
(define (translate-peq x p) (map-peq (lambda (v) (+ x v)) p))


;; Scales a parametric equation.  
(define (scale-peq s p) (map-peq (lambda (v) (* s v)) p))


;; Defines a simple parametric curve.
(define (simple-curve n k fs tmin tmax)
  (if (null? fs) '()
      (let ((c (simple-peq n k (car fs) tmin tmax)))
	(cons c (simple-curve n k (cdr fs) tmin tmax))
	)))
  
  
;; Samples a curve at the given point
(define (sample-curve s)
  (let ((scs (map sample-peq s)))
    (lambda (t) (map (lambda (sc) (sc t)) scs))
    ))


;; Inverse curve sample
(define (inverse-sample-curve s)
  (let ((scs (map inverse-sample-peq s)))
    (lambda (t) (map (lambda (sc) (sc t)) scs))
    ))


;; Samples a curve at the given points
(define (sample-curve* s)
  (let ((scs (map sample-peq* s)))
    (lambda (ts) (map (lambda (sc) (sc ts)) scs))
    ))


;; Inverse curve sample
(define (inverse-sample-curve* s)
  (let ((scs (map sample-peq* s)))
    (lambda (ts) (map (lambda (sc) (sc ts)) scs))
    ))


;; Linear curve of the form c1 * x + c2
;; Argument coeffs supplies c1 and c2 for the different dimensions
(define (linear-curve n k coeffs tmin tmax )
  (simple-curve n k (map (lambda (s) (lambda (x) (+ (* (car s) x) (cadr s)))) coeffs) tmin tmax))


;; Translates a parametric curve.
(define (translate-curve xs c)
  (map (lambda (x p) (translate-peq x p)) xs c))

			 
;; Scales a parametric curve.
(define (scale-curve xs c)
  (map (lambda (x p) (scale-peq x p)) xs c))

			 
;; Obtain the bounding box of a curve
(define (bbox-curve c) (append (map peq-ymin c) (map peq-ymax c)))


;; Samples a parametric curve at regular intervals in the range xmin..xmax inclusive.
(define (iterate-curve c n) (map (lambda (p) (iterate-peq p n)) c))
  

;; Samples a parametric curve at regular intervals in the given ranges.
(define (range-curve c n rs) (map (lambda (p r) (range-peq p n (car r) (cadr r))) c rs))
  

)