;;; Scheme Recursive Art Contest Entry
;;;
;;; Please do not include your name or personal info in this file.
;;;
;;; Title: Porifera
;;;
;;; Description:
;;;    The bouncing light rays.
;;;    Travel through the Menger sponge.
;;;    For your eyes to see.

;; NOTE: uncomment the block below in order to run in Racket,
;;       which is much faster than the 61A interpreter

;;;; #lang scheme
;;;; (define racket false) ;;change this to true for racket
;;;; (require graphics/graphics)
;;;; (define (speed x) '())
;;;; (define (setposition x y) '())
;;;; (define drawcolor "black")
;;;; (define (color c) (set! drawcolor c))
;;;; (define (penup) '())
;;;; (define (pendown) '())
;;;; (define (pu) '())
;;;; (define (pd) '())
;;;; (define (begin_fill) '())
;;;; (define (end_fill) '())
;;;; (define (rgb r g b) (make-rgb r g b))
;;;; (define (exitonclick) '())
;;;; (define (update) '())
;;;; (define (setpos x y) '())
;;;; (define w '())
;;;; (define (bgcolor c)
;;;;   (begin
;;;;     (open-graphics)
;;;;     ;; nothing appears to happen, but the library is initialized...
;;;; 
;;;;     (set! w (open-viewport "practice" draw-size draw-size))
;;;; 
;;;;     ((draw-viewport w) "black")))
;;;; (define (hideturtle) '())
;;;; (define (tracer a b) '())
;;;; (define nil '())
;;;; (define-syntax if
;;;;   (syntax-rules ()
;;;;     [(if c a) (if c a nil)]
;;;;     [(if c a b) (cond (c a) (else b))]))
;;;; (define (drawpixel x y)
;;;;   ((draw-pixel w) (make-posn (+ (/ draw-size 2) x) (-(/ draw-size 2) y)) drawcolor))

;; Comment this out if using racket
(define (drawpixel x y)
  (pu)
  (setpos (- (* pixelscale x) (/ pixelscale 2)) (- (* pixelscale y) (/ pixelscale 2)))
  (begin_fill)
  (setpos (- (* pixelscale x) (/ pixelscale 2)) (+ (* pixelscale y) (/ pixelscale 2)))
  (setpos (+ (* pixelscale x) (/ pixelscale 2)) (+ (* pixelscale y) (/ pixelscale 2)))
  (setpos (+ (* pixelscale x) (/ pixelscale 2)) (- (* pixelscale y) (/ pixelscale 2)))
  (end_fill))

;;END Racket block

;; color parameters
(define sponge-reflectivity '(0.9 0.9 0.9))
(define ground-reflectivity '(0.8 0.8 0.8))
(define center-light-color '(10 0 0))
(define sky-color '(0.05 0.05 0.05))
(define ambient-color (list 0.01 0.01 0.01))
;;

(define pi 3.14159265358979323846264338327950)

(define (cddr x) (cdr (cdr x)))
(define (cdddr x) (cddr (cdr x)))
(define (cddddr x) (cdddr (cdr x)))
(define (cdddddr x) (cddddr (cdr x)))
(define (cddddddr x) (cdddddr (cdr x)))

(define (cadr x) (car (cdr x)))
(define (caddr x) (car (cddr x)))
(define (cadddr x) (car (cdddr x)))
(define (caddddr x) (car (cddddr x)))
(define (cadddddr x) (car (cdddddr x)))
(define (caddddddr x) (car (cddddddr x)))

(define x car)
(define y cadr)
(define z caddr)

(define (clamp-rgb c) (map (lambda (x) (clamp x 0 1)) c))

(define (round x)
  (floor (+ x 0.5)))

(define (sign x) (if (> x 0)
                     +
                     -))

(define (safe/ x y) (if (= y 0) 0 (/ x y)))

(define (min x h)
  (if (> x h)
      h
      x))
(define (max x l)
  (if (< x l)
      l
      x))
(define (clamp x l h) (min (max x l) h))

(define (apply-to-list func list)
  (if (null? list)
      nil
      (begin (func (car list))
             (apply-to-list func (cdr list)))))

(define (apply-stack s x) (if (null? s) x (apply-stack (cdr s) ((car s) x))))

(define (map func list)
  (if (null? list)
      nil
      (cons (func (car list))
            (map func (cdr list)))))

(define (list-ref list i)
  (if (= i 0)
      (car list)
      (list-ref (cdr list) (- i 1))))

(define (list-set list i val)
  (if (= i 0)
      (cons val (cdr list))
      (cons (car list)
            (list-set (cdr list) (- i 1) val))))

(define (swizzle v is)
  (map (lambda (i) (list-ref v (abs i))) is))

(define (scale v s)
  (map
   (lambda (x) (* x s))
   v))

(define (component-wise func v1 v2)
  (if (null? v1)
      nil
      (cons (func (car v1) (car v2))
            (component-wise func (cdr v1) (cdr v2)))))

(define (add v1 v2)
  (component-wise + v1 v2))


(define (sub v1 v2)
  (component-wise - v1 v2))

(define (apply-to-element func e lst)
  (cond ((null? (car lst)) nil)
        ((= e 0) (cons (func (car lst)) (cdr lst)))
        (else (cons (car lst) (apply-to-element func (- e 1) (cdr lst))))))

(define (dot-helper v1 v2 prod)
  (if (null? v1)
      prod
      (dot-helper (cdr v1) (cdr v2)
                  (+ (* (car v1)
                        (car v2))
                     prod))))

(define (dot v1 v2)
  (dot-helper v1 v2 0))

(define (norm v)
  (sqrt (dot v v)))

(define (normalize v)
  (scale v (/ 1.0 (norm v))))

(define (projection a b)
  (scale b (/ (dot a b) (dot b b))))

(define (rejection a b)
  (sub a (projection a b)))

(define (reflection a b)
  (sub a (scale (projection a b) 2)))

(define (count-negatives v)
  (if (null? v)
      0
      (if (< (car v) 0)
          (+ 1 (count-negatives (cdr v)))
          (count-negatives (cdr v)))))

;;; returns a transformation which rotates a point about the y axis
(define (roty v theta)
  (list (- (* (x v) (cos theta)) (* (z v) (sin theta)))
        (y v)
        (+ (* (x v) (sin theta)) (* (z v) (cos theta)))))

;;; returns a transformation which rotates a point about the x axis
(define (rotx v theta)
  (list (x v)
        (- (* (y v) (cos theta)) (* (z v) (sin theta)))
        (+ (* (y v) (sin theta)) (* (z v) (cos theta)))))

(define draw-size 512)
(define pixelscale 1)

(define scale-x (* (sqrt 2) (/ 3 draw-size)))
(define scale-y (* (sqrt 2) (/ 3 draw-size)))

(define epsilon 0.0001)

(define (center-epsilon? i)
  (and (< (abs (x i)) (+ 0.5 epsilon))
       (< (abs (y i)) (+ 0.5 epsilon))))

(define (=epsilon x y)
  (< (abs (- x y)) 0.001))

;; casts a ray inside a cube, starting from the -z side,
;; returns the integer coordinates of side which the ray hits
(define (cubecast p0 d)
  (let ((i0 (map round p0))
        (tz (/ 1.0 (abs (z d)))))
    (if (center-epsilon? p0)
        (begin
          (define p1 (add p0 (scale d tz)))
          (define i1 (map round p1))
          (if (center-epsilon? p1)
              (list (list 0 0 ((sign (z d)) 1)) tz) ;; the ray hits on the z side
              ;; check whether the ray will hit in x or y first
              (let ((ti (lambda (x) (safe/ (- ((sign (x i1)) 0.5) (x p0)) (x d)))))
                (let ((tx (ti x))
                      (ty (ti y)))
                  (define side (cond
                                ((=epsilon (x i1) 0) false)
                                ((=epsilon (y i1) 0) true)
                                ((< tx 0) false)
                                ((< ty 0) true)
                                (else (< tx ty))))
                  (if side
                      (list (list ((sign (x i1)) 1) 0 0) tx)
                      (list (list 0 ((sign (y i1)) 1) 0) ty))))))
        ;; the ray is starting off the cube
        (list (append i0 '(0)) 0))))

;; parameters are the same as spongecast
(define (rotatedspongecast yc ty p i is d t r c trans level reflections)
  (let ((py (list-ref p yc))
        (dy (list-ref d yc))
        (swizzle-factor (cond ((= yc 1) '(0 2 1))
                              ((= yc 0) '(2 1 0)))))
    (define (forward-transformation v) (swizzle v swizzle-factor))
    (cond
     ;; the ray is going away from the subcube
     ((>= ((sign dy) py) 1.5)
      (if (null? is) ;; we are at the top level
          (list t r c reflections
                ((car trans) p)
                ((car trans) (forward-transformation
                              (list 0 0 ((sign dy) 1))))
                ((car trans) d))
          (spongecast (add (scale ((car trans) p) (/ 1 3)) i)
                      (add ((car trans)
                            (forward-transformation (list 0 0 ((sign dy) 1))))
                           (car is))
                      (cdr is)
                      ((car trans) d)
                      (/ t 3)
                      r c
                      (cdr trans)
                      (+ 1 level)
                      reflections)))
     (else
      (define newp (forward-transformation
                    (add p (scale d ty))))
      (spongecast newp
                  (map round (add newp (list 0 0 ((sign dy) 0.5))))
                  is
                  (forward-transformation d)
                  (+ t ty)
                  r c
                  (cons (lambda (v) ((car trans) (forward-transformation v)))
                        (cdr trans))
                  level
                  reflections)))))

(define max-reflections 10)
(define max-level 4)

(define (point-light p d c light_radius)
  (let ((radius (norm (rejection p d))))
    (if (> radius light_radius)
        '(0 0 0)
        (scale c (* 0.5 (- light_radius radius))))))

;; returns the information needed to calculate the output color if there was a hit
(define (spongecast p ; position in the current level
                    i ; grid position in the current level
                    is ; stack of grid positions from upper levels
                    d ; direction of the light ray
                    t ; distance the ray has traveled divided by the norm of d
                    r ; total reflectance between the current ray and the eye
                    oc ; total color
                    trans ; stack of symmetry tranforms used
                    level ; the level of the menger sponge
                    reflections) ; the number of times the ray has been reflected
  ;; place a light source in the center
  (define c
    (if (and (null? is) (< (dot i i) 0.1))
        (add oc (component-wise
                 * r
                 (point-light (sub p i) d center-light-color 0.5)))
        oc))
  
  (cond
   ;; the ray is passing out the front or back
   ((> (abs (z i)) 1)
    (if (null? is) ;; we are at the top level, return
        (list t r c reflections
              ((car trans) p)
              ((car trans) (list 0 0 ((sign (z d)) 1)))
              ((car trans) d))
        ;; otherwise step up one level and continue
        (spongecast (add (scale ((car trans) p) (/ 1 3)) (car is))
                    (add ((car trans) (list 0 0 ((sign (z i)) 1)))
                         (car is))
                    (cdr is)
                    ((car trans) d)
                    (/ t 3)
                    r c
                    (cdr trans)
                    (+ 1 level)
                    reflections)))
   
   ;; the ray is off of the z face and hitting one of the sides
   ((or (> (abs (x i)) 1) (> (abs (y i)) 1))
    (define (ti x) (safe/ (- ((sign (x d)) -1.5) (x p)) (x d)))
    (define tx (ti x))
    (define ty (ti y))
    
    (define side (cond
                  ((<= (abs (x i)) 1) false)
                  ((<= (abs (y i)) 1) true)
                  (else (> tx ty))))
    
    (if side
        (rotatedspongecast 0 tx p i is d t r c trans level reflections)
        (rotatedspongecast 1 ty p i is d t r c trans level reflections)))
   
   ;; go forward
   (else
    (if (> (dot i i) 1) ;; this checks whether the subcube is filled at this level
        ;; subcube is filled
        (let ((p0 (scale (sub p i) 3))
              (reflection-normal (list 0 0 ((sign (z d)) -1))))
          (if (= level 1)
              (let ((nd (reflection d reflection-normal))
                    (nr (component-wise * r sponge-reflectivity))
                    (nc (add c (component-wise * r ambient-color))))
                (if (< reflections max-reflections)
                    ;; reflect if we are under the limit
                    (spongecast p
                                (add i reflection-normal)
                                is
                                nd
                                t
                                nr
                                nc
                                trans
                                level
                                (+ reflections 1))
                    ;; otherwise return what we have
                    (list t nr nc reflections
                          nil
                          '(0 0 0)
                          nil)))
              ;; go into the smaller subcube
              (spongecast p0
                          (list-set (map round (add p0 (list 0 0 ((sign (z d)) 0.5))))
                                    2 ((sign (z d)) -1))
                          (cons i is)
                          d
                          (* 3 t)
                          r c
                          (cons (lambda (v) v) trans)
                          (- level 1)
                          reflections)))
        ;; subcube is not filled
        (let ((di (cubecast (swizzle (sub p i) '(0 1)) d)))
          (let ((newp (add p (scale d (cadr di))))
                (newi (add i (car di)))
                (newt (+ t (cadr di))))
            (cond
             ;; the ray is going off one of the sides
             ((and (=epsilon 1 (abs (x (car di)))) (=epsilon 0 (y (car di))))
              (rotatedspongecast 0 0 newp newi is d newt r c trans level reflections))
             ((and (=epsilon 1 (abs (y (car di)))) (=epsilon 0 (x (car di))))
              (rotatedspongecast 1 0 newp newi is d newt r c trans level reflections))
             
             ;; else we propogate normally
             (else
              (spongecast newp newi is d newt r c trans level reflections)))))))))

(define (random i)
  (sin (* 13219839 i)))

(define (randvec i) ;; not completely uniform, but close enough
  (normalize (list (random i) (random (+ 1 i)) (random (+ 2 i)))))

(define (avg-intensity p d t r c level reflections)
  (define dcube (normalize (scale p -1))) ;; look towards the cube
  
  (define (cast-ray s d c samples mix spread)
    (if (>= s samples)
        c
        (let ((nd (add d (scale (randvec (+ s (norm c))) spread))))
          (cast-ray (+ s 1)
                    nd
                    (add c (scale (scene p nd t r c level reflections)
                                  (/ mix samples)))
                    samples mix spread))))
  (clamp-rgb (cast-ray
              0
              '(0 1 0)
              (cast-ray
               0 dcube ;; most of the interesting light is from the dcube direction
               (cast-ray
                0
                (reflection d '(0 1 0))
                c 1 0.25 0)
               20 0.5 0.075)
              10 0.25 0.1)))

(define (scene p d t r c level reflections)
  (let ((tz (safe/ (- ((sign (z d)) -1.5) (z p)) (z d))))
    (let ((p0 (add p (scale d tz))))
      (define subcube (spongecast p0
                                  (list-set (map round p0) 2 -1)
                                  nil
                                  d
                                  (+ t tz)
                                  r c
                                  (list (lambda (v) v))
                                  level
                                  reflections))
      (let ((nt (car subcube))
            (nr (cadr subcube))
            (nc (caddr subcube))
            (nreflections (cadddr subcube))
            (np (caddddr subcube))
            (di (cadddddr subcube))
            (nd (caddddddr subcube)))
        (if (and (> (norm nr) 0.01) ;; stop if too little light is reaching the camera
                 (< nreflections max-reflections) (not (equal? di '(0 0 0))))
            (let ((ty (safe/ (- -1.5 (y np)) (y nd))))
              (define p-yplane (add np (scale nd ty)))
              (define nnc (add nc (component-wise * nr
                                                  (scale sky-color
                                                         (max 0
                                                              (dot (normalize nd)
                                                                   '(0 2 0)))))))
              (cond ((< (y nd) 0)
                     (avg-intensity p-yplane
                                    d
                                    (+ nt ty)
                                    (component-wise * nr ground-reflectivity)
                                    nc
                                    level (+ 7 nreflections)))
                    (else (clamp-rgb nnc))))
            (clamp-rgb nc))))))

(define (singlesample x y)
  (scene (roty (rotx (list (* x scale-x) (* y scale-y) -5.0)
                     (/ pi 12))
               (/ pi 8))
         (roty (rotx (list (* x scale-x) (* y scale-y) 5.0)
                     (/ pi 12))
               (/ pi 8))
         0
         (list 1 1 1)
         (list 0 0 0)
         max-level 0))

(define (multisample x y)
  (scale
   (add
    (add (singlesample (+ x 0.333) (+ y 0.333))
         (singlesample (+ x 0.333) (+ y 0.333)))
    (add
     (singlesample (+ x 0.333) (+ y 0.333))
     (singlesample (+ x 0.333) (+ y 0.333))))
   0.25))

(define (render x y)
  (color (apply rgb (multisample x y)))
  (drawpixel x y)
  (update))

(define (fill-recursively size x y)
  (if (= size 1)
      (render x y)
      (let ((newsize (quotient size 2))
            (offset (/ size 4)))
        (fill-recursively newsize (- x offset) (- y offset))
        (fill-recursively newsize (+ x offset) (+ y offset))
        (fill-recursively newsize (+ x offset) (- y offset))
        (fill-recursively newsize (- x offset) (+ y offset)))))

(define (draw)
  (bgcolor "#000000")
  (color (rgb 1 1 1))
  (speed 0)
  (hideturtle)
  (tracer 0 0)
  (fill-recursively draw-size 0 0)
  (exitonclick))

; Please leave this last line alone.  You may add additional procedures above
; this line.
(draw)