psf.lisp

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