utah-teapot.rkt

#lang racket/base
(require racket/match racket/list racket/port racket/string 
         flomat)

;;;
;;; Homogeneous coordinates for 3 dimensions
;;;

; Follow the tutorial
;   http://www.opengl-tutorial.org/beginners-tutorials/tutorial-3-matrices/

(define (point x y z) (column x y z 1))
(define (dir   x y z) (column x y z 0))
(define (normalize p) (column (xc p) (yc p) (zc p)))

(define (xc p)
  (define x (ref p 0 0))
  (define w (ref p 3 0))
  (if (= w 1.0) x (/ x w)))

(define (yc p)
  (define y (ref p 1 0))
  (define w (ref p 3 0))
  (if (= w 1.0) y (/ y w)))

(define (zc p)
  (define z (ref p 2 0))
  (define w (ref p 3 0))
  (if (= w 1.0) z (/ z w)))

;;;
;;; Translation
;;;

(define (translation/xyz x y z)
  (flomat: [[1 0 0 x]
            [0 1 0 y]
            [0 0 1 z]
            [0 0 0 1]]))

(define (translation . args)
  (match args
    [(list p)     (translation/xyz (xc p) (yc p) (zc p))]
    [(list x y z) (translation/xyz x y z)]
    [_ (error 'translation)]))

;;;
;;; Scaling
;;;

(define (scaling/xyz sx sy sz)
  (flomat: [[sx 0  0  0]
            [0  sy 0  0]
            [0  0  sz 0]
            [0  0  0  1]]))

(define (scaling . args)
  (match args
    [(list s)        (scaling/xyz s  s  s)]
    [(list sx sy sz) (scaling/xyz sx sy sz)]
    [_ (error 'scaling)]))

;;;
;;; Rotation
;;;

(define (rotation-z γ)
  (define cγ (cos γ)) (define sγ (sin γ))
  (flomat: 
   [[cγ (- sγ) 0 0]
    [sγ    cγ  0 0]
    [0     0   1 0]
    [0     0   0 1]]))

(define (rotation-y β)
  (define cβ (cos β)) (define sβ (sin β))
  (flomat: 
   [[   cβ  0    sβ  0]
    [   0   1    0   0]
    [(- sβ) 0    cβ  0]
    [   0   0    0   1]]))

(define (rotation-x α)
  (define cα (cos α)) (define sα (sin α))
  (flomat: 
   [[1 0      0   0]
    [0 cα  (- sα) 0]
    [0 sα     cα  0]
    [0 0      0   1]]))


;                        x y z
(define (rotation/angles γ β α) ; unfortunate naming ...
  (define cα (cos α)) (define sα (sin α))
  (define cβ (cos β)) (define sβ (sin β))
  (define cγ (cos γ)) (define sγ (sin γ))
  (flomat: 
   [[(* cα cβ)  (- (* cα sβ sγ) (* sα cγ))  (+ (* cα sβ cγ) (* sα sγ)) 0]
    [(* sα cβ)  (+ (* sα sβ sγ) (* cα cγ))  (- (* sα sβ cγ) (* cα sγ)) 0]
    [   (- sβ)        (* cβ sγ)                   (* cβ cγ)            0]
    [0          0                           0                          1]]))

(define (angle v w)
  (acos (/ (dot v w) (* (norm v) (norm w)))))

(define (rotation . args)
  (match args
    [(list (? number?) (? number?) (? number?))
     (apply rotation/angles args)]
    [(list (? flomat? direction)) ; a vector
     (define α (angle (dir 1 0 0) direction))
     (define β (angle (dir 0 1 0) direction))
     (define γ (angle (dir 0 0 1) direction))
     (rotation/angles α β γ)]))

;;;
;;; Cube
;;;

(define-values (cube-vertices cube-faces)
  ; Model space coordinates (i.e. (0,0,0) is the center.
  (let* ([T (λ (p) ; put center in origo
              (times (translation -0.5 -0.5 -0.5) p))] 

         [A (T (point 0 0 0))]
         [B (T (point 1 0 0))]
         [C (T (point 1 0 1))]
         [D (T (point 0 0 1))]
         
         [E (T (point 0 1 0))]
         [F (T (point 1 1 0))]
         [G (T (point 1 1 1))]
         [H (T (point 0 1 1))])
    (define vertices  (list A B C D E F G H))
    (define-values (a b c d e f g h) (values 0 1 2 3 4 5 6 7))
    (define face1     (list a b c d))
    (define face2     (list e f g h))
    (define face3     (list b f g c))
    (define face4     (list a e h d))
    (define face5     (list a b f e))
    (define face6     (list d c g h))
    (define face7     (list a a c c))
    (define faces     (list face1 face2 face3 face4 face5 face6 #;face7))
    (values vertices faces)))

(define cube (apply augment cube-vertices))


;;;
;;; Model Space
;;;

; The model center is in local coordinates in origo.
; To put the model in world coordinates with the proper scaling and rotation,
; we get:

(define (model position direction scale)
  (define S (scaling     scale))
  (define R (rotation    direction))
  (displayln (list 'position position))
  (define T (translation position))
  (displayln (list 'S S))
  (displayln (list 'R R))
  (displayln (list 'T T))
  (times T R S))

;;; 
;;; View Space (aka Camera Space)
;;;

; (define (view position direction)
  


;;;
;;; Projections
;;;

; The camera is pointed in direction (0,0,-1).
; It sees a cuboid   width x height x (far - near).

(define (ortographic-projection w h near far)
  ; From View Space into Ortho Projected Space
  ; The result is with in the [-1;1]³ cube
  (define Δz (- far near))
  (define Σz (+ far near))
  (flomat: 
   [[(/ w) 0      0         0]
    [0     (/ h)  0         0]
    [0     0      (/ -2 Δz) (- (/ Σz Δz))]
    [0     0      0         1]]))

(define (perspective-projection left right top bottom near far)
  (define-values (l r t b n f) (values left right top bottom near far))
  (define w (- r l))
  (define h (- t b))  
  (flomat: 
   [[(/ (* 2 n) w)    0              (/ (+ r l) w)         0]
    [0                (/ (* 2 n) h)  (/ (+ t b) h)         0]
    [0                0              (/ (+ f n) (- n f))   (/ (* 2 f n) (- n f))]
    [0                0              -1                    0]]))

;;;
;;; Face Sorting
;;;

(define (col! A j)
  (define r (nrows A))
  (sub! A 0 j r 1))

(define (sort-faces S faces)
  ; the columns in S are points
  ; todo: the norm below assumes w=1, which might not be the case...
  (define indices (for/list   ([j (ncols S)]) j))
  (define norms   (for/vector ([j (ncols S)]) (norm (col! S j))))
  (define (dist f) (apply + (map (λ (i) (vector-ref norms i)) f)))
  (sort faces > #:key dist))


;;;
;;; OBJ Parser
;;;

(define (read-obj [port (current-input-port)])
  (define lines         (for/list ([line (in-lines port)])
                          (with-input-from-string line
                            (λ ()
                              (if (string-prefix? line "#")
                                  (list 'comment line)
                                  (for/list ([x (in-port read)]) x))))))
  (define vertice-lines (filter (λ (l) (and (pair? l) (eq? (first l) 'v))) lines))
  (define face-lines    (filter (λ (l) (and (pair? l) (eq? (first l) 'f))) lines))
  (define nverts        (length vertice-lines))
  (define nfaces        (length face-lines))
  (define S             (make-flomat 4 nverts))
  (for ([v vertice-lines] [j (in-naturals)])
    (match v
      [(list 'v x y z)
       (mset! S 0 j x)
       (mset! S 1 j y)
       (mset! S 2 j z)
       (mset! S 3 j 1)]))
  (define (convert-triple t)
    ; input: 1/2/3
    ; we only use the vertex number (and ignore the texture vertex and
    ; vertex normal indices.)
    (first (map sub1 (map string->number (string-split (symbol->string t) "/")))))
    
  (define faces
    (for/list ([fl face-lines])
      (match fl
        [(list 'f (? number? a) ...) (map sub1 a)]
        [(list 'f (? symbol? a) ...) (map convert-triple a)])))        
  (values S faces))
                          



;;;
;;; Example
;;;

;;;
;;; MetaPict Rendering
;;;

;; (require (except-in metapict shape norm mat? dot dir angle))

;; (define pi (* 4 (atan 1)))

;; (define (polygon pts)
;;   (define (line p q) (bez p p q q))
;;   (define bezs (let loop ([pts pts])
;;                  (match pts
;;                    ['()           '()]
;;                    [(list  a)     '()]
;;                    [(list* a b d) (cons (line a b) (loop (cons b d)))])))
;;   (curve: #t bezs))
                     



;(set-curve-pict-size 100 100)
#;(apply
   beside 
   (for/list ([v (in-range 0 (+ pi (/ pi 8)) (/ pi 8))])
     (with-window (window -1.1 1.1 -1.1 1.1)
       (define d     (dir (cos v) (sin v) 0.5))
       ;                    pos,          dir, scale
       (define M     (model (point 0 0 -10) d   4))    ; model -> world
       ; (define P     (ortographic-projection 3        3        -2  2)) ; world -> device
       (define P     (perspective-projection -1.5 1.5 -1.5 1.5  -4  4 )) ; world -> device
       (define PM    (times P M))
       (define S     (times PM cube)) ; columns in curve are points
       (define vs    (for/vector ([j (ncols S)]) (define p (col S j)) (pt (xc p) (yc p))))
       (define faces  (sort-faces S cube-faces))
       (for/draw ([f (in-list faces)])    
         (define pts (for/list ([i f]) (vector-ref vs i)))
         (define c (curve* (append (add-between pts --) (list -- 'cycle))))
         (draw (color "white" (filldraw c))
               c)))))


(define-values (teapot teapot-faces) (read-obj (open-input-file "assets/teapot.obj")))
;(define-values (teapot teapot-faces) (read-obj (open-input-file "assets/teapot-low.obj")))

(define xmax (for/fold ([m -inf.0]) ([j (ncols teapot)]) (max m (ref teapot 0 j))))
(define ymax (for/fold ([m -inf.0]) ([j (ncols teapot)]) (max m (ref teapot 1 j))))
(define zmax (for/fold ([m -inf.0]) ([j (ncols teapot)]) (max m (ref teapot 2 j))))
(define xmin (for/fold ([m +inf.0]) ([j (ncols teapot)]) (min m (ref teapot 0 j))))
(define ymin (for/fold ([m +inf.0]) ([j (ncols teapot)]) (min m (ref teapot 1 j))))
(define zmin (for/fold ([m +inf.0]) ([j (ncols teapot)]) (min m (ref teapot 2 j))))
(list xmin xmax)
(list ymin ymax)
(list zmin zmax)
(define width  (* 2 (max (abs xmin) (abs xmax))))
(define height (* 2 (max (abs ymin) (abs ymax))))
(define depth  (* 2 (max (abs zmin) (abs zmax))))


; (set-curve-pict-size 600 600)

#;(with-window (window -1.1 1.1 -1.1 1.1)
  ; (define d     (dir (cos v) (sin v) 0.5))
  (define d     (dir 0 1 0))
  ;                    pos,          dir, scale
  (define M     (model (point 0 0 (- -10 7)) d   1))    ; model -> world
  (define P     (ortographic-projection   width     height   -10  0.5)) ; world -> device
  ; (define P     (perspective-projection -3.5 3.5 -3.5 3.5  -4  4 )) ; world -> device
  (define PM    (times P M))
  (define S     (times PM teapot)) ; columns in curve are points
  (define vs    (for/vector ([j (ncols S)]) (define p (col S j)) (pt (xc p) (yc p))))
  (define faces  (sort-faces S teapot-faces))
  (for/draw ([f (in-list faces)])    
    (define pts (for/list ([i f]) (vector-ref vs i)))
    (define c (polygon pts))
    (draw (color "white" (filldraw c))
          c)))

(require (except-in pict3d trace dir))

; (define d     (dir (cos v) (sin v) 0.5))
(define d     (dir -1 -1 0))
;                    pos,           dir, scale
(define M     (model (point 0 0 -4) d    1))                          ; model -> world
(displayln (list 'M M))
(define P     (ortographic-projection   width     height   -10  0.5)) ; world -> device
; (define P     (perspective-projection -3.5 3.5 -3.5 3.5  -4  4 ))   ; world -> device
(displayln (list 'P P))
(define PM    (times P M))
(displayln (list 'PM PM))

(define S     (times PM teapot)) ; columns in curve are points
(define vs    (for/vector ([j (ncols S)])
               (define p (col S j))
               (define x (xc p))
               (define y (yc p))
               (define z (zc p))
               (if (member +nan.0 (list x y z))
                   (begin (display ".") (pos 0 0 0))
                   (pos x y z))))
(define faces (sort-faces S teapot-faces)) ; no need to sort fact in pict3d

(define p
 (combine
  (for/list ([f (in-list teapot-faces)])
    (define pts (for/list ([i f]) (vector-ref vs i)))
    (case (length pts)
      [(3) (apply triangle pts #:back? #t)]
      [(4) (apply quad     pts #:back? #f)]
      [else (error)]))
  (light (pos 0 2 2) (emitted "white" 2))))

(current-pict3d-fov 80)
(current-pict3d-add-grid? #t)
(current-pict3d-add-wireframe 'color)
(pict3d->bitmap p 600 600)