nobutaka/smallpt.scm

## smallpt.scm
;;;
;;; Scheme version of smallpt.
;;; see http://www.kevinbeason.com/smallpt/
;;;


(define pi 3.141592653589793238462643383279)


;;;
;;; random-real
;;;
(define random-init 1)

(define rand
  (let ((x random-init))
    (lambda ()
      (set! x (rand-update x))
      (modulo (inexact->exact (floor (/ x 65536))) 32768))))

(define (rand-update x)
  (+ (* 1103515245 x) 12345))

(define (random-real)
  (/ (rand) 32767.0))


;;;
;;; Vec
;;;
(define (make-vec x y z)
  (list x y z))

(define (x-of vec)
  (car vec))
(define (y-of vec)
  (cadr vec))
(define (z-of vec)
  (caddr vec))

(define (apply-elem op a b)
  (make-vec (op (x-of a) (x-of b))
            (op (y-of a) (y-of b))
            (op (z-of a) (z-of b))))
(define (add a b)
  (apply-elem + a b))
(define (sub a b)
  (apply-elem - a b))
(define (mult a b)
  (apply-elem * a b))

(define (scalar-mult a b)
  (make-vec (* (x-of a) b)
            (* (y-of a) b)
            (* (z-of a) b)))

(define (norm a)
  (scalar-mult a (/ (sqrt (dot a a)))))

(define (dot a b)
  (+ (* (x-of a) (x-of b))
     (* (y-of a) (y-of b))
     (* (z-of a) (z-of b))))

(define (cross a b)
  (make-vec (- (* (y-of a) (z-of b))
               (* (z-of a) (y-of b)))
            (- (* (z-of a) (x-of b))
               (* (x-of a) (z-of b)))
            (- (* (x-of a) (y-of b))
               (* (y-of a) (x-of b)))))


;;;
;;; Ray
;;;
(define (make-ray o d) (cons o d))
(define (origin r) (car r))
(define (dir r) (cdr r))


;;; material types is symbol
;;;  diff, spec, refr


;;;
;;; Sphere
;;;
(define (make-sphere rad p e c refl)
  (list rad p e c refl))

(define (rad s) (car s))
(define (position s) (cadr s))
(define (emission s) (caddr s))
(define (color s) (cadddr s))
(define (refl s) (car (cddddr s)))

(define (intersect-sphere s r)
  (let* ((op (sub (position s) (origin r)))
         (eps 1e-4)
         (b (dot op (dir r)))
         (det (+ (- (* b b)
                    (dot op op))
                 (* (rad s) (rad s)))))
    (if (< det 0)
        0
        (let* ((det1 (sqrt det))
               (t (- b det1)))
          (if (> t eps)
              t
              (let ((t1 (+ b det1)))
                (if (> t1 eps) t1 0)))))))


;;;
;;; definition of cornell box
;;;
(define spheres (vector (make-sphere 1e5  (make-vec (+ 1e5 1) 40.8 81.6)   (make-vec 0 0 0) (make-vec .75 .25 .25) 'diff)
                        (make-sphere 1e5  (make-vec (+ -1e5 99) 40.8 81.6) (make-vec 0 0 0) (make-vec .25 .25 .75) 'diff)
                        (make-sphere 1e5  (make-vec 50 40.8 1e5)           (make-vec 0 0 0) (make-vec .75 .75 .75) 'diff)
                        (make-sphere 1e5  (make-vec 50 40.8 (+ -1e5 170))  (make-vec 0 0 0) (make-vec 0 0 0) 'diff)
                        (make-sphere 1e5  (make-vec 50 1e5 81.6)           (make-vec 0 0 0) (make-vec .75 .75 .75) 'diff)
                        (make-sphere 1e5  (make-vec 50 (+ -1e5 81.6) 81.6) (make-vec 0 0 0) (make-vec .75 .75 .75) 'diff)
                        (make-sphere 16.5 (make-vec 27 16.5 47)            (make-vec 0 0 0) (scalar-mult (make-vec 1 1 1) .999) 'spec)
                        (make-sphere 16.5 (make-vec 73 16.5 78)            (make-vec 0 0 0) (scalar-mult (make-vec 1 1 1) .999) 'refr)
                        (make-sphere 600  (make-vec 50 (- 681.6 .27) 81.6) (make-vec 12 12 12) (make-vec 0 0 0) 'diff)))


;; redefinition of clamp
(define (clamp x)
  (if (< x 0)
      0
      (if (> x 1) 1 x)))

(define (toint x)
  (inexact->exact (floor (+ (* (expt (clamp x) (/ 2.2)) 255) .5))))

(define (intersect r)
  (let ((n (vector-length spheres))
        (inf 1e20))
    (define (loop i t id)
      (if (< i 0)
          (list (< t inf) t id)
          (let ((d (intersect-sphere (vector-ref spheres i) r)))
            (cond ((and (not (zero? d)) (< d t)) (loop (- i 1) d i))
                  (else (loop (- i 1) t id))))))
    (loop (- n 1) inf 0)))


(define (radiance r depth)
  (let* ((result (intersect r))
         (t (cadr result))
         (id (caddr result)))
    (if (not (car result))
        (make-vec 0 0 0)
        (let* ((obj (vector-ref spheres id))
               (x (add (origin r) (scalar-mult (dir r) t)))
               (n (norm (sub x (position obj))))
               (nl (if (< (dot n (dir r)) 0) n (scalar-mult n -1)))
               (f (color obj))
               (p (if (and (> (x-of f) (y-of f)) (> (x-of f) (z-of f))) (x-of f) (if (> (y-of f) (z-of f)) (y-of f) (z-of f))))
               (depth1 (+ depth 1))
               (return-emission #f)
               (refl-ray-dir (sub (dir r)
                                  (scalar-mult (scalar-mult n 2)
                                               (dot n (dir r))))))
          (if (> depth1 5)
              (if (< (random-real) p)
                  (set! f (scalar-mult f (/ p)))
                  (set! return-emission #t)))
          (if return-emission
              (emission obj)
              (cond
                ((eq? (refl obj) 'diff)
                 (let* ((r1 (* 2 pi (random-real)))
                        (r2 (random-real))
                        (r2s (sqrt r2))
                        (w nl)
                        (u (norm (cross (if (> (abs (x-of w)) .1) (make-vec 0 1 0) (make-vec 1 0 0)) w)))
                        (v (cross w u))
                        (d (norm (add (add (scalar-mult u (* (cos r1) r2s)) (scalar-mult v (* (sin r1) r2s))) (scalar-mult w (sqrt (- 1 r2)))))))
                   (add (emission obj) (mult f (radiance (make-ray x d) depth1)))))
                ((eq? (refl obj) 'spec)
                 (add (emission obj)
                      (mult f (radiance (make-ray x refl-ray-dir)
                                        depth1))))
                (else
                 (let* ((refl-ray (make-ray x refl-ray-dir))
                        (into (> (dot n nl) 0))
                        (nc 1)
                        (nt 1.5)
                        (nnt (if into (/ nc nt) (/ nt nc)))
                        (ddn (dot (dir r) nl))
                        (cos2t (- 1 (* nnt nnt (- 1 (* ddn ddn))))))
                   (if (< cos2t 0)
                       (add (emission obj)
                            (mult f (radiance refl-ray depth1)))
                       (let* ((tdir (norm (sub (scalar-mult (dir r)
                                                            nnt)
                                               (scalar-mult n
                                                            (* (if into 1 -1)
                                                               (+ (* ddn nnt) (sqrt cos2t)))))))
                              (a (- nt nc))
                              (b (+ nt nc))
                              (R0 (/ (* a a)
                                     (* b b)))
                              (c (- 1 (if into (- ddn) (dot tdir n))))
                              (Re (+ R0
                                     (* (- 1 R0)
                                        c
                                        c
                                        c
                                        c
                                        c)))
                              (Tr (- 1 Re))
                              (P (+ .25
                                    (* .5 Re)))
                              (RP (/ Re P))
                              (TP (/ Tr (- 1 P))))
                         (add (emission obj)
                              (mult f (if (> depth1 2)
                                          (if (< (random-real) P)
                                              (scalar-mult (radiance refl-ray depth1) RP)
                                              (scalar-mult (radiance (make-ray x tdir) depth1) TP))
                                          (add (scalar-mult (radiance (make-ray x tdir) depth1) Tr)
                                               (scalar-mult (radiance refl-ray depth1) Re)))))))))))))))


(define (main args)
  (let* ((w 32)
         (h 32)
         (samps 1)
         (cam (make-ray (make-vec 50 52 295.6) (norm (make-vec 0 -0.042612 -1))))
         (cx (make-vec (/ (* w .5135) h) 0 0))
         (cy (scalar-mult (norm (cross cx (dir cam))) .5135))
         (c (make-vector (* w h) (make-vec 0 0 0))))
    (define (sample x y sx sy)
      (define (make-sample-ray)
        (let* ((r1 (* 2 (random-real)))
               (dx (if (< r1 1) (- (sqrt r1) 1) (- 1 (sqrt (- 2 r1)))))
               (r2 (* 2 (random-real)))
               (dy (if (< r2 1) (- (sqrt r2) 1) (- 1 (sqrt (- 2 r2)))))
               (d (norm (add (add (scalar-mult cx (- (/ (+ (/ (+ sx .5 dx) 2) x) w) .5))
                                  (scalar-mult cy (- (/ (+ (/ (+ sy .5 dy) 2) y) h) .5)))
                             (dir cam)))))
          (make-ray (add (origin cam) (scalar-mult d 140)) d)))
      (define (accmulate-radiance s r)
        (if (= s 0)
            r
            (accmulate-radiance (- s 1) (add r (scalar-mult (radiance (make-sample-ray) 0) (/ samps))))))
      (accmulate-radiance samps (make-vec 0 0 0)))
    (let loop-rows ((y 0))
      (when (< y h)
            (display (format "\rRendering (~d spp) ~d%" (* samps 4) (/ (* 100.0 y) (- h 1))) (standard-error-port))
            (let loop-cols ((x 0))
              (when (< x w)
                    (let ((i (+ (* (- h y 1) w) x)))
                      (let loop-subpixel-rows ((sy 0))
                        (when (< sy 2)
                              (let loop-subpixel-cols ((sx 0))
                                (when (< sx 2)
                                      (let ((r (sample x y sx sy)))
                                        (vector-set! c i (add (vector-ref c i) (scalar-mult (make-vec (clamp (x-of r)) (clamp (y-of r)) (clamp (z-of r))) .25))))
                                      (loop-subpixel-cols (+ sx 1))))
                              (loop-subpixel-rows (+ sy 1)))))
                    (loop-cols (+ x 1))))
            (loop-rows (+ y 1))))
    (display (format "P3\n~d ~d\n~d\n" w h 255))
    (let loop-pixels ((i 0))
      (when (< i (* w h))
            (display (format "~d ~d ~d " (toint (x-of (vector-ref c i))) (toint (y-of (vector-ref c i))) (toint (z-of (vector-ref c i)))))
            (loop-pixels (+ i 1))))))
	;;;
	;;; Scheme version of smallpt.
	;;; see http://www.kevinbeason.com/smallpt/
	;;;


	(define pi 3.141592653589793238462643383279)


	;;;
	;;; random-real
	;;;
	(define random-init 1)

	(define rand
	(let ((x random-init))
	(lambda ()
	(set! x (rand-update x))
	(modulo (inexact->exact (floor (/ x 65536))) 32768))))

	(define (rand-update x)
	(+ (* 1103515245 x) 12345))

	(define (random-real)
	(/ (rand) 32767.0))


	;;;
	;;; Vec
	;;;
	(define (make-vec x y z)
	(list x y z))

	(define (x-of vec)
	(car vec))
	(define (y-of vec)
	(cadr vec))
	(define (z-of vec)
	(caddr vec))

	(define (apply-elem op a b)
	(make-vec (op (x-of a) (x-of b))
	(op (y-of a) (y-of b))
	(op (z-of a) (z-of b))))
	(define (add a b)
	(apply-elem + a b))
	(define (sub a b)
	(apply-elem - a b))
	(define (mult a b)
	(apply-elem * a b))

	(define (scalar-mult a b)
	(make-vec (* (x-of a) b)
	(* (y-of a) b)
	(* (z-of a) b)))

	(define (norm a)
	(scalar-mult a (/ (sqrt (dot a a)))))

	(define (dot a b)
	(+ (* (x-of a) (x-of b))
	(* (y-of a) (y-of b))
	(* (z-of a) (z-of b))))

	(define (cross a b)
	(make-vec (- (* (y-of a) (z-of b))
	(* (z-of a) (y-of b)))
	(- (* (z-of a) (x-of b))
	(* (x-of a) (z-of b)))
	(- (* (x-of a) (y-of b))
	(* (y-of a) (x-of b)))))


	;;;
	;;; Ray
	;;;
	(define (make-ray o d) (cons o d))
	(define (origin r) (car r))
	(define (dir r) (cdr r))


	;;; material types is symbol
	;;; diff, spec, refr


	;;;
	;;; Sphere
	;;;
	(define (make-sphere rad p e c refl)
	(list rad p e c refl))

	(define (rad s) (car s))
	(define (position s) (cadr s))
	(define (emission s) (caddr s))
	(define (color s) (cadddr s))
	(define (refl s) (car (cddddr s)))

	(define (intersect-sphere s r)
	(let* ((op (sub (position s) (origin r)))
	(eps 1e-4)
	(b (dot op (dir r)))
	(det (+ (- (* b b)
	(dot op op))
	(* (rad s) (rad s)))))
	(if (< det 0)
	0
	(let* ((det1 (sqrt det))
	(t (- b det1)))
	(if (> t eps)
	t
	(let ((t1 (+ b det1)))
	(if (> t1 eps) t1 0)))))))


	;;;
	;;; definition of cornell box
	;;;
	(define spheres (vector (make-sphere 1e5 (make-vec (+ 1e5 1) 40.8 81.6) (make-vec 0 0 0) (make-vec .75 .25 .25) 'diff)
	(make-sphere 1e5 (make-vec (+ -1e5 99) 40.8 81.6) (make-vec 0 0 0) (make-vec .25 .25 .75) 'diff)
	(make-sphere 1e5 (make-vec 50 40.8 1e5) (make-vec 0 0 0) (make-vec .75 .75 .75) 'diff)
	(make-sphere 1e5 (make-vec 50 40.8 (+ -1e5 170)) (make-vec 0 0 0) (make-vec 0 0 0) 'diff)
	(make-sphere 1e5 (make-vec 50 1e5 81.6) (make-vec 0 0 0) (make-vec .75 .75 .75) 'diff)
	(make-sphere 1e5 (make-vec 50 (+ -1e5 81.6) 81.6) (make-vec 0 0 0) (make-vec .75 .75 .75) 'diff)
	(make-sphere 16.5 (make-vec 27 16.5 47) (make-vec 0 0 0) (scalar-mult (make-vec 1 1 1) .999) 'spec)
	(make-sphere 16.5 (make-vec 73 16.5 78) (make-vec 0 0 0) (scalar-mult (make-vec 1 1 1) .999) 'refr)
	(make-sphere 600 (make-vec 50 (- 681.6 .27) 81.6) (make-vec 12 12 12) (make-vec 0 0 0) 'diff)))


	;; redefinition of clamp
	(define (clamp x)
	(if (< x 0)
	0
	(if (> x 1) 1 x)))

	(define (toint x)
	(inexact->exact (floor (+ (* (expt (clamp x) (/ 2.2)) 255) .5))))

	(define (intersect r)
	(let ((n (vector-length spheres))
	(inf 1e20))
	(define (loop i t id)
	(if (< i 0)
	(list (< t inf) t id)
	(let ((d (intersect-sphere (vector-ref spheres i) r)))
	(cond ((and (not (zero? d)) (< d t)) (loop (- i 1) d i))
	(else (loop (- i 1) t id))))))
	(loop (- n 1) inf 0)))


	(define (radiance r depth)
	(let* ((result (intersect r))
	(t (cadr result))
	(id (caddr result)))
	(if (not (car result))
	(make-vec 0 0 0)
	(let* ((obj (vector-ref spheres id))
	(x (add (origin r) (scalar-mult (dir r) t)))
	(n (norm (sub x (position obj))))
	(nl (if (< (dot n (dir r)) 0) n (scalar-mult n -1)))
	(f (color obj))
	(p (if (and (> (x-of f) (y-of f)) (> (x-of f) (z-of f))) (x-of f) (if (> (y-of f) (z-of f)) (y-of f) (z-of f))))
	(depth1 (+ depth 1))
	(return-emission #f)
	(refl-ray-dir (sub (dir r)
	(scalar-mult (scalar-mult n 2)
	(dot n (dir r))))))
	(if (> depth1 5)
	(if (< (random-real) p)
	(set! f (scalar-mult f (/ p)))
	(set! return-emission #t)))
	(if return-emission
	(emission obj)
	(cond
	((eq? (refl obj) 'diff)
	(let* ((r1 (* 2 pi (random-real)))
	(r2 (random-real))
	(r2s (sqrt r2))
	(w nl)
	(u (norm (cross (if (> (abs (x-of w)) .1) (make-vec 0 1 0) (make-vec 1 0 0)) w)))
	(v (cross w u))
	(d (norm (add (add (scalar-mult u (* (cos r1) r2s)) (scalar-mult v (* (sin r1) r2s))) (scalar-mult w (sqrt (- 1 r2)))))))
	(add (emission obj) (mult f (radiance (make-ray x d) depth1)))))
	((eq? (refl obj) 'spec)
	(add (emission obj)
	(mult f (radiance (make-ray x refl-ray-dir)
	depth1))))
	(else
	(let* ((refl-ray (make-ray x refl-ray-dir))
	(into (> (dot n nl) 0))
	(nc 1)
	(nt 1.5)
	(nnt (if into (/ nc nt) (/ nt nc)))
	(ddn (dot (dir r) nl))
	(cos2t (- 1 (* nnt nnt (- 1 (* ddn ddn))))))
	(if (< cos2t 0)
	(add (emission obj)
	(mult f (radiance refl-ray depth1)))
	(let* ((tdir (norm (sub (scalar-mult (dir r)
	nnt)
	(scalar-mult n
	(* (if into 1 -1)
	(+ (* ddn nnt) (sqrt cos2t)))))))
	(a (- nt nc))
	(b (+ nt nc))
	(R0 (/ (* a a)
	(* b b)))
	(c (- 1 (if into (- ddn) (dot tdir n))))
	(Re (+ R0
	(* (- 1 R0)
	c
	c
	c
	c
	c)))
	(Tr (- 1 Re))
	(P (+ .25
	(* .5 Re)))
	(RP (/ Re P))
	(TP (/ Tr (- 1 P))))
	(add (emission obj)
	(mult f (if (> depth1 2)
	(if (< (random-real) P)
	(scalar-mult (radiance refl-ray depth1) RP)
	(scalar-mult (radiance (make-ray x tdir) depth1) TP))
	(add (scalar-mult (radiance (make-ray x tdir) depth1) Tr)
	(scalar-mult (radiance refl-ray depth1) Re)))))))))))))))


	(define (main args)
	(let* ((w 32)
	(h 32)
	(samps 1)
	(cam (make-ray (make-vec 50 52 295.6) (norm (make-vec 0 -0.042612 -1))))
	(cx (make-vec (/ (* w .5135) h) 0 0))
	(cy (scalar-mult (norm (cross cx (dir cam))) .5135))
	(c (make-vector (* w h) (make-vec 0 0 0))))
	(define (sample x y sx sy)
	(define (make-sample-ray)
	(let* ((r1 (* 2 (random-real)))
	(dx (if (< r1 1) (- (sqrt r1) 1) (- 1 (sqrt (- 2 r1)))))
	(r2 (* 2 (random-real)))
	(dy (if (< r2 1) (- (sqrt r2) 1) (- 1 (sqrt (- 2 r2)))))
	(d (norm (add (add (scalar-mult cx (- (/ (+ (/ (+ sx .5 dx) 2) x) w) .5))
	(scalar-mult cy (- (/ (+ (/ (+ sy .5 dy) 2) y) h) .5)))
	(dir cam)))))
	(make-ray (add (origin cam) (scalar-mult d 140)) d)))
	(define (accmulate-radiance s r)
	(if (= s 0)
	r
	(accmulate-radiance (- s 1) (add r (scalar-mult (radiance (make-sample-ray) 0) (/ samps))))))
	(accmulate-radiance samps (make-vec 0 0 0)))
	(let loop-rows ((y 0))
	(when (< y h)
	(display (format "\rRendering (~d spp) ~d%" (* samps 4) (/ (* 100.0 y) (- h 1))) (standard-error-port))
	(let loop-cols ((x 0))
	(when (< x w)
	(let ((i (+ (* (- h y 1) w) x)))
	(let loop-subpixel-rows ((sy 0))
	(when (< sy 2)
	(let loop-subpixel-cols ((sx 0))
	(when (< sx 2)
	(let ((r (sample x y sx sy)))
	(vector-set! c i (add (vector-ref c i) (scalar-mult (make-vec (clamp (x-of r)) (clamp (y-of r)) (clamp (z-of r))) .25))))
	(loop-subpixel-cols (+ sx 1))))
	(loop-subpixel-rows (+ sy 1)))))
	(loop-cols (+ x 1))))
	(loop-rows (+ y 1))))
	(display (format "P3\n~d ~d\n~d\n" w h 255))
	(let loop-pixels ((i 0))
	(when (< i (* w h))
	(display (format "~d ~d ~d " (toint (x-of (vector-ref c i))) (toint (y-of (vector-ref c i))) (toint (z-of (vector-ref c i)))))
	(loop-pixels (+ i 1))))))