vol/draw.lisp

   1 (in-package :vol)
   2
   3 ;; A circular window in the center of the BFP with radius Rap gives
   4 ;; rise to rays up to the angle beta into sample space. The radius of
   5 ;; the focal sphere is n*f. Therefore one can write
   6 ;; sin(beta)=Rap/(n*f). Changing illumination direction of the grating
   7 ;; will shear the intensity image. In order to generate an image of
   8 ;; limited coherence one has to convolve each plane with a disk. The
   9 ;; radius of the disk is: Rd(z)=z*sin(beta) with defocus z.
  10 ;; Eliminate sin(beta):
  11 ;; Rd(z)=abs(z*Rap/(n*f))
  12 ;; for a 63x objective we get f=2.61, with n=1.515
  13
  14 (def-generator (draw-disk (type))
  15   `(defun ,name (radius y x)
  16      (declare (single-float radius)
  17               (fixnum y x)
  18               (values (simple-array ,long-type 2) &optional))
  19      (let ((result (make-array (list y x) :element-type ',long-type))
  20            (xh (floor x 2))
  21            (yh (floor y 2))
  22            (r2 (* radius radius)))
  23        (do-region ((j i) (y x))
  24          (let* ((ii (- i xh))
  25                 (jj (- j yh))
  26                 (rr2 (+ (* ii ii) (* jj jj))))
  27            (when (< rr2 r2)
  28              (setf (aref result j i) (* ,(case type
  29                                                (ub8 255)
  30                                                (otherwise (coerce 1 long-type))))))))
  31        result)))
  32
  33 #+nil
  34 (def-draw-disk-type ub8)
  35
  36 (defmacro def-draw-disk-functions (types)
  37   (let* ((specifics nil)
  38          (name (format-symbol "draw-disk")))
  39     (loop for type in types do
  40          (let ((def-name (format-symbol "def-~a-type" name)))
  41            (push `(,def-name ,type) specifics)))
  42     `(progn ,@specifics)))
  43
  44 (def-draw-disk-functions (ub8 sf df csf cdf))
  45
  46 #+nil
  47 (write-pgm "/home/martin/tmp/disk.pgm"
  48            (draw-disk-ub8 33.0 100 100))
  49
  50 #+nil
  51 (write-pgm "/home/martin/tmp/disk.pgm"
  52            (normalize2-cdf/ub8-realpart
  53             (fftshift2 (convolve2-circ
  54                         (draw-disk 12d0 100 200)
  55                         (draw-disk 12d0 100 200)))))
  56
  57 (sb-alien:define-alien-routine ("j1" bessel-j1-df)
  58     sb-alien:double
  59   (arg sb-alien:double))
  60
  61 (sb-alien:define-alien-routine ("j1f" bessel-j1-sf)
  62     sb-alien:float
  63   (arg sb-alien:float))
  64
  65 ;; isi.ssl.berkeley.edu/~tatebe/whitepapers/FT%20of%20Uniform%20Disk.pdf
  66 (def-generator (draw-uniform-disk-precise (type))
  67   (let* ((rtype (ecase type
  68                  (csf 'sf)
  69                  (cdf 'df)))
  70         (rlong-type (get-long-type rtype)))
  71     `(defun ,name (radius y x)
  72       (declare (single-float radius)
  73                (fixnum y x)
  74                (values (simple-array ,long-type 2) &optional))
  75       (let ((a (make-array (list y x)
  76                            :element-type ',long-type))
  77             (xh (floor x 2))
  78             (yh (floor y 2))
  79             (one ,(coerce 1 rlong-type))
  80             (tiny ,(coerce 1e-12 rlong-type)))
  81         (do-region ((j i) (y x))
  82           (let* ((xx (/ (* one (- i xh)) x))
  83                  (yy (/ (* one (- j yh)) y))
  84                  (f (sqrt (+ (* xx xx) (* yy yy)))))
  85             (setf (aref a j i) (if (< f tiny)
  86                                    (complex (* ,(coerce pi rlong-type) radius))
  87                                    (complex (/ (,(ecase type
  88                                                         (csf `bessel-j1-sf)
  89                                                         (cdf `bessel-j1-df))
  90                                                  (* ,(coerce (* 2 pi) rlong-type)
  91                                                     f radius))
  92                                                f))))))
  93         (fftshift (ift (fftshift a)))))))
  94
  95 (def-draw-uniform-disk-precise-type csf)
  96 (def-draw-uniform-disk-precise-type cdf)
  97
  98 (def-generator (draw-unit-energy-disk-precise (type))
  99   `(defun ,name (radius y x)
 100      (declare (single-float radius)
 101               (fixnum y x)
 102               (values (simple-array ,long-type 2) &optional))
 103      (let* ((disk (,(format-symbol "draw-uniform-disk-precise-~a" type)
 104                     radius y x))
 105             (sum (reduce #'+ (sb-ext:array-storage-vector disk))))
 106        (s* (/ (realpart sum)) disk))))
 107
 108 (def-draw-unit-energy-disk-precise-type csf)
 109 (def-draw-unit-energy-disk-precise-type cdf)
 110
 111 #+NIL
 112 (write-pgm "/home/martin/tmp/disk.pgm"
 113            (normalize-2-csf/ub8-realpart
 114             (draw-unit-energy-disk-precise-csf 12.3 300 300)))
 115
 116 (defun square-sf (x)
 117   (declare (single-float x)
 118            (values single-float &optional))
 119   (* x x))
 120
 121 (def-generator (draw-sphere (type))
 122  `(defun ,name (radius z y x)
 123     (declare (single-float radius)
 124              (fixnum z y x)
 125              (values (simple-array ,long-type 3)
 126                      &optional))
 127    (let ((sphere (make-array (list z y x)
 128                              :element-type ',long-type)))
 129      (let ((xh (floor x 2))
 130            (yh (floor y 2))
 131            (zh (floor z 2))
 132            (radius2 (* radius radius)))
 133        (do-region ((k j i) (z y x))
 134          (let ((r2 (+ (square-sf (* 1s0 (- i xh)))
 135                       (square-sf (* 1s0 (- j yh)))
 136                       (square-sf (* 1s0 (- k zh))))))
 137            (setf (aref sphere k j i)
 138                  (if (< r2 radius2)
 139                      ,(coerce 1 long-type)
 140                      ,(coerce 0 long-type))))))
 141      sphere)))
 142
 143 (defmacro def-draw-sphere-functions (types)
 144   `(progn ,@(loop for type in types collect
 145                  `(def-draw-sphere-type ,type))))
 146
 147 (def-draw-sphere-functions (ub8 sf df csf cdf))
 148
 149 #+nil
 150 (count-non-zero-ub8 (draw-sphere-ub8 1.0 41 58 70))
 151 #+nil
 152 (let ((a (draw-sphere-ub8 1.0
 153                           4 4 4
 154                           ;;3 3 3
 155                           ;;41 58 70
 156                           ))
 157       (res ()))
 158   (destructuring-bind (z y x)
 159       (array-dimensions a)
 160     (do-region ((k j i) (z y x))
 161       (when (eq 1 (aref a k j i))
 162         (setf res (list k j i))))
 163     res))
 164
 165 (def-generator (draw-oval (type))
 166  `(defun ,name (radius z y x)
 167    (declare (single-float radius) (fixnum z y x)
 168             (values (simple-array ,long-type 3) &optional))
 169    (let ((sphere (make-array (list z y x) :element-type ',long-type))
 170          (scale-z 5.0))
 171      (do-region ((k j i) (z y x))
 172        (let ((r (sqrt (+ (square-sf (- i (* .5 x)))
 173                          (square-sf (- j (* .5 y)))
 174                          (square-sf (* scale-z (- k (* .5 z))))))))
 175          (setf (aref sphere k j i)
 176                (if (< r radius)
 177                    ,(coerce 1 long-type)
 178                    ,(coerce 0 long-type)))))
 179      sphere)))
 180
 181 (defmacro def-draw-oval-functions (types)
 182   `(progn ,@(loop for type in types collect
 183                  `(def-draw-oval-type ,type))))
 184
 185 (def-draw-oval-functions (ub8 sf df csf cdf))