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 (defun draw-disk (radius y x)
  15   (declare (my-float radius)
  16            (fixnum y x)
  17            (values (simple-array (complex my-float) 2) &optional))
  18   (let ((result (make-array (list y x) :element-type '(complex my-float)))
  19         (xh (floor x 2))
  20         (yh (floor y 2))
  21         (r2 (* radius radius)))
  22     (do-rectangle (j i 0 y 0 x)
  23       (let* ((ii (- i xh))
  24              (jj (- j yh))
  25              (rr2 (+ (* ii ii) (* jj jj))))
  26         (when (< rr2 r2)
  27           (setf (aref result j i) (complex (/ r2))))))
  28     result))
  29 #+nil
  30 (write-pgm "/home/martin/tmp/disk.pgm"
  31            (normalize2-cdf/ub8-realpart
  32             (fftshift2 (convolve2-circ
  33                         (draw-disk 12d0 100 200)
  34                         (draw-disk 12d0 100 200)))))
  35
  36 (sb-alien:define-alien-routine ("j1" bessel-j1-df)
  37     sb-alien:double
  38   (arg sb-alien:double))
  39
  40 (sb-alien:define-alien-routine ("j1f" bessel-j1-sf)
  41     sb-alien:float
  42   (arg sb-alien:float))
  43
  44 ;; isi.ssl.berkeley.edu/~tatebe/whitepapers/FT%20of%20Uniform%20Disk.pdf
  45 (defun draw-unit-intensity-disk-precise (radius y x)
  46   (declare (my-float radius)
  47            (fixnum y x)
  48            (values (simple-array (complex my-float) 2) &optional))
  49   (let ((a (make-array (list y x)
  50                        :element-type '(complex my-float)))
  51         (xh (floor x 2))
  52         (yh (floor y 2))
  53         (bessel-j1 (etypecase one
  54                      (single-float #'bessel-j1-sf)
  55                      (double-float #'bessel-j1-df)))
  56         (tiny (coerce 1d-12 'my-float)))
  57     (do-rectangle (j i 0 y 0 x)
  58       (let* ((xx (/ (* one (- i xh)) x))
  59              (yy (/ (* one (- j yh)) y))
  60              (f (sqrt (+ (* xx xx) (* yy yy)))))
  61         (setf (aref a j i) (if (< f tiny)
  62                                (complex (* my-pi radius))
  63                                (complex (/ (bessel-j1 (* two my-pi f radius))
  64                                            f))))))
  65     (fftshift2 (ift2 (fftshift2 a)))))
  66
  67 (defun draw-unit-energy-disk-precise (radius y x)
  68   (declare (my-float radius)
  69            (fixnum y x)
  70            (values (simple-array (complex my-float) 2) &optional))
  71   (let* ((disk (draw-unit-intensity-disk-precise radius y x))
  72          (sum (reduce #'+ (sb-ext:array-storage-vector disk))))
  73     (s*2 (/ (realpart sum)) disk)))
  74
  75 #+NIL
  76 (write-pgm "/home/martin/tmp/disk.pgm"
  77            (normalize2-cdf/ub8-realpart (draw-unit-energy-disk-precise 12.3d0 300 300)))
  78
  79
  80 (defun draw-sphere-ub8 (radius z y x)
  81   (declare (my-float radius)
  82            (fixnum z y x)
  83            (values (simple-array (unsigned-byte 8) 3)
  84                    &optional))
  85   (let ((sphere (make-array (list z y x)
  86                             :element-type '(unsigned-byte 8))))
  87     (let ((xh (floor x 2))
  88           (yh (floor y 2))
  89           (zh (floor z 2))
  90           (radius2 (* radius radius)))
  91      (do-box (k j i 0 z 0 y 0 x)
  92        (let ((r2 (+ (expt (* one (- i xh)) 2)
  93                     (expt (* one (- j yh)) 2)
  94                     (expt (* one (- k zh)) 2))))
  95          (setf (aref sphere k j i)
  96                (if (< r2 radius2)
  97                    1 0)))))
  98     sphere))
  99 #+nil
 100 (count-non-zero-ub8 (draw-sphere-ub8 one 41 58 70))
 101 #+nil
 102 (let ((a (draw-sphere-ub8 one
 103                           4 4 4
 104                           ;;3 3 3
 105                           ;;41 58 70
 106                           ))
 107       (res ()))
 108  (destructuring-bind (z y x)
 109      (array-dimensions a)
 110    (do-box (k j i 0 z 0 y 0 x)
 111      (when (eq 1 (aref a k j i))
 112        (setf res (list k j i))))
 113    res))
 114
 115 (defun draw-oval-ub8 (radius z y x)
 116   (declare (my-float radius)
 117            (fixnum z y x)
 118            (values (simple-array (unsigned-byte 8) 3)
 119                    &optional))
 120   (let ((sphere (make-array (list z y x)
 121                             :element-type '(unsigned-byte 8)))
 122         (scale-z (coerce 5d0 'my-float)))
 123     (do-box (k j i 0 z 0 y 0 x)
 124       (let ((r (sqrt (+ (square (- i (* half x)))
 125                         (square (- j (* half y)))
 126                         (square (* scale-z (- k (* half z))))))))
 127         (setf (aref sphere k j i)
 128               (if (< r radius)
 129                   1
 130                   0))))
 131     sphere))