raytrace/raytrace.lisp

   1 (in-package :raytrace)
   2
   3 (declaim (optimize (speed 2) (safety 3) (debug 3)))
   4
   5 (define-condition one-solution () ())
   6 (define-condition no-solution () ())
   7
   8 (defun quadratic-roots (a b c)
   9   (declare (double-float a b c)
  10            (values double-float double-float &optional))
  11   "Find the two roots of ax^2+bx+c=0 and return them as multiple
  12   values."
  13   ;; see numerical recipes sec 5.6, p. 251 on how to avoid roundoff
  14   ;; error
  15   (let ((det2 (- (* b b) (* 4 a c))))
  16     (unless (<= 0d0 det2)
  17       (error 'no-solution))
  18     (let* ((pdet2 det2)
  19            (q (* .5d0 (+ b (* (signum b) (sqrt pdet2)))))
  20            (aa (abs a))
  21            (aq (abs q)))
  22       (declare ((double-float 0d0) pdet2))
  23       (cond ((and (< aq 1d-12) (< aa 1d-12)) (error 'no-solution))
  24             ((or (< aq 1d-12) (< aa 1d-12)) (error 'one-solution))
  25             (t (values (/ q a) (/ c q)))))))
  26
  27 #+nil ;; two solution
  28 (quadratic-roots 1d0 2d0 -3d0)
  29 #+nil ;; one solution
  30 (quadratic-roots 0d0 1d0 0d0)
  31 #+nil ;; no solution
  32 (quadratic-roots 0d0 -0d0 1d0)
  33
  34 (defclass sphere ()
  35   ((center :accessor center :initarg :center :initform (v) :type vec)
  36    (radius :accessor radius :initarg :radius :initform 1d0 :type double-float)))
  37
  38
  39 (defmethod print-object ((sphere sphere) stream)
  40   (with-slots (center radius) sphere
  41     (format stream "#<sphere radius: ~4f center: <~4f ~4f ~4f>>"
  42             radius (vec-x center) (vec-y center) (vec-z center))))
  43
  44 (defmethod ray-sphere-intersection-length ((ray ray) center radius)
  45   (declare (vec center)
  46            (double-float radius)
  47            (values double-float &optional))
  48   ;; (c-x)^2=r^2 defines the sphere, substitute x with the rays p+alpha a,
  49   ;; the raydirection should have length 1, solve the quadratic equation,
  50   ;; the distance between the two solutions is the distance that the ray
  51   ;; travelled through the sphere
  52   (check-direction-norm ray)
  53   (with-slots ((start vector::start) (direction vector::direction)) ray
  54     (let* ((l (v- center start))
  55            (c (- (v. l l) (* radius radius)))
  56            (b (* -2d0 (v. l direction))))
  57       (handler-case
  58           (multiple-value-bind (x1 x2)
  59               (quadratic-roots 1d0 b c)
  60             (abs (- x1 x2)))
  61         (one-solution () 0d0)
  62         (no-solution () 0d0)))))
  63
  64 #+nil
  65 (ray-sphere-intersection-length (v 0d0 .1d0 -12d0) (v 0d0 0d0 1d0) (v) 3d0)
  66
  67 (defmethod ray-spheres-intersection ((ray ray) (model sphere-algebraic-model)
  68                                      illuminated-sphere-index)
  69   (declare (fixnum illuminated-sphere-index)
  70            (values double-float &optional))
  71   (with-slots (centers-mm radii-mm) model
  72     (let ((sum 0d0))
  73       (loop for c in centers-mm and r in radii-mm and i from 0 do
  74            (unless (eq i illuminated-sphere-index)
  75             (incf sum (ray-sphere-intersection-length ray c r))))
  76       sum)))
  77
  78
  79
  80
  81