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