;;;; 3dr.lisp
 
 (in-package #:3dr)
 
 ;;; "3dr" goes here. Hacks and glory await!
 
 (defstruct (point (:type vector)) x y z)
 
 (defclass plane ()
   ((cells :initarg :cells :accessor cells :type (array integer (* *)))
    (distance :initarg :distance :accessor distance :type integer)))
 

 ;;(declaim (ftype (function ((vector integer 3) integer) (vector rational 3)) intercept-coordinate))

 (defun intercept-coordinate (point distance &optional out-point)

   (declare (inline) (optimize (speed 3)))

   (when (null out-point)

                                         ;(break)

     (setf out-point (vector 0 0 0)))
   (let ((a (elt point 0)) (b (elt point 1)) (c (elt point 2)))
     (setf (elt out-point 0) distance)
     (setf (elt out-point 1) (/ (* b distance) a))
     (setf (elt out-point 2) (/ (* c distance) a))
     out-point))
 
 (defun run-plane (plane point-cb)

   (declare (optimize (speed 3)))

   (let ((a (distance plane))
         (points (cells plane)))
     (destructuring-bind (b-bound c-bound) (array-dimensions points)
       (lparallel:pdotimes (b (1- b-bound))
         (let ((the-point (make-point :x a :y (- b (/ b-bound 2)) :z 0))
               (out-point (make-point :x 0 :y 0 :z 0)))
           (dotimes (c (1- c-bound))
             (setf (point-z the-point) (- c (/ c-bound 2)))
             (setf (aref points b c)
                   (funcall point-cb the-point out-point))))))))
 
 (declaim (ftype (function ((vector rational 3) (vector rational 3)) float) c-distance))
 (defun c-distance (point-1 point-2)
   ;(declare (optimize (speed 0) (safety 3) (debug 3)))
   (when (not (and (point-x point-1) (point-x point-2)))
     (break))
 
   ;(declare (inline) (optimize (speed 3) (safety 1) (debug 0) (space 0)))
   (sqrt (+ (expt (- (aref point-2 0) (aref point-1 0)) 2)
            (expt (- (aref point-2 1) (aref point-1 1)) 2)
            (expt (- (aref point-2 2) (aref point-1 2)) 2))))
 
 (defun extract-projected-point (point obj-dist &optional out-point)
   ;(declare (inline) (optimize (speed 3)))
   (intercept-coordinate point obj-dist out-point))
 
 (defclass shape ()
   ((origin :initarg :origin :accessor origin)
    (radius :initarg :radius :accessor radius)
    (args :initarg :args :accessor args)
    (render-func :initarg :render-func :accessor render-func)))
 
 (defgeneric get-render-cb (shape))
 (defmethod get-render-cb ((shape shape))
   (with-slots (origin radius render-func args) shape
     (apply render-func origin radius args)))
 
 (defun make-neg-sphere-cb (origin radius)
   (let ((obj-dist (point-x origin))
         (r_sqr (expt radius 2))) 
     (lambda (point &optional out-point extracted-point)
       ;(declare (inline) (optimize (speed 3)))
       (let* ((point (or extracted-point (extract-projected-point point obj-dist out-point)))
              (d_o (c-distance origin point)))
         (if (<= d_o radius)
           (* -1 (ceiling (sqrt (- r_sqr (expt d_o 2))))) 
           0)))))
 
 ;; d_o + r : upright :: upright : r - d_o
 ;; upright^2 = r^2 - d_o^2
 ;; upright = sqrt ( r^2 - d_o^2 )
 (defun make-sphere-cb (origin radius)
   (let ((obj-dist (point-x origin))
         (r_sqr (expt radius 2)))
     (lambda (point &optional out-point extracted-point)

       (declare (inline) (optimize (speed 3)))

       (let* ((out-point (or out-point (make-point :x 0 :y 0 :z 0)))
              (point (or extracted-point (extract-projected-point point obj-dist out-point))))
         (let ((d_o (c-distance origin point)))
           (values (if (<= d_o radius)

                       (ceiling (* (/ (sqrt (- r_sqr
                                               (expt d_o
                                                     2)))
                                      radius)
                                   #. (expt 2 7))) 
                       0)

                   point))))))
 

 (defun make-interesting-cb (origin radius)
   (let ((obj-dist (point-x origin))
         (r_sqr (expt radius 3)))
     (lambda (point &optional out-point extracted-point)
       (declare (inline) (optimize (speed 3)))
       (let* ((out-point (or out-point (make-point :x 0 :y 0 :z 0)))
              (point (or extracted-point (extract-projected-point point obj-dist out-point))))
         (let ((d_o (c-distance origin point)))
           (values (if (<= d_o radius)
                       (ceiling (* (/ (expt (- r_sqr
                                               (expt d_o
                                                     3))
                                            0.33)
                                      radius)
                                   #. (expt 2 7))) 
                       0)
                   point))))))
 

 (defun make-donut-cb (origin radius width)
   (let ((pos-cb (make-sphere-cb origin radius))
         (neg-cb (make-neg-sphere-cb origin (- radius width)))
         (obj-dist (point-x origin)))
     (lambda (point &optional out-point)

       (declare (inline) (optimize (speed 3)))

       (let* ((out-point (or out-point (make-point :x 0 :y 0 :z 0))) 
              (extracted-point (extract-projected-point point obj-dist out-point)))
         (+ (funcall pos-cb point out-point extracted-point)
            (funcall neg-cb point out-point extracted-point))))))
 
 (defun combine-cbs (&rest cbs)
   (lambda (point &optional out-point)
     (let ((out-point (or out-point (make-point :x 0 :y 0 :z 0))))
       (loop for cb in cbs sum (funcall cb point out-point)))))
 
 (defun combine-shapes (&rest shapes)
   (let* ((shape-cbs (mapcar #'get-render-cb
                             (lparallel:psort shapes #'<
                                              :key (lambda (shape)
                                                     (with-slots (origin radius) shape
                                                       (-
                                                         (c-distance #(0 0 0)
                                                                     origin)
                                                         radius))))))
          (origins (mapcar 'origin shapes))
          (radii (mapcar 'radius shapes))
          (light-scaling-factor (lparallel:pmapcar
                                  (lambda (origin radius)
                                    (abs
                                      (/ (- (c-distance #(0 100 100)
                                                        origin)
                                            radius)
                                         100)))
                                  origins
                                  radii)))
     (lambda (point &optional out-point)
       (let ((out-point (or out-point (make-point :x 0 :y 0 :z 0))))
         (loop for shape-cb in shape-cbs
               for scale-factor in light-scaling-factor
               for val = (ceiling
                           (* scale-factor
                              (funcall shape-cb point out-point)))
               until (< 0 val)
               finally (return val))))))
 
 (defun array-to-pgm (arr)

   ;;(declare (optimize (speed 3)))

   (with-output-to-string (s)
     (destructuring-bind (height width) (array-dimensions arr)
       (format s "P2~%~3d ~3d~%~d~%"
               width height
               (1+ (reduce #'max (make-array (* width height) :displaced-to arr)))))
     (loop with (x-bound y-bound) = (array-dimensions arr)

        for x from 0 to (1- x-bound)
        do (loop for y from 0 to (1- y-bound)
              for val = (aref arr x y)
              do (princ val s)
              do (princ #\space s))
        do (terpri s))))

 
 (defmacro defshape (render-func origin radius &rest args)
   `(make-instance 'shape
                   :render-func (function ,render-func)
                   :origin (copy-seq ,origin)
                   :radius ,radius
                   :args (list ,@args)))
 
 (defun main (&rest args)
   (declare (ignore args))

   ;;(sb-profile:profile intercept-coordinate c-distance)

   (loop for img-num from 0 to 1
      do
        (let ((the-plane (make-instance 'plane
                                        :cells (make-array '(1024 1024))
                                        :distance 50))
              (lparallel:*kernel* (lparallel:make-kernel 7)))
 
          (labels ((100- (x) (- x 100))
                   (random-coord () (100- (random 200))))
            (time
             (run-plane the-plane
                        (apply #'combine-shapes
                               (loop for x from 0 to 100 by 20
                                  append (loop for y from 0 to 100 by 20
                                            collect
                                              (defshape make-sphere-cb
                                                  (vector 30 x y)
                                                10)))))))
 
          (with-open-file (s (format nil "/tmp/spheres.~d.pgm" img-num)
                             :direction :output :if-exists :supersede)
            (write-sequence (array-to-pgm (cells the-plane))
                            s))
          nil))
                                         ;(sb-profile:report)

   )
 
 ;(loop with (x-bound y-bound) = (array-dimensions (cells the-plane))
 ;      for x from 0 to (1- x-bound)
 ;      do (loop for y from 0 to (1- y-bound)
 ;               for val = (aref (cells the-plane) x y)
 ;               when (= val 0) do (princ #\space)
 ;               unless (= val 0) do (princ val)
 ;               do (princ #\space))
 ;      do (terpri))