psf.lisp

   1 #.
   2 (require :vol)
   3
   4 (defpackage :psf
   5   (:use :cl :vol)
   6   (:export #:init
   7            #:intensity-psf
   8            #:electric-field-psf
   9            #:get-uv
  10            #:abs2))
  11
  12 (in-package :psf)
  13
  14 (declaim (optimize (speed 2) (debug 3) (safety 3)))
  15
  16 (progn
  17   (sb-alien:define-alien-routine j0
  18      sb-alien:double-float
  19    (x sb-alien:double-float))
  20
  21  (sb-alien:define-alien-routine j1
  22      sb-alien:double-float
  23    (x sb-alien:double-float))
  24
  25  (sb-alien:define-alien-routine jn
  26      sb-alien:double-float
  27    (n sb-alien:int)
  28    (x sb-alien:double-float))
  29
  30  (sb-alien:define-alien-routine j0f
  31      sb-alien:single-float
  32    (x sb-alien:single-float))
  33
  34  (sb-alien:define-alien-routine j1f
  35      sb-alien:single-float
  36    (x sb-alien:single-float))
  37
  38  (sb-alien:define-alien-routine jnf
  39      sb-alien:single-float
  40    (n sb-alien:int)
  41    (x sb-alien:single-float)))
  42
  43 (deftype my-float-helper ()
  44   `single-float)
  45
  46 (deftype my-float (&optional (low '*) (high '*))
  47   `(single-float ,(if (eq low '*)
  48                       '*
  49                       (coerce low 'my-float-helper))
  50                  ,(if (eq high '*)
  51                       '*
  52                       (coerce high 'my-float-helper))))
  53
  54 (defun abs2 (z)
  55   (declare ((complex my-float) z)
  56            (values my-float &optional))
  57   (let* ((x (realpart z))
  58          (y (imagpart z)))
  59     (+ (* x x) (* y y))))
  60
  61
  62 (declaim (type fixnum n)
  63          (type my-float /sin-alpha /sin-alpha2  sin-alpha2)
  64          (type (my-float -1 1) sin-alpha)
  65          (type (simple-array my-float 1) ac as as^2 ac2))
  66
  67 (progn
  68  (defparameter n 31)
  69  (defparameter ac (make-array n :element-type 'my-float))
  70  (defparameter as (make-array n :element-type 'my-float))
  71  (defparameter as^2 (make-array n :element-type 'my-float))
  72  (defparameter ac2 (make-array n :element-type 'my-float))
  73  (defparameter sin-alpha (coerce .4 'my-float))
  74  (defparameter sin-alpha2 (* sin-alpha sin-alpha))
  75  (defparameter /sin-alpha (/ sin-alpha))
  76  (defparameter /sin-alpha2 (/ sin-alpha2)))
  77
  78 (defun init (&key (numerical-aperture (coerce 1.2 'my-float))
  79              (immersion-index (coerce 1.515 'my-float))
  80              (integrand-evaluations 31))
  81   (declare (my-float numerical-aperture immersion-index)
  82            (fixnum integrand-evaluations)
  83            (values null &optional))
  84   (setf n integrand-evaluations
  85         ac (make-array n :element-type 'my-float)
  86         as (make-array n :element-type 'my-float)
  87         as^2 (make-array n :element-type 'my-float)
  88         ac2 (make-array n :element-type 'my-float)
  89         sin-alpha (/ numerical-aperture immersion-index)
  90         sin-alpha2 (* sin-alpha sin-alpha)
  91         /sin-alpha (/ sin-alpha)
  92         /sin-alpha2 (/ sin-alpha2))
  93   (let* ((alpha (asin sin-alpha))
  94          (dalpha (/ alpha n)))
  95
  96     (dotimes (iter n)
  97       (let* ((theta (the (my-float 0d0 2d0) (* dalpha iter)))
  98             (ct (cos theta))
  99              (st (sin theta))
 100              (ct2 (sqrt ct))
 101             (st^2 (* st st)))
 102         (declare ((my-float 0 1) ct)) ;; for the sqrt
 103         (setf (aref ac iter) ct
 104               (aref as iter) st
 105               (aref as^2 iter) st^2
 106               (aref ac2 iter) ct2)))))
 107 #+nil
 108 (init)
 109
 110 ;; k = omega/c = 2 pi / lambda
 111 ;; u = k z sin(alpha)^2
 112 ;; v = k sqrt(x^2+y^2) sin(alpha)
 113
 114 ;; for small angular aperture alpha << 1
 115 ;; u ~ k (a/R)^2 z
 116 ;; v ~ k (a/R) sqrt(x^2+y^2)
 117 ;; with a .. radius of the exit pupil and
 118 ;;      R .. distance between exit pupil and focal plane
 119
 120 (defun transform-xyz-to-uv (x y z &key (immersion-index (coerce 1.515 'my-float))
 121                             (numerical-aperture (coerce 1.38 'my-float))
 122                             (wavelength (coerce .480 'my-float)))
 123   "Return the cylindrical coordinates u and v of a given 3D point. All
 124 lengths in micrometer."
 125   (declare (my-float x y z immersion-index numerical-aperture wavelength)
 126            (values my-float my-float &optional))
 127   (let* ((k (/ (* #.(coerce 2 'my-float) (coerce pi 'my-float) immersion-index)
 128                wavelength))
 129          (sina (/ numerical-aperture
 130                   immersion-index))
 131          (u (* k z sina sina))
 132          (v (* k (sqrt (+ (* x x) (* y y))) sina)))
 133    (values u v)))
 134
 135 #+nil
 136 (transform-xyz-to-uv 0.0 1.0 1.0)
 137
 138 (defun intermediate-integrals-point (u v)
 139   (declare (my-float u v)
 140            (values (complex my-float) (complex my-float)
 141                    (complex my-float) &optional))
 142   (let* ((zero #.(coerce 0 '(complex my-float)))
 143          (i0 zero)
 144          (i1 zero)
 145          (i2 zero))
 146      (let ((scale 1))
 147        (dotimes (iter n)
 148          (let* ((s (aref as iter))
 149                 (c (aref ac iter))
 150                 (c2 (aref ac iter))
 151                 (vv (* v s /sin-alpha))
 152                 (uu (* u c /sin-alpha2))
 153                 (e (exp (complex 0 uu)))
 154                 (cmul0 (* scale e (* c2 s (+ 1 c) (j0 vv))))
 155                 (cmul1 (* scale e (* c2 (aref as^2 iter) (j1 vv))))
 156                 (cmul2 (* scale e (* c2 s (- 1 c) (jn 2 vv)))))
 157            (incf i0 cmul0)
 158            (incf i1 cmul1)
 159            (incf i2 cmul2)
 160            (cond ((eq iter 0) (setf scale 2))
 161                  ((eq iter (- n 2)) (setf scale 1))
 162                  ((eq scale 2) (setf scale 4))
 163                  ((eq scale 4) (setf scale 2))))))
 164      (let* ((s (* n 3d0)))
 165        (values (* s i0) (* s i1) (* s i2)))))
 166
 167 (defun energy-density (u v)
 168   (declare (my-float u v)
 169            (values my-float &optional))
 170   (multiple-value-bind (i0 i1 i2)
 171       (intermediate-integrals-point u v)
 172     (+ (abs2 i0)
 173        (abs2 i1)
 174        (abs2 i2))))
 175
 176 (defun energy-density-cyl (&optional (nu 100) (nv 100) (du #.(coerce .1 'my-float))
 177                            (dv du))
 178   (declare (fixnum nu nv)
 179            (my-float du dv)
 180            (values (simple-array my-float 2) &optional))
 181   (let* ((res (make-array (list nu nv)
 182                           :element-type 'my-float)))
 183     (dotimes (iu nu)
 184       (let ((u (* iu du)))
 185        (dotimes (iv nv)
 186          (let ((v (* iv dv)))
 187            (setf (aref res iu iv) (energy-density u v))))))
 188     res))
 189
 190 (defun intermediate-integrals-cyl (&optional (nu 100) (nv 100)
 191                                    (du (coerce .1 'my-float)) (dv du))
 192   (declare (fixnum nu nv)
 193            (my-float du dv)
 194            (values (simple-array (complex my-float) 2)
 195                    (simple-array (complex my-float) 2)
 196                    (simple-array (complex my-float) 2) &optional))
 197   (let* ((dims (list nu nv))
 198          (ii0 (make-array dims :element-type '(complex my-float)))
 199          (ii1 (make-array dims :element-type '(complex my-float)))
 200          (ii2 (make-array dims :element-type '(complex my-float))))
 201     (dotimes (iu nu)
 202       (let ((u (* iu du)))
 203        (dotimes (iv nv)
 204          (let ((v (* iv dv)))
 205            (multiple-value-bind (i0 i1 i2)
 206              (intermediate-integrals-point u v)
 207              (setf (aref ii0 iu iv) i0
 208                    (aref ii1 iu iv) i1
 209                    (aref ii2 iu iv) i2))))))
 210     (values ii0 ii1 ii2)))
 211
 212 #+nil
 213 (let* ((x 10)
 214        (y 10)
 215        (z 10)
 216        (size-x 3d0)
 217        (size-z 1d0)
 218        (nradius (1+ (ceiling (* (sqrt 2d0) (max x y)))))
 219        (nz (1+ (ceiling z 2))))
 220   (multiple-value-bind (u v)
 221       (transform-xyz-to-uv size-x 0 (* .5d0 size-z))
 222     (multiple-value-bind (i0 i1 i2)
 223         (intermediate-integrals-cyl nz nradius (/ u nz) (/ v nradius))
 224       (vol::interpolate2-cdf i0 2d0 2d0))))
 225
 226 ;; size radius is either the extend in x or y (in um) depending on what is bigger
 227 (defun electric-field-psf
 228     (z y x size-z size-radius &key (numerical-aperture (coerce 1.38 'my-float))
 229      (immersion-index (coerce 1.515 'my-float))
 230      (wavelength (coerce .480 'my-float))
 231      (integrand-evaluations 31))
 232   (declare (fixnum z y x integrand-evaluations)
 233            (my-float size-z size-radius numerical-aperture immersion-index
 234                      wavelength)
 235            (values (simple-array (complex my-float) 3)
 236                    (simple-array (complex my-float) 3)
 237                    (simple-array (complex my-float) 3) &optional))
 238   (init :numerical-aperture numerical-aperture
 239         :immersion-index immersion-index
 240         :integrand-evaluations integrand-evaluations)
 241   (let* ((nradius (1+ (ceiling (* (sqrt 2d0) (max x y)))))
 242          (nz (ceiling z 2))
 243          (dims (list z y x))
 244          (e0 (make-array dims :element-type '(complex my-float)))
 245          (e1 (make-array dims :element-type '(complex my-float)))
 246          (e2 (make-array dims :element-type '(complex my-float))))
 247     (multiple-value-bind (u v)
 248         (transform-xyz-to-uv 0 (* (sqrt #.(coerce 2 'my-float)) size-radius)
 249                              (* #.(coerce .5 'my-float) size-z)
 250                              :numerical-aperture numerical-aperture
 251                              :immersion-index immersion-index
 252                              :wavelength wavelength)
 253       (multiple-value-bind (i0 i1 i2)
 254           (intermediate-integrals-cyl nz nradius
 255                                       (/ (* 1d0 u) nz) (/ (* 1d0 v) nradius))
 256         (let ((rad-a (make-array (list y x) :element-type 'my-float))
 257               (cphi-a (make-array (list y x) :element-type 'my-float))
 258               (c2phi-a (make-array (list y x) :element-type 'my-float))
 259               (s2phi-a (make-array (list y x) :element-type 'my-float)))
 260           (do-region ((j i) (y x))
 261             (let* ((ii (- i (floor x 2)))
 262                    (jj (- j (floor y 2)))
 263                    (radius (sqrt (+ (* 1d0 ii ii) (* jj jj))))
 264                    (phi (atan (* 1d0 jj) ii)))
 265               (setf (aref rad-a j i) radius
 266                     (aref cphi-a j i) (cos phi)
 267                     (aref c2phi-a j i) (cos (* 2d0 phi))
 268                     (aref s2phi-a j i) (sin (* 2d0 phi)))))
 269           (let ((del (if (eq 1 (mod z 2))
 270                          #.(coerce 0 'my-float)
 271                          #.(coerce .5 'my-float)))) ;; add .5 if z is even
 272            (do-region ((k j i) (nz y x))
 273              (let* ((zi (- nz k (- 1d0 del)))
 274                     (r (aref rad-a j i))
 275                     (v0 (interpolate2 i0 zi r))
 276                     (v1 (interpolate2 i1 zi r))
 277                     (v2 (interpolate2 i2 zi r)))
 278                (setf (aref e0 k j i) (* (complex 0d0 -1d0)
 279                                         (+ v0 (* v2 (aref c2phi-a j i))))
 280                      (aref e1 k j i) (* (complex 0d0 -1d0)
 281                                         v2 (aref s2phi-a j i))
 282                      (aref e2 k j i) (* -2d0 v1 (aref cphi-a j i))))))
 283           (do-box (k j i 0 nz 0 y 0 x)
 284             (setf (aref e0 (- z k 1) j i) (- (conjugate (aref e0 k j i)))
 285                   (aref e1 (- z k 1) j i) (- (conjugate (aref e1 k j i)))
 286                   (aref e2 (- z k 1) j i) (conjugate (aref e2 k j i)))))))
 287     (values e0 e1 e2)))
 288
 289 #+nil
 290 (defparameter *e0* (electric-field 10 10 10 3d0 3d0))
 291
 292 (declaim (ftype (function (my-float my-float &key
 293                                         (:numerical-aperture my-float)
 294                                         (:immersion-index my-float)
 295                                         (:nz fixnum)
 296                                         (:nradius fixnum)
 297                                         (:wavelength my-float)
 298                                         (:integrand-evaluations fixnum))
 299                           (values (simple-array my-float 2) &optional))
 300                 intensity-psf-cyl))
 301 (defun intensity-psf-cyl (z radius &key (numerical-aperture 1.38d0)
 302                           (immersion-index 1.515d0)
 303                           (nz 50) (nradius 50) (wavelength .480d0)
 304                           (integrand-evaluations 31))
 305   "Calculate 2D cylindrical PSF with axial extend Z micrometers and
 306 transversal extend RADIUS micrometers."
 307   (psf:init :numerical-aperture numerical-aperture
 308             :immersion-index immersion-index
 309             :integrand-evaluations integrand-evaluations)
 310   (multiple-value-bind (u v)
 311       (transform-xyz-to-uv 0 radius z
 312                   :numerical-aperture numerical-aperture
 313                   :immersion-index immersion-index
 314                   :wavelength wavelength)
 315     (let* ((a (energy-density-cyl nz nradius (/ u nz) (/ v nradius))))
 316       a)))
 317 #+nil
 318 (time
 319  (progn (cyl-psf 3d0 1.5d0)
 320         nil))
 321
 322 (declaim (ftype (function (fixnum fixnum fixnum my-float my-float
 323                                   &key (:numerical-aperture my-float)
 324                                   (:immersion-index my-float)
 325                                   (:integrand-evaluations fixnum)
 326                                   (:wavelength my-float))
 327                           (values (simple-array (complex my-float) 3) &optional))
 328                 intensity-psf))
 329 (defun intensity-psf (z y x size-z size-radius &key (numerical-aperture 1.38d0)
 330                       (immersion-index 1.515d0) (integrand-evaluations 31)
 331                       (wavelength .480d0))
 332    "Calculate an intensity point spread function for an aplanatic microobjective with the given NUMERICAL-APERTURE, IMMERSION-INDEX and WAVELENGTH. Distances in micrometer."
 333    (let* ((psf (make-array (list z y x)
 334                            :element-type '(complex my-float)))
 335           (nz (1+ (ceiling z 2)))
 336           (nradius (1+ (ceiling (* (sqrt 2d0) (max x y)))))
 337           (cyl (intensity-psf-cyl
 338                 (* .5d0 size-z)
 339                 (* (sqrt 2d0)
 340                    size-radius)
 341                 :nz nz
 342                 :nradius nradius
 343                 :numerical-aperture numerical-aperture
 344                 :immersion-index immersion-index
 345                 :integrand-evaluations integrand-evaluations
 346                 :wavelength wavelength)))
 347      (let ((rad-a (make-array (list y x) :element-type 'my-float)))
 348        (do-rectangle (j i 0 y 0 x)
 349          (let* ((ii (- i (floor x 2)))
 350                 (jj (- j (floor y 2)))
 351                 (radius (sqrt (+ (* 1d0 ii ii) (* jj jj)))))
 352            (setf (aref rad-a j i) radius)))
 353        (let ((del (if (eq 1 (mod z 2)) ;; add .5 when z is even
 354                       0d0
 355                       .5d0)))
 356          (do-box (k j i 0 nz 0 y 0 x)
 357            (setf (aref psf k j i)
 358                  (complex (vol::interpolate2-df cyl
 359                                                 (-  nz k (- 1d0 del)) (aref rad-a j i))))))
 360        (do-box (k j i 0 nz 0 y 0 x)
 361          (setf (aref psf (- z k 1) j i) (aref psf k j i))))
 362      psf))
 363
 364
 365 #+nil
 366 (time (init))
 367 #+nil
 368 (time (progn (integ-all)
 369              nil))
 370
 371 ;; cylindrical coordinates:
 372 ;; v is rho
 373 ;; u is z
 374
 375 #+nil ;; print out a cross section of the psf transversal 1.5 um and
 376       ;; axial 3 um wide.
 377 (let ((numerical-aperture 1.38d0)
 378       (immersion-index 1.515d0))
 379   (init :numerical-aperture numerical-aperture
 380         :immersion-index immersion-index)
 381   (multiple-value-bind (u v)
 382       (transform-xyz-to-uv 0 1.5 3
 383                            :numerical-aperture numerical-aperture
 384                            :immersion-index immersion-index)
 385    (let* ((nu 10)
 386           (nv 20)
 387           (a (energy-density-cyl nu nv  (/ u nu) (/ v nv)))
 388           (max (aref a 0 0))
 389           (scale (/ 99 max)))
 390      (destructuring-bind (uu vv)
 391          (array-dimensions a)
 392        (format t "~%")
 393        (dotimes (u uu)
 394          (dotimes (v vv)
 395            (format t "~2d" (truncate (* scale (aref a u v)))))
 396          (format t "~%"))))))