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