youz/aobench.lisp

## aobench.lisp
(declaim (optimize (speed 3) (debug 0) (safety 0)))

(defparameter image-width  256)
(defparameter image-height 256)
(defparameter nsubsamples 2)
(defparameter nao-samples 8)

;; macros
(eval-when (:compile-toplevel :load-toplevel :execute)

  #+:xyzzy(require 'cmu_loop)
  #+:xyzzy(defmacro the (type sym) sym)

  (defmacro with-gensyms (syms &body body)
    `(let ,(mapcar #'(lambda (s) `(,s (gensym))) syms)
       ,@body))

  (defmacro mvlet ((&rest binds) &rest body)
    (let ((vars (caar binds))
          (values (cadar binds)))
      `(multiple-value-bind ,vars ,values
         ,@(if (cdr binds)
               `((mvlet ,(cdr binds) ,@body))
             body))))

  (defmacro defvaccessor (syms float)
    (flet ((defref (sym idx)
             (if float
                 `(defmacro ,sym (v) `(the double-float (svref ,v ,,idx)))
               `(defmacro ,sym (v) `(svref ,v ,,idx)))))
      `(progn
         ,@(loop
             for sym in syms
             for i from 0 below (length syms)
             collect (defref sym i)))))

  (defmacro dotimes% ((sym times) &body form)
    `(loop for ,sym from 0 below ,times do ,@form))


  ;; vector
  (defvaccessor (vx vy vz) t)

  (defmacro vv (v) `(values (vx ,v) (vy ,v) (vz ,v)))
  (defmacro sq (x) `(* ,x ,x))
  (defmacro len (x y z) `(sqrt (+ (sq ,x) (sq ,y) (sq ,z))))

  (defmacro vecv (form)
    `(multiple-value-call #'vector ,form))

  (defmacro vset (v form)
    (with-gensyms (x y z)
      `(multiple-value-bind (,x ,y ,z) ,form
         (setf (vx ,v) ,x (vy ,v) ,y (vz ,v) ,z)
         ,v)))

  (defmacro vsetn (r &rest forms)
    `(setf ,@(loop
                for i from 0 below (length forms)
                for f in forms
                append `((svref ,r ,i) ,f))))

  (defmacro vcall (f &rest vectors)
    (let* ((binds (mapcar (lambda (form) (list (gensym) form)) vectors))
           (args (mapcan (lambda (b) (let ((s (car b))) `((vx ,s) (vy ,s) (vz, s)))) binds)))
      `(let ,binds (,f ,@args))))

  )

(defun vadd-v (ax ay az bx by bz c)
  (declare (type double-float ax ay az bx by bz c))
  (values (+ ax (* c bx))
          (+ ay (* c by))
          (+ az (* c bz))))

(defun vsub-v (ax ay az bx by bz)
  (declare (type double-float ax ay az bx by bz))
  (values (- ax bx) (- ay by) (- az bz)))

(defun vcross-v (ax ay az bx by bz)
  (declare (type double-float ax ay az bx by bz))
  (values (- (* ay bz) (* az by))
          (- (* az bx) (* ax bz))
          (- (* ax by) (* ay bx))))

(defun vdot-v (ax ay az bx by bz)
  (declare (type double-float ax ay az bx by bz))
  (+ (* ax bx) (* ay by) (* az bz)))


(defmacro vadd (a b c)
  `(vadd-v (vx ,a) (vy ,a) (vz ,a)
           (vx ,b) (vy ,b) (vz ,b) ,c))

(defmacro vsub (a b)
  `(vsub-v (vx ,a) (vy ,a) (vz ,a)
           (vx ,b) (vy ,b) (vz ,b)))

(defmacro vcross (a b)
  `(vcross-v (vx ,a) (vy ,a) (vz ,a)
             (vx ,b) (vy ,b) (vz ,b)))

(defmacro vdot (a b)
  `(vdot-v (vx ,a) (vy ,a) (vz ,a)
           (vx ,b) (vy ,b) (vz ,b)))

#+:nil
(defun vnormalize (a)
  (let ((d (len (vx a) (vy a) (vz a))))
    (when (> d 1d-17)
      (vset a (values (/ (vx a) d) (/ (vy a) d) (/ (vz a) d))))
    a))

(defun normalize (x y z)
  (declare (type double-float x y z))
  (let ((d (len x y z)))
    (declare (type double-float d))
    (if (> d 1d-17)
        (values (/ x d) (/ y d) (/ z d))
      (values x y z))))


;; geometry

(defun make-isect ()
  (vector 1.0d30 nil
          (vector 0d0 0d0 0d0)
          (vector 0d0 0d0 0d0)))

(defvaccessor (rox roy roz rdx rdy rdz) t)
(defvaccessor (is-dist is-hit is-p is-n) nil)

(defmacro init-isect (i)
  `(setf (is-dist ,i) 1.0d30
         (is-hit ,i) nil
         (vx (is-p ,i)) 0d0
         (vy (is-p ,i)) 0d0
         (vz (is-p ,i)) 0d0
         (vx (is-n ,i)) 0d0
         (vy (is-n ,i)) 0d0
         (vz (is-n ,i)) 0d0))


(defun sphere (center radius)
  (declare (type (simple-vector 3) center)
           (type double-float radius))
  #'(lambda (ray isect)
      (declare (type simple-vector ray isect))
      (let* ((rsx (- (rox ray) (vx center)))
             (rsy (- (roy ray) (vy center)))
             (rsz (- (roz ray) (vz center)))
             (b (vdot-v rsx rsy rsz (rdx ray) (rdy ray) (rdz ray)))
             (d (+ (sq (the double-float b))
                   (- (the double-float (vdot-v rsx rsy rsz rsx rsy rsz)))
                   (* radius radius))))
        (declare (type double-float b d))
        (when (> d 0d0)
          (let ((dist (- (- b) (sqrt d))))
            (declare (type double-float dist))
            (when (< 0d0 dist (the double-float (is-dist isect)))
              (setf (is-dist isect) dist
                    (is-hit isect) t)
              (vset (is-p isect)
                    (vadd-v (rox ray) (roy ray) (roz ray)
                            (rdx ray) (rdy ray) (rdz ray) dist))
              (vset (is-n isect)
                    (multiple-value-call #'normalize
                      (vsub (is-p isect) center)))))))))

(defun plane (p n)
  #'(lambda (ray isect)
      (declare (type simple-vector ray isect))
      (let ((d (-(vdot p n)))
            (v (vdot-v (rdx ray) (rdy ray) (rdz ray)
                       (vx n) (vy n) (vz n))))
        (declare (type double-float d v))
        (when (> (abs v) 1d-17)
          (let ((dist (/ (- (+ (vdot-v (rox ray) (roy ray) (roz ray)
                                       (vx n) (vy n) (vz n))
                               d)) v)))
            (declare (type double-float dist))
            (when (< 0d0 dist (the double-float (is-dist isect)))
              (setf (is-dist isect) dist
                    (is-hit isect) t)
              (vset (is-n isect) (vv n))
              (vset (is-p isect)
                    (vadd-v (rox ray) (roy ray) (roz ray)
                            (rdx ray) (rdy ray) (rdz ray) dist))))))))

(defun defscene ()
  (list (sphere (vector -2d0 0d0 -3.5d0) 0.5d0)
        (sphere (vector -0.5d0 0d0 -3d0) 0.5d0)
        (sphere (vector 1d0 0d0 -2.2d0) 0.5d0)
        (plane (vector 0d0 -0.5d0 0d0) (vector 0d0 1d0 0d0))))


;; ao

(let ((b0 (vector 0d0 0d0 0d0))
      (b1 (vector 0d0 0d0 0d0))
      (b2 (vector 0d0 0d0 0d0))
      (isect (make-isect))
      (occ-isect (make-isect))
      (ray (vector 0d0 0d0 0d0 0d0 0d0 0d0))
      (occ-ray (vector 0d0 0d0 0d0 0d0 0d0 0d0)))

  (declare (type simple-vector b0 b1 b2 isect occ-isect ray occ-ray))

  (defun clamp (f)
    (declare (type double-float f))
    (let ((i (* f 255.5)))
      (cond ((< i 0) 0)
            ((> i 255) 255)
            (t (round i)))))

  (defun ortho-basis ()
    (let ((n (is-n isect)))
      (setf (vx b0) 0d0 (vy b0) 0d0 (vz b0) 0d0)
      (cond ((< -0.6 (vy n) 0.6) (setf (vy b0) 1d0))
            ((< -0.6 (vz n) 0.6) (setf (vz b0) 1d0))
            (t (setf (vx b0) 1d0)))
      (vset b0 (multiple-value-call #'normalize (vcross b0 n)))
      (vset b1 (multiple-value-call #'normalize (vcross n b0)))
      (setq b2 n)))

  (defun ambient-occlusion (scene)
    (let* ((ntheta nao-samples)
           (nphi nao-samples)
           (eps 1d-4)
           (occlusion 0))
      (declare (type fixnum ntheta nphi occlusion))
      (multiple-value-bind (px py pz) (vadd (is-p isect) (is-n isect) eps)
        (ortho-basis)
        (dotimes% (j nphi)
          (dotimes% (i ntheta)
            (let* ((r (random 1d0))
                   (phi (* 6.283185307179586d0 (random 1d0)))
                   (sqr (sqrt (- 1d0 (the double-float r)))))
              (declare (type double-float r phi sqr))
              (let ((x (* (cos phi) sqr))
                    (y (* (sin phi) sqr))
                    (z (sqrt r)))
                (declare (type double-float x y z))
                (vsetn occ-ray
                       px py pz
                       (+ (* x (vx b0)) (* y (vx b1)) (* z (vx b2)))
                       (+ (* x (vy b0)) (* y (vy b1)) (* z (vy b2)))
                       (+ (* x (vz b0)) (* y (vz b1)) (* z (vz b2))))
                (init-isect occ-isect)
                (dolist (f scene)
                  (funcall f occ-ray occ-isect))
                (when (is-hit occ-isect)
                  (incf occlusion)))))))
      (/ (- (* ntheta nphi) occlusion)
         (* ntheta nphi)
         1d0)))


  (defun render (scene w h nsubs)
    (declare (type fixnum w h nsubs))
    (let ((image (make-array (list w h) :element-type '(unsigned-byte 8))))
      (dotimes% (y h)
        (dotimes% (x w)
          (let ((rad 0d0))
            (declare (type double-float rad))
            ;; subsampling
            (dotimes% (v nsubs)
              (dotimes% (u nsubs)
                (let* ((px (/ (+ x (/ u (float nsubs)) (- (/ w 2d0))) (/ w 2d0)))
                       (py (- (/ (+ y (/ v (float nsubs)) (- (/ h 2d0))) (/ h 2d0)))))
                  (multiple-value-bind (rx ry rz) (normalize px py -1d0)
                    (vsetn ray 0d0 0d0 0d0 rx ry rz))
                  (init-isect isect)
                  (loop for f in scene do
                    (funcall f ray isect))
                  (when (is-hit isect)
                    (let ((col (ambient-occlusion scene)))
                      (incf rad (the double-float col)))))))
            (setf (aref image x y) (clamp (/ rad (sq nsubs)))))))
      image)))


;; write-file

(defun write-pnm (file image w h)
  (declare (type fixnum w h))
  (with-open-file (s file :direction :output :if-exists :overwrite :if-does-not-exist :create)
    (format s "P2~%~D ~D~%~D~%" w h 255)
    (dotimes% (y h)
      (dotimes% (x w)
        (let ((p (aref image x y)))
          (format s "~D~%" p))))
    t))


(defun main (&optional (fn "test2.pgm") (w image-width) (h image-height))
  (let ((img (render (defscene) w h nsubsamples)))
    (write-pnm fn img w h)))

(time (main))
	(declaim (optimize (speed 3) (debug 0) (safety 0)))

	(defparameter image-width 256)
	(defparameter image-height 256)
	(defparameter nsubsamples 2)
	(defparameter nao-samples 8)

	;; macros
	(eval-when (:compile-toplevel :load-toplevel :execute)

	#+:xyzzy(require 'cmu_loop)
	#+:xyzzy(defmacro the (type sym) sym)

	(defmacro with-gensyms (syms &body body)
	`(let ,(mapcar #'(lambda (s) `(,s (gensym))) syms)
	,@body))

	(defmacro mvlet ((&rest binds) &rest body)
	(let ((vars (caar binds))
	(values (cadar binds)))
	`(multiple-value-bind ,vars ,values
	,@(if (cdr binds)
	`((mvlet ,(cdr binds) ,@body))
	body))))

	(defmacro defvaccessor (syms float)
	(flet ((defref (sym idx)
	(if float
	`(defmacro ,sym (v) `(the double-float (svref ,v ,,idx)))
	`(defmacro ,sym (v) `(svref ,v ,,idx)))))
	`(progn
	,@(loop
	for sym in syms
	for i from 0 below (length syms)
	collect (defref sym i)))))

	(defmacro dotimes% ((sym times) &body form)
	`(loop for ,sym from 0 below ,times do ,@form))


	;; vector
	(defvaccessor (vx vy vz) t)

	(defmacro vv (v) `(values (vx ,v) (vy ,v) (vz ,v)))
	(defmacro sq (x) `(* ,x ,x))
	(defmacro len (x y z) `(sqrt (+ (sq ,x) (sq ,y) (sq ,z))))

	(defmacro vecv (form)
	`(multiple-value-call #'vector ,form))

	(defmacro vset (v form)
	(with-gensyms (x y z)
	`(multiple-value-bind (,x ,y ,z) ,form
	(setf (vx ,v) ,x (vy ,v) ,y (vz ,v) ,z)
	,v)))

	(defmacro vsetn (r &rest forms)
	`(setf ,@(loop
	for i from 0 below (length forms)
	for f in forms
	append `((svref ,r ,i) ,f))))

	(defmacro vcall (f &rest vectors)
	(let* ((binds (mapcar (lambda (form) (list (gensym) form)) vectors))
	(args (mapcan (lambda (b) (let ((s (car b))) `((vx ,s) (vy ,s) (vz, s)))) binds)))
	`(let ,binds (,f ,@args))))

	)

	(defun vadd-v (ax ay az bx by bz c)
	(declare (type double-float ax ay az bx by bz c))
	(values (+ ax (* c bx))
	(+ ay (* c by))
	(+ az (* c bz))))

	(defun vsub-v (ax ay az bx by bz)
	(declare (type double-float ax ay az bx by bz))
	(values (- ax bx) (- ay by) (- az bz)))

	(defun vcross-v (ax ay az bx by bz)
	(declare (type double-float ax ay az bx by bz))
	(values (- (* ay bz) (* az by))
	(- (* az bx) (* ax bz))
	(- (* ax by) (* ay bx))))

	(defun vdot-v (ax ay az bx by bz)
	(declare (type double-float ax ay az bx by bz))
	(+ (* ax bx) (* ay by) (* az bz)))


	(defmacro vadd (a b c)
	`(vadd-v (vx ,a) (vy ,a) (vz ,a)
	(vx ,b) (vy ,b) (vz ,b) ,c))

	(defmacro vsub (a b)
	`(vsub-v (vx ,a) (vy ,a) (vz ,a)
	(vx ,b) (vy ,b) (vz ,b)))

	(defmacro vcross (a b)
	`(vcross-v (vx ,a) (vy ,a) (vz ,a)
	(vx ,b) (vy ,b) (vz ,b)))

	(defmacro vdot (a b)
	`(vdot-v (vx ,a) (vy ,a) (vz ,a)
	(vx ,b) (vy ,b) (vz ,b)))

	#+:nil
	(defun vnormalize (a)
	(let ((d (len (vx a) (vy a) (vz a))))
	(when (> d 1d-17)
	(vset a (values (/ (vx a) d) (/ (vy a) d) (/ (vz a) d))))
	a))

	(defun normalize (x y z)
	(declare (type double-float x y z))
	(let ((d (len x y z)))
	(declare (type double-float d))
	(if (> d 1d-17)
	(values (/ x d) (/ y d) (/ z d))
	(values x y z))))


	;; geometry

	(defun make-isect ()
	(vector 1.0d30 nil
	(vector 0d0 0d0 0d0)
	(vector 0d0 0d0 0d0)))

	(defvaccessor (rox roy roz rdx rdy rdz) t)
	(defvaccessor (is-dist is-hit is-p is-n) nil)

	(defmacro init-isect (i)
	`(setf (is-dist ,i) 1.0d30
	(is-hit ,i) nil
	(vx (is-p ,i)) 0d0
	(vy (is-p ,i)) 0d0
	(vz (is-p ,i)) 0d0
	(vx (is-n ,i)) 0d0
	(vy (is-n ,i)) 0d0
	(vz (is-n ,i)) 0d0))


	(defun sphere (center radius)
	(declare (type (simple-vector 3) center)
	(type double-float radius))
	#'(lambda (ray isect)
	(declare (type simple-vector ray isect))
	(let* ((rsx (- (rox ray) (vx center)))
	(rsy (- (roy ray) (vy center)))
	(rsz (- (roz ray) (vz center)))
	(b (vdot-v rsx rsy rsz (rdx ray) (rdy ray) (rdz ray)))
	(d (+ (sq (the double-float b))
	(- (the double-float (vdot-v rsx rsy rsz rsx rsy rsz)))
	(* radius radius))))
	(declare (type double-float b d))
	(when (> d 0d0)
	(let ((dist (- (- b) (sqrt d))))
	(declare (type double-float dist))
	(when (< 0d0 dist (the double-float (is-dist isect)))
	(setf (is-dist isect) dist
	(is-hit isect) t)
	(vset (is-p isect)
	(vadd-v (rox ray) (roy ray) (roz ray)
	(rdx ray) (rdy ray) (rdz ray) dist))
	(vset (is-n isect)
	(multiple-value-call #'normalize
	(vsub (is-p isect) center)))))))))

	(defun plane (p n)
	#'(lambda (ray isect)
	(declare (type simple-vector ray isect))
	(let ((d (-(vdot p n)))
	(v (vdot-v (rdx ray) (rdy ray) (rdz ray)
	(vx n) (vy n) (vz n))))
	(declare (type double-float d v))
	(when (> (abs v) 1d-17)
	(let ((dist (/ (- (+ (vdot-v (rox ray) (roy ray) (roz ray)
	(vx n) (vy n) (vz n))
	d)) v)))
	(declare (type double-float dist))
	(when (< 0d0 dist (the double-float (is-dist isect)))
	(setf (is-dist isect) dist
	(is-hit isect) t)
	(vset (is-n isect) (vv n))
	(vset (is-p isect)
	(vadd-v (rox ray) (roy ray) (roz ray)
	(rdx ray) (rdy ray) (rdz ray) dist))))))))

	(defun defscene ()
	(list (sphere (vector -2d0 0d0 -3.5d0) 0.5d0)
	(sphere (vector -0.5d0 0d0 -3d0) 0.5d0)
	(sphere (vector 1d0 0d0 -2.2d0) 0.5d0)
	(plane (vector 0d0 -0.5d0 0d0) (vector 0d0 1d0 0d0))))


	;; ao

	(let ((b0 (vector 0d0 0d0 0d0))
	(b1 (vector 0d0 0d0 0d0))
	(b2 (vector 0d0 0d0 0d0))
	(isect (make-isect))
	(occ-isect (make-isect))
	(ray (vector 0d0 0d0 0d0 0d0 0d0 0d0))
	(occ-ray (vector 0d0 0d0 0d0 0d0 0d0 0d0)))

	(declare (type simple-vector b0 b1 b2 isect occ-isect ray occ-ray))

	(defun clamp (f)
	(declare (type double-float f))
	(let ((i (* f 255.5)))
	(cond ((< i 0) 0)
	((> i 255) 255)
	(t (round i)))))

	(defun ortho-basis ()
	(let ((n (is-n isect)))
	(setf (vx b0) 0d0 (vy b0) 0d0 (vz b0) 0d0)
	(cond ((< -0.6 (vy n) 0.6) (setf (vy b0) 1d0))
	((< -0.6 (vz n) 0.6) (setf (vz b0) 1d0))
	(t (setf (vx b0) 1d0)))
	(vset b0 (multiple-value-call #'normalize (vcross b0 n)))
	(vset b1 (multiple-value-call #'normalize (vcross n b0)))
	(setq b2 n)))

	(defun ambient-occlusion (scene)
	(let* ((ntheta nao-samples)
	(nphi nao-samples)
	(eps 1d-4)
	(occlusion 0))
	(declare (type fixnum ntheta nphi occlusion))
	(multiple-value-bind (px py pz) (vadd (is-p isect) (is-n isect) eps)
	(ortho-basis)
	(dotimes% (j nphi)
	(dotimes% (i ntheta)
	(let* ((r (random 1d0))
	(phi (* 6.283185307179586d0 (random 1d0)))
	(sqr (sqrt (- 1d0 (the double-float r)))))
	(declare (type double-float r phi sqr))
	(let ((x (* (cos phi) sqr))
	(y (* (sin phi) sqr))
	(z (sqrt r)))
	(declare (type double-float x y z))
	(vsetn occ-ray
	px py pz
	(+ (* x (vx b0)) (* y (vx b1)) (* z (vx b2)))
	(+ (* x (vy b0)) (* y (vy b1)) (* z (vy b2)))
	(+ (* x (vz b0)) (* y (vz b1)) (* z (vz b2))))
	(init-isect occ-isect)
	(dolist (f scene)
	(funcall f occ-ray occ-isect))
	(when (is-hit occ-isect)
	(incf occlusion)))))))
	(/ (- (* ntheta nphi) occlusion)
	(* ntheta nphi)
	1d0)))



	(defun render (scene w h nsubs)
	(declare (type fixnum w h nsubs))
	(let ((image (make-array (list w h) :element-type '(unsigned-byte 8))))
	(dotimes% (y h)
	(dotimes% (x w)
	(let ((rad 0d0))
	(declare (type double-float rad))
	;; subsampling
	(dotimes% (v nsubs)
	(dotimes% (u nsubs)
	(let* ((px (/ (+ x (/ u (float nsubs)) (- (/ w 2d0))) (/ w 2d0)))
	(py (- (/ (+ y (/ v (float nsubs)) (- (/ h 2d0))) (/ h 2d0)))))
	(multiple-value-bind (rx ry rz) (normalize px py -1d0)
	(vsetn ray 0d0 0d0 0d0 rx ry rz))
	(init-isect isect)
	(loop for f in scene do
	(funcall f ray isect))
	(when (is-hit isect)
	(let ((col (ambient-occlusion scene)))
	(incf rad (the double-float col)))))))
	(setf (aref image x y) (clamp (/ rad (sq nsubs)))))))
	image)))


	;; write-file

	(defun write-pnm (file image w h)
	(declare (type fixnum w h))
	(with-open-file (s file :direction :output :if-exists :overwrite :if-does-not-exist :create)
	(format s "P2~%~D ~D~%~D~%" w h 255)
	(dotimes% (y h)
	(dotimes% (x w)
	(let ((p (aref image x y)))
	(format s "~D~%" p))))
	t))


	(defun main (&optional (fn "test2.pgm") (w image-width) (h image-height))
	(let ((img (render (defscene) w h nsubsamples)))
	(write-pnm fn img w h)))

	(time (main))