aobench CommonLisp ver.

aobench.lisp

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

;; vector
(defmacro vx (v) `(svref ,v 0))
(defmacro vy (v) `(svref ,v 1))
(defmacro vz (v) `(svref ,v 2))

(defun vadd (a b)
  (vector (+ (vx a) (vx b))
          (+ (vy a) (vy b))
          (+ (vz a) (vz b))))

(defun vsub (a b)
  (vector (- (vx a) (vx b))
          (- (vy a) (vy b))
          (- (vz a) (vz b))))

(defun vcross (a b)
  (vector (- (* (vy a) (vz b)) (* (vz a) (vy b)))
          (- (* (vz a) (vx b)) (* (vx a) (vz b)))
          (- (* (vx a) (vy b)) (* (vy a) (vx b)))))

(defun vdot (a b)
  (+ (* (vx a) (vx b))
     (* (vy a) (vy b))
     (* (vz a) (vz b))))

(defun sq (x) (* x x))

(defun vlen (a)
  (sqrt (+ (sq (vx a)) (sq (vy a)) (sq (vz a)))))

(defun vnormalize (a)
  (let ((d (vlen a)))
    (if (> d 1.d-17)
        (vector (/ (vx a) d) (/ (vy a) d) (/ (vz a) d))
        a)))


;; geometry
(defstruct ray
  org dir)

(defstruct (isect (:conc-name is-))
  (dist 1.0d30)
  (hit nil)
  (p (vector 0.0 0.0 0.0))
  (n (vector 0.0 0.0 0.0)))

(defun sphere (center radius)
  #'(lambda (ray isect)
      (let* ((rs (vsub (ray-org ray) center))
             (b (vdot rs (ray-dir ray)))
             (c (- (vdot rs rs) (* radius radius)))
             (d (- (sq b) c)))
        (when (> d 0.0)
          (let ((dist (- (- b) (sqrt d))))
            (when (< 0.0 dist (is-dist isect))
              (setf (is-dist isect) dist
                    (is-hit isect) t
                    (is-p isect)
                    (let ((ro (ray-org ray)) (rd (ray-dir ray)))
                      (vector (+ (vx ro) (* (vx rd) dist))
                              (+ (vy ro) (* (vy rd) dist))
                              (+ (vz ro) (* (vz rd) dist))))
                    (is-n isect) (vnormalize (vsub (is-p isect) center)))))))))

(defun plane (p n)
  #'(lambda (ray isect)
      (let ((d (-(vdot p n)))
            (v (vdot (ray-dir ray) n)))
        (when (> (abs v) 1.d-17)
          (let ((dist (/ (- (+ (vdot (ray-org ray) n) d)) v)))
            (when (< 0.0 dist (is-dist isect))
		(setf (is-dist isect) dist
		      (is-hit isect) t
		      (is-n isect) n
		      (is-p isect)
		      (let ((ro (ray-org ray)) (rd (ray-dir ray)))
			(vector (+ (vx ro) (* (vx rd) dist))
				(+ (vy ro) (* (vy rd) dist))
				(+ (vz ro) (* (vz rd) dist)))))))))))

(defun defscene ()
  (list (sphere (vector -2.0 0.0 -3.5) 0.5)
        (sphere (vector -0.5 0.0 -3.0) 0.5)
        (sphere (vector 1.0 0.0 -2.2) 0.5)
        (plane (vector 0.0 -0.5 0.0) (vector 0.0 1.0 0.0))))

(defun ortho-basis (n)
  (let* ((v (cond
	      ((< -0.6 (vy n) 0.6)
	       (vector 0.0 1.0 0.0))
	      ((< -0.6 (vz n) 0.6)
	       (vector 0.0 0.0 1.0))
	      (t (vector 1.0 0.0 0.0))))
         (s (vnormalize (vcross v n))))
    (values s (vnormalize (vcross n s)) n)))

(defun random-real ()
  (/ (random 1d16) 1d16))

(defun ambient-occlusion (scene isect)
  (let* ((ntheta nao-samples)
         (nphi nao-samples)
         (eps 0.0001)
         (occlusion 0.0)
         (p (vector (+ (vx (is-p isect)) (* eps (vx (is-n isect))))
                    (+ (vy (is-p isect)) (* eps (vy (is-n isect))))
                    (+ (vz (is-p isect)) (* eps (vz (is-n isect)))))))
    (multiple-value-bind (b0 b1 b2) (ortho-basis (is-n isect))
      (dotimes (j nphi)
        (dotimes (i ntheta)
          (let* ((r (random-real)) (phi (* 2.0 pi (random-real)))
                 (x (* (cos phi) (sqrt (- 1 r))))
                 (y (* (sin phi) (sqrt (- 1 r))))
                 (z (sqrt r))
                 (newdir (vector (+ (* 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)))))
                 (newray (make-ray :org p :dir newdir))
                 (occ-isect (make-isect)))
            (mapc #'(lambda (f) (funcall f newray occ-isect)) scene)
            (when (is-hit occ-isect)
              (incf occlusion 1.0))))))
    (/ (- (* ntheta nphi) occlusion) (* ntheta nphi))))


(defun clamp (f)
  (let ((i (* f 255.5)))
    (cond ((< i 0) 0)
          ((> i 255) 255)
          (t (round i)))))

(defun render (scene w h nsubs)
  (let ((image (make-array (list w h) :element-type 'vector)))
    (dotimes (y h image)
      (dotimes (x w)
        (let ((rad 0.0))
          ;; subsampling
          (dotimes (v nsubs)
            (dotimes (u nsubs)
              (let* ((px (/ (+ x (/ u nsubs) (- (/ w 2))) (/ w 2)))
                     (py (- (/ (+ y (/ v nsubs) (- (/ h 2))) (/ h 2))))
                     (eye (vnormalize (vector px py -1.0)))
                     (newray (make-ray :org (vector 0.0 0.0 0.0) :dir eye))
                     (isect (make-isect)))
                (mapc #'(lambda (f) (funcall f newray isect)) scene)
                (when (is-hit isect)
                  (let ((col (ambient-occlusion scene isect)))
                    (incf rad col))))))
          (setf (aref image x y) (clamp (/ rad (sq nsubs)))))))))


(defun write-pnm (file image 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 "test.pgm"))
  (let ((img (render (defscene) image-width image-height nsubsamples)))
    (write-pnm fn img image-width image-height)))