raytrace.lisp

   1 #.(progn
   2     (require :asdf)
   3     (require :vector)
   4 #+nil    (require :simple-gl)
   5     (require :simplex-anneal))
   6
   7 ;; for i in `cat raytrace.lisp|grep defun|grep -v '^;'|cut -d " " -f2`;
   8 ;; do echo \#:$i ;done
   9
  10
  11 (defpackage :raytrace
  12 #+nil  (:shadowing-import-from :cl close get special)
  13 (:export #:quadratic-roots
  14          #:ray-sphere-intersection-length
  15          #:direction
  16          #:ray-spheres-intersection
  17          #:sphere
  18          #:center
  19          #:radius
  20          #:make-sphere
  21          #:sphere-center
  22          #:sphere-radius)
  23 (:use :cl #+nil :gl #+nil :glut :vector :simplex-anneal))
  24
  25 (in-package :raytrace)
  26
  27 (declaim (optimize (speed 2) (safety 3) (debug 3)))
  28
  29 (declaim (ftype (function (double-float double-float double-float)
  30                           (values (or null double-float)
  31                                   &optional double-float))
  32                 quadratic-roots))
  33 (defun quadratic-roots (a b c)
  34   "Find the two roots of ax^2+bx+c=0 and return them as multiple
  35   values."
  36   ;; see numerical recipes sec 5.6, p. 251 on how to avoid roundoff
  37   ;; error
  38   (let ((det2 (- (* b b) (* 4 a c))))
  39     #+nil (let ((tiny 1d-12))
  40       (when (or (< (abs a) tiny)
  41                 (< (abs c) tiny)
  42                 (< (abs b) tiny))
  43        ;; i get division by zero when centered sphere
  44         ;; note: I didn't implement all the corner cases I think
  45        (return-from quadratic-roots nil #+nil (values 0d0 0d0))))
  46     (when (<= 0d0 det2) ;; we don't want complex results
  47       (let* ((pdet2 det2)
  48              (q (* .5d0 (+ b (* (signum b) (sqrt pdet2)))))
  49              (x1 (/ q a))
  50              (x2 (/ c q)))
  51         (declare ((double-float 0d0) pdet2)) ;; its positive
  52         (values x1 x2)))))
  53
  54 #+nil
  55 (quadratic-roots 1d0 2d0 -3d0)
  56 #+nil
  57 (quadratic-roots 0d0 -0d0 0d0)
  58
  59 (declaim (ftype (function (vec vec vec double-float)
  60                           (values (or null double-float) &optional))
  61                 ray-sphere-intersection-length))
  62 (defun ray-sphere-intersection-length
  63     (ray-start ray-direction sphere-center sphere-radius)
  64   ;; (c-x)^2=r^2 defines the sphere, substitute x with the rays p+alpha a,
  65   ;; the raydirection should have length 1, solve the quadratic equation,
  66   ;; the distance between the two solutions is the distance that the ray
  67   ;; travelled through the sphere
  68   (let* ((l (v- sphere-center ray-start))
  69          (c (- (v. l l) (* sphere-radius sphere-radius)))
  70          (a (normalize ray-direction))
  71          (b (* -2d0 (v. l a))))
  72     (multiple-value-bind (x1 x2)
  73         (quadratic-roots 1d0 b c)
  74       (when x1
  75         (abs (- x1 x2))))))
  76
  77 #+nil
  78 (ray-sphere-intersection-length (v 0d0 .1d0 -12d0) (v 0d0 0d0 1d0) (v) 3d0)
  79
  80 (defun direction (theta-degree phi-degree)
  81   "Convert spherical coordinates into cartesian."
  82   (let* ((theta (* theta-degree (/ pi 180d0)))
  83          (st (sin theta))
  84          (phi (* phi-degree (/ pi 180d0))))
  85     #+nil    (declare ((double-float 0d0 #.pi) theta) ;; ranges were a bit complicated
  86                       ((double-float 0d0 #.pi) phi))
  87     (v (* st (cos phi))
  88        (* st (sin phi))
  89        (* (cos theta)))))
  90
  91 #+nil
  92 (direction 45 0)
  93
  94
  95 (defstruct sphere
  96   (center (v) :type vec)
  97   (radius 0d0 :type double-float))
  98
  99 (declaim (ftype (function (vec vec (simple-array sphere 1) fixnum)
 100                           (values double-float &optional))
 101                 ray-spheres-intersection))
 102 (defun ray-spheres-intersection (ray-position ray-dir spheres
 103                                  illuminated-sphere-index)
 104  (let ((sum 0d0))
 105    (dotimes (i (length spheres))
 106      (unless (eq i illuminated-sphere-index)
 107        (with-slots (center radius)
 108            (aref spheres i)
 109         (let ((len (ray-sphere-intersection-length
 110                     ray-position ray-dir
 111                     center radius)))
 112           (when len
 113             (incf sum len))))))
 114    sum))
 115
 116 #+nil
 117 (defparameter centers
 118   '(( 6 87 157) ( 7 69 111) ( 7 87 196) ( 7 88 66) ( 7 92 33) ( 7 137 224)
 119     ( 7 144 149) ( 8 41 163) ( 8 61 201) ( 8 110 123)( 8 126 180)
 120     ( 8 131 99)( 9 34 107) ( 9 81 238) ( 10 166 208)( 11 58 77)  ( 11 111 245)
 121     ( 11 143 75)  ( 11 164 112)( 11 167 173)  ( 12 72 142) ( 12 99 76)
 122     ( 12 100 156) ( 12 122 44)  ( 12 144 248) ( 13 108 200)   ( 13 136 137)
 123     ( 14 138 221)   ( 15 50 219) ( 15 113 106)   ( 15 175 145)
 124     ( 16 56 107)  ( 16 77 189) ( 17 37 140)   ( 17 45 179)
 125     ( 17 89 51)   ( 17 162 187)( 18 88 240)    ( 18 129 171)
 126     ( 18 130 64) ( 18 163 188)( 19 151 99) ( 20 61 68)
 127     ( 21 81 139)  ( 21 86 97)( 21 136 210) ( 21 151 141)
 128     ( 22 48 109)( 22 53 151) ( 22 99 219)  ( 23 67 206)
 129     ( 23 98 187)   ( 23 107 72)( 25 107 146)( 25 113 111)))
 130 #+nil
 131 (progn
 132  (defparameter *central-sphere* 22)
 133  (defparameter *spheres*
 134    (make-array (length centers)
 135                :element-type 'sphere
 136                :initial-contents
 137                (let* ((q (elt centers *central-sphere*))
 138                       (s (/ 70d0))
 139                       (cen (v (* s (third q))
 140                               (* s (second q))
 141                               (* 5d0 s (first q)))))
 142                  (loop for c in centers collect
 143                       (make-sphere :center
 144                                    (v- (v (* s (third c))
 145                                           (* s (second c))
 146                                           (* 5d0 s (first c)))
 147                                        cen)
 148                                    :radius (* s 22d0)))))
 149    "center diameter"))
 150
 151
 152 #+nil
 153 (progn
 154    (declaim (ftype (function ((array double-float (2)))
 155                              (values double-float &optional))
 156                    merit-function))
 157    (defun merit-function (vec)
 158      (ray-spheres-intersection
 159       (v .1d0 .2d0 0d0)
 160       (normalize (direction (aref vec 0) (aref vec 1)))
 161       *spheres*
 162       *central-sphere*)))
 163
 164 #+nil
 165 (let ((start (make-array 2 :element-type 'double-float
 166                          :initial-contents (list 100d0 100d0))))
 167   (with-open-file (*standard-output* "/dev/shm/anneal.log"
 168                                      :direction :output
 169                                      :if-exists :supersede)
 170     (anneal (make-simplex start 1d0)
 171             #'merit-function
 172             :start-temperature 12d0)))
 173
 174 ;; scan over the full parameters space
 175 #+nil
 176 (progn
 177  (with-open-file (s "/dev/shm/o.dat" :direction :output :if-exists :supersede)
 178    (let ((ntheta 60)
 179          (nphi 90))
 180      (dotimes (theta ntheta)
 181        (dotimes (phi  nphi)
 182          (let ((a (* 90 (/ theta ntheta)))
 183                (b (* 180 (/ phi nphi))))
 184           (format s "~f ~f ~f~%" a b
 185                   (ray-spheres-intersection (v) (direction a b)))))
 186       (terpri s))))
 187  (with-open-file (s "/dev/shm/p1.gp" :direction :output
 188                     :if-exists :supersede)
 189    (format s "set term posts; set outp \"/dev/shm/p~2,'0d.ps\";set hidden
 190 set title \"nucleus nr. ~d\"
 191 unset key
 192 splot \"/dev/shm/o.dat\" u 1:2:3 w l
 193 #pause -1
 194 " *central-sphere* *central-sphere*)))