apotome.lisp

   1 (in-package :run)
   2
   3
   4 ;; calculate e field
   5 (let* ((dx .1d0)
   6              (dz .1d0)
   7              (x 100)
   8              (y x)
   9              (z 100))
  10         (multiple-value-bind (ex ey ez)
  11        (psf:electric-field-psf z y x (* dz z) (* dx x) :integrand-evaluations 200)
  12           (defparameter *kex* (fftshift3 (ft3 ex)))
  13           (defparameter *key* (fftshift3 (ft3 ey)))
  14           (defparameter *kez* (fftshift3 (ft3 ez)))
  15           nil))
  16 #+Nil
  17 (progn
  18  (save-stack-ub8 "/home/martin/tmp/kex" (normalize3-cdf/ub8-realpart *kex*))
  19  (save-stack-ub8 "/home/martin/tmp/key" (normalize3-cdf/ub8-realpart *key*))
  20  (save-stack-ub8 "/home/martin/tmp/kez" (normalize3-cdf/ub8-realpart *kez*))
  21  )
  22 #+NIL
  23 (write-pgm "/home/martin/tmp/kex.pgm"
  24            (normalize2-cdf/ub8-abs (cross-section-xz *kex*)))
  25
  26 ;; calculate volume containing a grating
  27 (destructuring-bind (z y x)
  28      (array-dimensions *kex*)
  29    (let* ((grat (make-array (list z y x) :element-type '(complex double-float))))
  30      (let ((k (floor z 2))
  31            (q 30))
  32        (do-rectangle (j i q (- y q) q (- x q))
  33          (setf (aref grat k j i) (complex (* 1d0 (mod i 8))))))
  34      (defparameter *kgrat* (fftshift3 (ft3 grat)))
  35      nil))
  36 #+nil
  37 (save-stack-ub8 "/home/martin/tmp/kgrat" (normalize3-cdf/ub8-abs *kgrat*))
  38
  39 ;; calculate image of grating
  40 (destructuring-bind (z y x)
  41     (array-dimensions *kex*)
  42  (let* ((ex (fftshift3 (ift3 (fftshift3 (.* *kgrat* *kex*)))))
  43         (ey (fftshift3 (ift3 (fftshift3 (.* *kgrat* *key*)))))
  44         (ez (fftshift3 (ift3 (fftshift3 (.* *kgrat* *kez*)))))
  45         (intens (make-array (array-dimensions ex)
  46                             :element-type '(complex double-float))))
  47    (do-box (k j i 0 z 0 y 0 x)
  48      (setf (aref intens k j i)
  49            (+ (* (conjugate (aref ex k j i)) (aref ex k j i))
  50               (* (conjugate (aref ey k j i)) (aref ey k j i))
  51               (* (conjugate (aref ez k j i)) (aref ez k j i)))))
  52    (defparameter *intens* intens)
  53    (save-stack-ub8 "/home/martin/tmp/intens-grat" (normalize3-cdf/ub8-realpart intens))
  54    (write-pgm "/home/martin/tmp/intens-grat.pgm" (normalize2-cdf/ub8-realpart (cross-section-xz intens)))))
  55
  56 ;; store laser psf
  57 (destructuring-bind (z y x)
  58     (array-dimensions *kex*)
  59  (let* ((ex (ift3 (fftshift3 *kex*)))
  60         (ey (ift3 (fftshift3 *key*)))
  61         (ez (ift3 (fftshift3 *kez*)))
  62         (intens (make-array (array-dimensions ex)
  63                             :element-type '(complex double-float))))
  64    (do-box (k j i 0 z 0 y 0 x)
  65      (setf (aref intens k j i)
  66            (+ (* (conjugate (aref ex k j i)) (aref ex k j i))
  67               (* (conjugate (aref ey k j i)) (aref ey k j i))
  68               (* (conjugate (aref ez k j i)) (aref ez k j i)))))
  69    (defparameter *intens-psf-laser* intens)
  70    (save-stack-ub8 "/home/martin/tmp/intens-psf-laser" (normalize3-cdf/ub8-realpart intens))))
  71 #+nil
  72 (write-pgm
  73  "/home/martin/tmp/intens-psf-laser.pgm"
  74  (normalize2-cdf/ub8-realpart (cross-section-xz *intens-psf-laser*)))
  75
  76 #+nil
  77 (time
  78  (save-stack-ub8 "/home/martin/tmp/kintens-grat"
  79                  (normalize3-cdf/ub8-abs (fftshift3 (ft3 *intens*)))))
  80
  81 ;; A circular window in the center of the BFP with radius Rap gives
  82 ;; rise to rays up to the angle beta into sample space. The radius of
  83 ;; the focal sphere is n*f. Therefore one can write
  84 ;; sin(beta)=Rap/(n*f). Changing illumination direction of the grating
  85 ;; will shear the intensity image. In order to generate an image of
  86 ;; limited coherence one has to convolve each plane with a disk. The
  87 ;; radius of the disk is: Rd(z)=z*sin(beta) with defocus z.
  88 ;; Eliminate sin(beta):
  89 ;; Rd(z)=abs(z*Rap/(n*f))
  90 ;; for a 63x objective we get f=2.61, with n=1.515
  91
  92 (defun draw-disk (radius y x)
  93   (declare (double-float radius)
  94            (fixnum y x)
  95            (values (simple-array (complex double-float) 2) &optional))
  96   (let ((result (make-array (list y x) :element-type '(complex double-float)))
  97         (xh (floor x 2))
  98         (yh (floor y 2))
  99         (r2 (* radius radius)))
 100     (do-rectangle (j i 0 y 0 x)
 101       (let* ((ii (- i xh))
 102              (jj (- j yh))
 103              (rr2 (+ (* ii ii) (* jj jj))))
 104         (when (< rr2 r2)
 105           (setf (aref result j i) (complex (/ r2))))))
 106     result))
 107 #+nil
 108 (write-pgm "/home/martin/tmp/disk.pgm"
 109            (normalize2-cdf/ub8-realpart
 110             (fftshift2 (convolve2-circ
 111                         (draw-disk 12d0 100 200)
 112                         (draw-disk 12d0 100 200)))))
 113
 114 (sb-alien:define-alien-routine ("j1" bessel-j1)
 115     sb-alien:double
 116   (arg sb-alien:double))
 117
 118 ;; isi.ssl.berkeley.edu/~tatebe/whitepapers/FT%20of%20Uniform%20Disk.pdf
 119 (defun draw-unit-intensity-disk-precise (radius y x)
 120   (declare (double-float radius)
 121            (fixnum y x)
 122            (values (simple-array (complex double-float) 2) &optional))
 123   (let ((a (make-array (list y x)
 124                        :element-type '(complex double-float)))
 125         (xh (floor x 2))
 126         (yh (floor y 2)))
 127     (do-rectangle (j i 0 y 0 x)
 128       (let* ((xx (/ (* 1d0 (- i xh)) x))
 129              (yy (/ (* 1d0 (- j yh)) y))
 130              (f (sqrt (+ (* xx xx) (* yy yy)))))
 131         (setf (aref a j i) (if (< f 1d-12)
 132                                (complex (* pi radius))
 133                                (complex (/ (bessel-j1 (* 2d0 pi f radius))
 134                                            f))))))
 135     (fftshift2 (ift2 (fftshift2 a)))))
 136
 137 (defun draw-unit-energy-disk-precise (radius y x)
 138   (declare (double-float radius)
 139            (fixnum y x)
 140            (values (simple-array (complex double-float) 2) &optional))
 141   (let* ((disk (draw-unit-intensity-disk-precise radius y x))
 142          (sum (reduce #'+ (sb-ext:array-storage-vector disk))))
 143     (s*2 (/ (realpart sum)) disk)))
 144
 145 #+NIL
 146 (write-pgm "/home/martin/tmp/disk.pgm"
 147            (normalize2-cdf/ub8-realpart (draw-unit-energy-disk-precise 12.3d0 300 300)))
 148
 149 (time (destructuring-bind (z y x)
 150      (array-dimensions *intens*)
 151    (let ((n 1.515d0)
 152          (f (lens:focal-length-from-magnification 63d0))
 153          (Rap 4d0)
 154          (dz .21d0)
 155          (result (make-array (array-dimensions *intens*)
 156                              :element-type '(complex double-float))))
 157      ;; center plane doesn't have to be convolved
 158      (let ((k (floor z 2)))
 159        (do-rectangle (j i 0 y 0 x)
 160          (setf (aref result k j i) (aref *intens* k j i))))
 161
 162      ;; convolve all other planes with disks of different sizes
 163      (dotimes (k z)
 164        (unless (eq k (floor z 2))
 165          (let* ((Rd-pixels (abs (* (- k (floor z 2)) dz (/ Rap (* n f)))))
 166                 (disk (draw-unit-energy-disk-precise Rd-pixels y x))
 167                 (img (make-array (list y x) :element-type '(complex double-float))))
 168            (write-pgm (format nil "/home/martin/tmp/disk/disk~3,'0d.pgm" k)
 169                       (normalize2-cdf/ub8-realpart disk))
 170            (do-rectangle (j i 0 y 0 x)
 171              (setf (aref img j i) (aref *intens* k j i)))
 172            (let ((conv (fftshift2 (convolve2-circ img disk))))
 173              (do-rectangle (j i 0 y 0 x)
 174                (setf (aref result k j i) (+ (aref disk j i) (aref conv j i))))))))
 175      (defparameter *intens-disk* result)
 176      (save-stack-ub8 "/home/martin/tmp/intens-grat-conv" (normalize3-cdf/ub8-abs result))
 177      (save-stack-ub8 "/home/martin/tmp/kintens-grat-conv"
 178                      (normalize3-cdf/ub8-realpart (ft3 result)))
 179      (write-pgm "/home/martin/tmp/intens-grat-conv.pgm"
 180                 (normalize2-cdf/ub8-realpart (cross-section-xz *intens-disk*))))))