;; An implementation of parametric curves.
;; 
;; This code is inspired by the Haskell rsagl library.
;;
;; Copyright 2012-2015 Ivan Raikov
;;
;; 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 line-segment
	 map-curve translate-curve scale-curve
         compose-curve
         sample-curve sample-curve*
	 iterate-curve fold-curve foldi-curve
         bbox-curve arc-length
	 )
	
	(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!
                      list->u32vector u32vector-ref)
		(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 and interpolating spline
(define-record-type peq (make-peq fn xmin xmax ymin ymax spl)
  peq? 
  (fn      peq-fn )
  (xmin    peq-xmin)
  (xmax    peq-xmax)
  (ymin    peq-ymin)
  (ymax    peq-ymax)
  (spl     peq-spl)
  )


(define-record-printer (peq x out)
  (fprintf out "#(peq fn=~A xmin=~A xmax=~A ymin=~A ymax=~A spl=~A)"
           (peq-fn x)
           (peq-xmin x)
           (peq-xmax x)
           (peq-ymin x)
           (peq-ymax x)
           (peq-spl x)
           ))



;; 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"))
		   ))
                )

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


;; 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)))
	   ))
    ))


;; 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)))
        (let ((res (bvsp-spline:evaluate (spline-n s) (spline-k s) 
                                         (spline-x s) (spline-y s) 
                                         (spline-d s) (spline-d2 s) 
                                         xpsv 0)))
          res
          ))
      ))
  )



;; 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))
         (delta (- xmax xmin))
	 (dx    (if (zero? delta) 0
                    (if (< n 2)
                        (error 'iterate-peq "number of iterations must be >= 2")
                        (/ (- xmax xmin) (- n 1)))))
         )
    (f (list-tabulate n (lambda (i) (+ xmin (* dx i)))))
    ))



;; 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 (compose-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 c)
  (let ((scs (map sample-peq c)))
    (lambda (t) (map (lambda (sc) (sc t)) scs))
    ))

;; Samples a curve at the given points
(define (sample-curve* c)
  (let ((scs (map sample-peq* c)))
    (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 coeffs tmin tmax )
  (simple-curve n 1 (map (lambda (s) 
                           (let ((c1 (car s)) (c2 (cadr s)))
                             (cond ((and (zero? c2) (zero? c1)) (lambda (x) 0.0))
                                   ((zero? c1) (lambda (x) c2))
                                   (else (lambda (x) (+ (* c1 x) c2))))))
                           coeffs)
                tmin tmax))


;; Line segment curve of the form (x1,xn) defined on the parameter range 0.0 .. 1.0
(define (line-segment n coeffs )
  (simple-curve n 1 (map (lambda (s) (lambda (x) (* x s))) coeffs) 0.0 1.0))


;; Maps the given functions to the parametric curve.
(define (map-curve fs c)
  (map (lambda (f p) (map-peq f p)) fs c))

;; Composes the parametric curves using the given functions.
(define (compose-curve fs c1 c2)
  (map (lambda (f p1 p2) (compose-peq f p1 p2)) fs c1 c2))


;; 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))


;; Folds a parametric curve at regular intervals in the range xmin..xmax inclusive.
(define (fold-curve c n f init) 
  (let* (
         (gs     (map sample-peq* c))
	 (xmins  (map peq-xmin c))
	 (xmaxs  (map peq-xmax c))
         (deltas (map - xmaxs xmins))
	 (dxs    (map (lambda (xmin xmax delta)
                        (if (zero? delta) 0
                            (if (< n 2)
                                (error 'fold-curve "number of iterations must be >= 2")
                                (/ (- xmax xmin) (- n 1)))))
                      xmins xmaxs deltas))
         (inds   (map (lambda (xmin dx) 
                        (list-tabulate n (lambda (i) (+ xmin (* dx i)))))
                      xmins dxs))
         (vs     (map (lambda (g ind) (g ind)) gs inds))
         )
    (let recur ((i 0) (init init))
      (if (< i n)
          (let* ((vsi   (map (lambda (vect) (f64vector-ref vect i)) vs))
                 (init1 (f vsi init)))
            (recur (+ 1 i) init1))
          init))
    ))
          

;; Like fold-curve, but F is of the form F(I,V,INIT)

;; Folds a parametric curve at regular intervals in the range xmin..xmax inclusive.
(define (foldi-curve c n f init) 
  (let* (
         (gs     (map sample-peq* c))
	 (xmins  (map peq-xmin c))
	 (xmaxs  (map peq-xmax c))
         (deltas (map - xmaxs xmins))
	 (dxs    (map (lambda (xmin xmax delta)
                        (if (zero? delta) 0
                            (if (< n 2)
                                (error 'fold-curve "number of iterations must be >= 2")
                                (/ (- xmax xmin) (- n 1)))))
                      xmins xmaxs deltas))
         (inds   (map (lambda (xmin dx) 
                        (list-tabulate n (lambda (i) (+ xmin (* dx i)))))
                      xmins dxs))
         (vs     (map (lambda (g ind) (g ind)) gs inds))
         )
    (let recur ((i 0) (init init))
      (if (< i n)
          (let* ((vsi   (map (lambda (vect) (f64vector-ref vect i)) vs))
                 (init1 (f i vsi init)))
            (recur (+ 1 i) init1))
          init))
    ))
          


;; Computes the arc length of the parametric curve given step dx
(define (arc-length c dx)
  (let* ((n (inexact->exact (round (/ 1.0 dx))))
         (v (iterate-curve c n)))
    (let recur ((i 1) (l 0.0) 
                (s (map (lambda (x) (f64vector-ref x 0)) v)))
      (if (< i n)
          (let ((s1 (map (lambda (x) (f64vector-ref x i)) v)))
            (recur (+ 1 i)
                   (+ l (sqrt (fold (lambda (x x1 l) (+ l (expt (- x1 x) 2))) 0.0 s s1)))
                   s1))
          l))
    ))
            
    
  

)