raytrace.lisp

   1 #.(progn
   2     (require :asdf)
   3 #+nil    (require :simple-gl)
   4     (require :simplex-anneal))
   5
   6 ;; for i in `cat raytrace.lisp|grep defun|grep -v '^;'|cut -d " " -f2`;
   7 ;; do echo \#:$i ;done
   8
   9
  10 (defpackage :raytrace
  11 #+nil  (:shadowing-import-from :cl close get special)
  12 (:export #:v
  13          #:v.
  14          #:v-
  15          #:v+
  16          #:quadratic-roots
  17          #:ray-sphere-intersection-length
  18          #:direction
  19          #:ray-spheres-intersection
  20          #:vec
  21          #:mat
  22          #:make-sphere)
  23 (:use :cl #+nil :gl #+nil :glut :simplex-anneal))
  24
  25 (in-package :raytrace)
  26
  27 (declaim (optimize (speed 2) (safety 3) (debug 3)))
  28
  29 (deftype vec ()
  30   `(simple-array double-float (3)))
  31 (deftype mat ()
  32   `(simple-array double-float (3 3)))
  33
  34 (declaim (ftype (function (&optional double-float double-float double-float)
  35                           (values vec &optional))
  36                 v))
  37 (defun v (&optional (x 0d0) (y 0d0) (z 0d0))
  38   (make-array 3 :element-type 'double-float
  39               :initial-contents (list x y z)))
  40
  41 (declaim (ftype (function (vec vec)
  42                           (values double-float &optional))
  43                 v.))
  44 (defun v. (a b)
  45   "Dot product between two vectors."
  46   (let ((sum 0d0))
  47     (declare (double-float sum))
  48     (dotimes (i 3)
  49       (incf sum (* (aref a i)
  50                    (aref b i))))
  51     sum))
  52
  53 (declaim (ftype (function (vec vec)
  54                           (values vec &optional))
  55                 v- v+))
  56 (defun v- (a b)
  57   "Subtract two vectors."
  58   (let ((result (v)))
  59     (dotimes (i 3)
  60       (setf (aref result i) (- (aref a i)
  61                                (aref b i))))
  62     result))
  63
  64 (defun v+ (a b)
  65   "Add two vectors."
  66   (let ((result (v)))
  67     (dotimes (i 3)
  68       (setf (aref result i) (+ (aref a i)
  69                                (aref b i))))
  70     result))
  71
  72 (declaim (ftype (function (vec double-float)
  73                           (values vec &optional))
  74                 v*))
  75 (defmethod v* (a scalar)
  76   "Multiply vector with scalar."
  77   (declare (double-float scalar)
  78            (vec a))
  79   (let* ((result (v)))
  80     (dotimes (i 3)
  81       (setf (aref result i) (* scalar (aref a i))))
  82     result))
  83
  84 (declaim (ftype (function (vec)
  85                           (values double-float &optional))
  86                 norm))
  87 (defmethod norm (a)
  88   "Length of a vector."
  89   (let ((l2 (v. a a)))
  90     (declare (type (double-float 0d0) l2)) ;; Otherwise warning with complex-double
  91     (sqrt l2)))
  92
  93 (declaim (ftype (function (vec)
  94                           (values vec &optional))
  95                 normalize))
  96 (defmethod normalize (a)
  97   "Rescale vector to unit length."
  98   (let ((len (norm a)))
  99     (if (eq len 0d0)
 100         (v 0d0 0d0 1d0)
 101         (v* a (/ len)))))
 102
 103 (declaim (ftype (function (double-float double-float double-float)
 104                           (values (or null double-float)
 105                                   &optional double-float))
 106                 quadratic-roots))
 107 (defun quadratic-roots (a b c)
 108   "Find the two roots of ax^2+bx+c=0 and return them as multiple
 109   values."
 110   ;; see numerical recipes sec 5.6, p. 251 on how to avoid roundoff
 111   ;; error
 112   (let ((det2 (- (* b b) (* 4 a c))))
 113     #+nil (let ((tiny 1d-12))
 114       (when (or (< (abs a) tiny)
 115                 (< (abs c) tiny)
 116                 (< (abs b) tiny))
 117        ;; i get division by zero when centered sphere
 118         ;; note: I didn't implement all the corner cases I think
 119        (return-from quadratic-roots nil #+nil (values 0d0 0d0))))
 120     (when (<= 0d0 det2) ;; we don't want complex results
 121       (let* ((pdet2 det2)
 122              (q (* .5d0 (+ b (* (signum b) (sqrt pdet2)))))
 123              (x1 (/ q a))
 124              (x2 (/ c q)))
 125         (declare ((double-float 0d0) pdet2)) ;; its positive
 126         (values x1 x2)))))
 127
 128 #+nil
 129 (quadratic-roots 1d0 2d0 -3d0)
 130 #+nil
 131 (quadratic-roots 0d0 -0d0 0d0)
 132
 133 (declaim (ftype (function (vec vec vec double-float)
 134                           (values (or null double-float) &optional))
 135                 ray-sphere-intersection-length))
 136 (defun ray-sphere-intersection-length
 137     (ray-start ray-direction sphere-center sphere-radius)
 138   ;; (c-x)^2=r^2 defines the sphere, substitute x with the rays p+alpha a,
 139   ;; the raydirection should have length 1, solve the quadratic equation,
 140   ;; the distance between the two solutions is the distance that the ray
 141   ;; travelled through the sphere
 142   (let* ((l (v- sphere-center ray-start))
 143          (c (- (v. l l) (* sphere-radius sphere-radius)))
 144          (a (normalize ray-direction))
 145          (b (* -2d0 (v. l a))))
 146     (multiple-value-bind (x1 x2)
 147         (quadratic-roots 1d0 b c)
 148       (when x1
 149         (abs (- x1 x2))))))
 150
 151 #+nil
 152 (ray-sphere-intersection-length (v 0d0 .1d0 -12d0) (v 0d0 0d0 1d0) (v) 3d0)
 153
 154 (defun direction (theta-degree phi-degree)
 155   (let* ((theta (* theta-degree (/ pi 180d0)))
 156          (st (sin theta))
 157          (phi (* phi-degree (/ pi 180d0))))
 158 #+nil    (declare ((double-float 0d0 #.pi) theta) ;; ranges were a bit complicated
 159              ((double-float 0d0 #.pi) phi))
 160     (v (* st (cos phi)) ;; convert spherical coordinates into cartesian
 161        (* st (sin phi))
 162        (* (cos theta)))))
 163
 164 #+nil
 165 (direction 45 0)
 166
 167 #+nil
 168 (monte-carlo)
 169
 170 #+nil (defparameter *rays* '())
 171
 172 #+nil (defparameter centers
 173   '(( 6 87 157) ( 7 69 111) ( 7 87 196) ( 7 88 66) ( 7 92 33) ( 7 137 224)
 174     ( 7 144 149) ( 8 41 163) ( 8 61 201) ( 8 110 123)( 8 126 180)
 175     ( 8 131 99)( 9 34 107) ( 9 81 238) ( 10 166 208)( 11 58 77)  ( 11 111 245)
 176     ( 11 143 75)  ( 11 164 112)( 11 167 173)  ( 12 72 142) ( 12 99 76)
 177     ( 12 100 156) ( 12 122 44)  ( 12 144 248) ( 13 108 200)   ( 13 136 137)
 178     ( 14 138 221)   ( 15 50 219) ( 15 113 106)   ( 15 175 145)
 179     ( 16 56 107)  ( 16 77 189) ( 17 37 140)   ( 17 45 179)
 180     ( 17 89 51)   ( 17 162 187)( 18 88 240)    ( 18 129 171)
 181     ( 18 130 64) ( 18 163 188)( 19 151 99) ( 20 61 68)
 182     ( 21 81 139)  ( 21 86 97)( 21 136 210) ( 21 151 141)
 183     ( 22 48 109)( 22 53 151) ( 22 99 219)  ( 23 67 206)
 184     ( 23 98 187)   ( 23 107 72)( 25 107 146)( 25 113 111)))
 185 #+nil (progn
 186  (defparameter *central-sphere* 22)
 187  (defparameter *spheres*
 188    (let* ((q (elt centers *central-sphere*))
 189           (s (/ 70d0))
 190           (cen (v (* s (third q))
 191                   (* s (second q))
 192                   (* 5d0 s (first q)))))
 193      (loop for c in centers collect
 194           (list (v- (v (* s (third c))
 195                        (* s (second c))
 196                        (* 5d0 s (first c)))
 197                     cen)
 198                 (* s 11d0))))
 199    "center diameter"))
 200
 201 (defstruct sphere
 202   (center (v) :type vec)
 203   (radius 0d0 :type double-float))
 204
 205 #+nil (dotimes (i 100)
 206   (setf *spheres* (append *spheres*
 207                           (list (list (v* (direction (random 180d0)
 208                                                      (random 360d0))
 209                                           (1+ (random 1d0)))
 210                                       (+ .15d0
 211                                          (random .18d0)))))))
 212 (declaim (ftype (function (vec vec (simple-array sphere 1) fixnum)
 213                           (values double-float &optional))
 214                 ray-spheres-intersection))
 215 (defun ray-spheres-intersection (ray-position ray-dir spheres
 216                                  illuminated-sphere-index)
 217  (let ((sum 0d0))
 218    (dotimes (i (length spheres))
 219      (unless (eq i illuminated-sphere-index)
 220        (with-slots (center radius)
 221            (aref spheres i)
 222         (let ((len (ray-sphere-intersection-length
 223                     ray-position ray-dir
 224                     center radius)))
 225           (when len
 226             (incf sum len))))))
 227    sum))
 228
 229 ;; (declaim (ftype (function (vec)
 230 ;;                        (values double-float &optional))
 231 ;;              merit-function))
 232 ;; (defun merit-function (vec)
 233 ;;   (ray-spheres-intersection
 234 ;;    (v .1d0 .2d0 0d0)
 235 ;;    (normalize (direction (aref vec 0) (aref vec 1)))))
 236
 237 #+nil
 238 (let ((start (make-array 2 :element-type 'double-float
 239                           :initial-contents (list 100d0 100d0))))
 240    (with-open-file (*standard-output* "/dev/shm/anneal.log"
 241                                       :direction :output
 242                                       :if-exists :supersede)
 243      (anneal (make-simplex start 1d0)
 244              #'merit-function
 245              :start-temperature 3d0)))
 246
 247 #+nil (RAY-SPHERE-INTERSECTION-LENGTH
 248  (v 0.0d0 0.0d0 0.0d0)
 249  (v 0.0d0 0.0d0 1.0d0)
 250  (v 0.0d0 0.0d0 0.0d0)
 251  0.15714285714285714d0)
 252
 253 ;; scan over the full parameters spcae
 254 #+nil
 255 (progn
 256  (with-open-file (s "/dev/shm/o.dat" :direction :output :if-exists :supersede)
 257    (let ((ntheta 60)
 258          (nphi 90))
 259      (dotimes (theta ntheta)
 260        (dotimes (phi  nphi)
 261          (let ((a (* 90 (/ theta ntheta)))
 262                (b (* 180 (/ phi nphi))))
 263           (format s "~f ~f ~f~%" a b
 264                   (ray-spheres-intersection (v) (direction a b)))))
 265       (terpri s))))
 266  (with-open-file (s "/dev/shm/p1.gp" :direction :output
 267                     :if-exists :supersede)
 268    (format s "set term posts; set outp \"/dev/shm/p~2,'0d.ps\";set hidden
 269 set title \"nucleus nr. ~d\"
 270 unset key
 271 splot \"/dev/shm/o.dat\" u 1:2:3 w l
 272 #pause -1
 273 " *central-sphere* *central-sphere*)))
 274
 275 #+nil
 276 (ray-spheres-intersection (v .1d0 .2d0 0d0)
 277                           (direction -1631995.127159034d0 -881148.5633918238d0))
 278
 279 ;; (direction -1631995.127159034d0 -881148.5633918238d0)
 280 ;; #(0.5991641440565072d0 -0.6787435869828116d0 -0.4246286278699771d0)
 281
 282 ;; smallest value in the plot
 283 ;; cat /dev/shm/o.dat |cut -d " " -f 3|sort -n|uniq -c|head
 284
 285
 286 ;; (declaim (ftype (function () (values (or null double-float) &optional))
 287 ;;              sum-random-dir))
 288 ;; (defun sum-random-dir ()
 289 ;;   (let ((start (v (random 1d0) (random 1d0) 1d0))
 290 ;;      (end (v 0d0 0d0 -1d0))
 291 ;;      (sum 0d0))
 292 ;;     (setf *rays* (append *rays* (list (list start end))))
 293 ;;     (dolist (s *spheres*)
 294 ;;       (let ((len (ray-sphere-intersection-length
 295 ;;                start (v- end start)
 296 ;;                (first s) (second s))))
 297 ;;      (when len
 298 ;;        (incf sum len))))
 299 ;;     sum))
 300
 301 ;; (declaim (ftype (function () (values double-float &optional))
 302 ;;              monte-carlo))
 303 ;; (defun monte-carlo ()
 304 ;;   ;; nr3 p. 422
 305 ;;   (let ((f 0d0)
 306 ;;      (f2 0d0)
 307 ;;      (rel-err 20d0)
 308 ;;      (val 0d0))
 309 ;;     (declare (double-float f f2 rel-err val))
 310 ;;     (do ((count 0))
 311 ;;      ((and (< 100 count)
 312 ;;            (< rel-err .02d0)))
 313 ;;       (declare (fixnum count))
 314 ;;       (let ((integ (sum-random-dir)))
 315 ;;      (when integ
 316 ;;        (incf f integ)
 317 ;;        (incf f2 (* integ integ))
 318 ;;        (incf count)
 319 ;;        (when (< 1 count)
 320 ;;          (setf val (/ f count))
 321 ;;          (setf rel-err (/ (sqrt (/ (- (/ f2 count) (* val val))
 322 ;;                                    (1+ count)))
 323 ;;                           val)))
 324 ;;       #+nil (format t "~a~%" (list count val rel-err)))))
 325 ;;     val))
 326
 327
 328 ;; (declaim (ftype (function (vec)
 329 ;;                        (values null &optional))
 330 ;;              vertex-v translate-v))
 331 ;; (defun vertex-v (vert)
 332 ;;   (vertex (aref vert 0) (aref vert 1) (aref vert 2))
 333 ;;   nil)
 334 ;; (defun translate-v (vert)
 335 ;;   (translate (aref vert 0) (aref vert 1) (aref vert 2))
 336 ;;   nil)
 337
 338 ;; (declaim (ftype (function (double-float vec)
 339 ;;                        (values null &optional))
 340 ;;              rotate-v))
 341 ;; (defun rotate-v (angle vert)
 342 ;;   (rotate angle (aref vert 0) (aref vert 1) (aref vert 2))
 343 ;;   nil)
 344
 345 ;; (defmacro with-light (&body body)
 346 ;;   `(progn
 347 ;;      (enable :light0 :light1 :lighting)
 348 ;;      (shade-model :flat)
 349 ;;      (material :front :specular #(1d0 1d0 1d0 1d0))
 350 ;;      (material :front :shininess 70d0)
 351 ;;      (light :light0 :position #(12d0 32d0 3d0 1d0))
 352 ;;      (light :light0 :ambient #(.3d0 .3d0 .3d0 1d0))
 353 ;;      (light :light1 :position #(12d0 -32d0 -3d0 1d0))
 354 ;;      (light :light1 :diffuse #(1d0 1d0 1d0 1d0))
 355 ;;      ,@body
 356 ;;      (disable :light0 :light1 :lighting)))
 357
 358 ;; (defparameter *rot* 0d0)
 359
 360 ;; (defun draw ()
 361 ;;   (if (< *rot* 360d0)
 362 ;;       (incf *rot*)
 363 ;;       (setf *rot* 0d0))
 364 ;;   (clear :depth-buffer-bit)
 365 ;;   (with-primitive :lines
 366 ;;     (color 1 0 0) (vertex 0 0 0) (vertex 1 0 0)
 367 ;;     (color 0 1 0) (vertex 0 0 0) (vertex 0 1 0)
 368 ;;     (color 0 0 1) (vertex 0 0 0) (vertex 0 0 1))
 369
 370 ;;   (rotate *rot* 0 0 1)
 371 ;;   (line-width 3)
 372 ;; #+nil  (enable :blend)
 373 ;;   (blend-func :src-alpha :one)
 374 ;;   (color 1 1 1 .2)
 375 ;;   (enable :depth-test)
 376 ;;   (with-light
 377 ;;       (let ((s 10))
 378 ;;      (dolist (e *spheres*)
 379 ;;       (with-pushed-matrix
 380 ;;         (translate-v (first e))
 381 ;;         (solid-sphere (second e) (* 2 s) s)))))
 382 ;;   (disable :depth-test)
 383 ;;   (color 1 1 1 .002)
 384 ;;   (with-primitive :lines
 385 ;;     (dolist (c *rays*)
 386 ;;       (vertex-v (first c))
 387 ;;       (vertex-v (second c))))
 388 ;;   (disable :blend)
 389 ;;   (scale 3 3 3)
 390 ;;   (with-primitive :lines
 391 ;;     (color 1 0 0) (vertex 0 0 0) (vertex 1 0 0)
 392 ;;     (color 0 1 0) (vertex 0 0 0) (vertex 0 1 0)
 393 ;;     (color 0 0 1) (vertex 0 0 0) (vertex 0 0 1)))
 394
 395 ;; (defun run ()
 396 ;;   (simple-gl:with-simple-gl
 397 ;;       (draw)))
 398
 399 ;; (eval-when (:execute)
 400 ;;  (run))