run.lisp

   1 #.(progn (require :asdf)
   2          (require :vol)
   3          (require :psf)
   4          (require :simplex-anneal)
   5          (require :raytrace)
   6          (require :lens))
   7
   8 (defpackage :run
   9   (:use :cl :vol :raytrace))
  10 (in-package :run)
  11 #||
  12 mkdir ~/tmp
  13 cp /home/martin/0519/MedianofCelegans-10-02-09-LSM700-t58.tif ~/tmp/med.tif
  14 cd ~/tmp/
  15 tiffsplit med.tif
  16 for i in *.tif ; do tifftopnm $i > `basename $i .tif`.pgm;done
  17 ||#
  18
  19 (deftype vec-i ()
  20   `(simple-array fixnum (3)))
  21
  22 (defstruct (vec-i (:type (vector fixnum)))
  23   (z 0) (y 0) (x 0))
  24
  25 (declaim (ftype (function (vec-i vec-i)
  26                           (values fixnum &optional))
  27                 v.-i))
  28 (defun v.-i (a b)
  29   (+ (* (vec-i-x a) (vec-i-x b))
  30      (* (vec-i-y a) (vec-i-y b))
  31      (* (vec-i-z a) (vec-i-z b))))
  32
  33 (declaim (ftype (function (vec-i vec-i)
  34                           (values vec-i &optional))
  35                 v--i v+-i))
  36 (defun v--i (a b)
  37   (make-vec-i :x (- (vec-i-x a) (vec-i-x b))
  38               :y (- (vec-i-y a) (vec-i-y b))
  39               :z (- (vec-i-z a) (vec-i-z b))))
  40 (defun v+-i (a b)
  41   (make-vec-i :x (+ (vec-i-x a) (vec-i-x b))
  42               :y (+ (vec-i-y a) (vec-i-y b))
  43               :z (+ (vec-i-z a) (vec-i-z b))))
  44
  45 (declaim (ftype (function (vec-i)
  46                           (values double-float &optional))
  47                 norm-i))
  48 (defun norm-i (a)
  49   (sqrt (* 1d0 (v.-i a a))))
  50
  51 (declaim (ftype (function (string (simple-array (complex double-float) 3)
  52                                   &key (:function function))
  53                           (values null &optional))
  54                 save-stack))
  55 (defun save-stack (fn vol &key (function #'realpart))
  56   (ensure-directories-exist fn)
  57   (destructuring-bind (z y x)
  58       (array-dimensions vol)
  59     (let ((b (make-array (list y x)
  60                          :element-type '(unsigned-byte 8))))
  61       (dotimes (k z)
  62         (do-rectangle (j i 0 y 0 x)
  63           (setf (aref b j i)
  64                 (clamp (floor (* 255 (funcall function (aref vol k j i)))))))
  65         (write-pgm b (format nil "~a/~3,'0d.pgm" fn k)))))
  66   nil)
  67
  68 (declaim (ftype (function (string (simple-array (unsigned-byte 8) 3))
  69                           (values null &optional))
  70                 save-stack-ub8))
  71 (defun save-stack-ub8 (fn vol)
  72   (ensure-directories-exist fn)
  73   (destructuring-bind (z y x)
  74       (array-dimensions vol)
  75     (let ((b (make-array (list y x)
  76                          :element-type '(unsigned-byte 8))))
  77       (dotimes (k z)
  78         (do-rectangle (j i 0 y 0 x)
  79           (setf (aref b j i)
  80                 (aref vol k j i)))
  81         (write-pgm b (format nil "~a/~3,'0d.pgm" fn k)))))
  82   nil)
  83
  84
  85 (declaim (ftype (function (fixnum fixnum string
  86                                   (simple-array (complex double-float) 3))
  87                           (values null &optional))
  88                 save-scaled-stack))
  89 (defun save-scaled-stack (h w fn vol)
  90     (destructuring-bind (z y x)
  91       (array-dimensions vol)
  92       (let* ((b (make-array (list y x)
  93                             :element-type '(unsigned-byte 8)))
  94              (bb (make-array (list h w)
  95                              :element-type '(unsigned-byte 8))))
  96        (dotimes (k z)
  97          (do-rectangle (j i 0 y 0 x)
  98            (setf (aref b j i)
  99                  (clamp (floor (* 255 (realpart (aref vol k j i)))))))
 100          (do-rectangle (j i 0 h 0 w)
 101            (setf (aref bb j i)
 102                  (floor (interpolate2 b (* (/ 1d0 h) j y) (* (/ 1d0 w) i x)))))
 103          (write-pgm bb (format nil "~a~3,'0d.pgm" fn k)))))
 104   nil)
 105
 106 (declaim (ftype (function ((simple-array (complex double-float) 3)
 107                            (simple-array (complex double-float) 3))
 108                           (values (simple-array (complex double-float) 3) &optional))
 109                 .*))
 110 (defun .* (vola volb)
 111   (let ((result (make-array (array-dimensions vola)
 112                             :element-type (array-element-type vola))))
 113    (destructuring-bind (z y x)
 114        (array-dimensions vola)
 115      (do-box (k j i 0 z 0 y 0 x)
 116        (setf (aref result k j i)
 117              (* (aref vola k j i)
 118                 (aref volb k j i)))))
 119    result))
 120
 121 (declaim (ftype (function (double-float (simple-array (complex double-float) 3))
 122                           (values (simple-array (complex double-float) 3) &optional))
 123                 s*))
 124
 125 (defun s* (s vol)
 126   (let* ((a (sb-ext:array-storage-vector vol))
 127          (n (length a)))
 128     (dotimes (i n)
 129       (setf (aref a i) (* s (aref a i)))))
 130   vol)
 131
 132
 133 (declaim (ftype (function ((simple-array (complex double-float) 3)
 134                            (simple-array (complex double-float) 3))
 135                           (values (simple-array (complex double-float) 3) &optional))
 136                 convolve3-circ convolve3))
 137 (defun convolve3-circ (vola volb)
 138   (let* ((da (array-dimensions vola))
 139          (db (array-dimensions volb))
 140          (compare-ab (map 'list #'(lambda (x y) (eq x y)) da db)))
 141     (when (some #'null compare-ab)
 142       (error "convolve3-circ expects both input arrays to have the same dimensions.")))
 143   (ift3 (.* (ft3 vola) (ft3 volb))))
 144
 145 ;; volb is the kernel
 146 (defun convolve3 (vola volb)
 147   (destructuring-bind (za ya xa)
 148       (array-dimensions vola)
 149     (destructuring-bind (zb yb xb)
 150         (array-dimensions volb)
 151       (let* ((biga (make-array (list (+ za zb)
 152                                      (+ ya yb)
 153                                      (+ xa xb))
 154                                :element-type '(complex double-float)))
 155              (bigb (make-array (array-dimensions biga)
 156                                :element-type '(complex double-float)))
 157              (fzb (floor zb 2))
 158              (fyb (floor yb 2))
 159              (fxb (floor xb 2))
 160              (fza (floor za 2))
 161              (fya (floor ya 2))
 162              (fxa (floor xa 2)))
 163         (do-box (k j i 0 za 0 ya 0 xa)
 164           (setf (aref biga (+ k fzb) (+ j fyb) (+ i fxb))
 165                 (aref vola k j i)))
 166         (do-box (k j i 0 zb 0 yb 0 xb)
 167           (setf (aref bigb (+ k fza) (+ j fya) (+ i fxa))
 168                 (aref volb k j i)))
 169         (let* ((conv (convolve3-circ biga (fftshift3 bigb)))
 170                (result (make-array (array-dimensions vola)
 171                                    :element-type '(complex double-float))))
 172           (do-box (k j i 0 za 0 ya 0 xa)
 173             (setf (aref result k j i)
 174                   (aref conv (+ k fzb) (+ j fyb) (+ i fxb))))
 175           result)))))
 176
 177 #+nil
 178 (let ((a (make-array (list 100 200 300)
 179                      :element-type '(complex double-float)))
 180       (b (make-array (list 10 200 30)
 181                      :element-type '(complex double-float))))
 182   (convolve3 a b)
 183   nil)
 184
 185 (declaim (ftype (function (double-float fixnum fixnum fixnum &key (:scale-z double-float))
 186                           (values (simple-array (complex double-float) 3) &optional))
 187                 draw-sphere))
 188 (defun draw-sphere (radius z y x &key (scale-z 1d0))
 189   (let ((sphere (make-array (list z y x)
 190                             :element-type '(complex double-float))))
 191     (do-box (k j i 0 z 0 y 0 x)
 192       (let ((r (sqrt (+ (square (- i (* .5d0 x)))
 193                         (square (- j (* .5d0 y)))
 194                         (square (* scale-z (- k (* .5d0 z))))))))
 195         (setf (aref sphere k j i)
 196               (if (< r radius)
 197                   (complex 1d0)
 198                   (complex 0d0)))))
 199     sphere))
 200
 201 (declaim (ftype (function ()
 202                           (values (simple-array vec-i 1)
 203                                   (simple-array double-float 1)
 204                                   cons
 205                                   &optional))
 206                 find-centers))
 207 (defvar *stack* ())
 208 ;; to find centers of cells do a convolution with a sphere
 209 (defun find-centers ()
 210   (let* ((stack-byte (read-stack "/home/martin/tmp/xa*.pgm"))
 211          (stack (make-array (array-dimensions stack-byte)
 212                             :element-type '(complex double-float))))
 213     (destructuring-bind (z y x)
 214         (array-dimensions stack)
 215       (dotimes (k z)
 216         (dotimes (j y)
 217           (dotimes (i x)
 218             (let ((v (+ (* .43745d0 k)
 219                         (aref stack-byte k j i))))
 220               (setf
 221                (aref stack-byte k j i) (clamp (floor v))
 222                (aref stack k j i) (complex v))))))
 223       #+nil(setf *stack* stack)
 224       ;; find centers of cells by convolving with sphere
 225       (let* ((sphere (draw-sphere 11d0 z y x :scale-z 5d0))
 226              (conv (convolve3-circ stack (fftshift3 sphere)))
 227              (conv-byte (make-array (list y x)
 228                                    :element-type '(unsigned-byte 8)))
 229              (centers (make-array 0 :element-type 'vec-i
 230                                   :initial-element (make-vec-i)
 231                                   :adjustable t
 232                                   :fill-pointer t))
 233              (center-heights (make-array 0 :element-type 'double-float
 234                                   :adjustable t
 235                                   :fill-pointer t)))
 236         (loop for k from 6 below (- z 3) do
 237              (loop for j from 1 below (1- y) do
 238                   (loop for i from 1 below (1- x) do
 239                        (let ((v (abs (aref conv k j i))))
 240                          (setf (aref conv-byte j i)
 241                                (if (and (< (abs (aref conv k j (1- i))) v)
 242                                         (< (abs (aref conv k j (1+ i))) v)
 243                                         (< (abs (aref conv k (1- j) i)) v)
 244                                         (< (abs (aref conv k (1+ j) i)) v)
 245                                         (< (abs (aref conv (1- k) j i)) v)
 246                                         (< (abs (aref conv (1+ k) j i)) v))
 247                                    (progn
 248                                      (vector-push-extend
 249                                       (make-vec-i :z k :y j :x i)
 250                                       centers)
 251                                      (vector-push-extend v center-heights)
 252                                      0)
 253                                    (clamp (floor (/ v (* 4e2 x y z)))))))))
 254              #+nil(write-pgm conv-byte
 255                         (format nil "/home/martin/tmp/conv~3,'0d.pgm" k)))
 256         (let ((c (make-array (length centers)
 257                              :element-type 'vec-i
 258                              :initial-contents centers))
 259               (ch (make-array (length center-heights)
 260                              :element-type 'double-float
 261                              :initial-contents center-heights)))
 262           (values c ch (array-dimensions stack)))))))
 263
 264
 265 (declaim (ftype (function (fixnum)
 266                           (values fixnum &optional))
 267                 ensure-even))
 268 (defun ensure-even (x)
 269   (if (eq 1 (mod x 2))
 270       (1+ x)
 271       x))
 272
 273 (declaim (ftype (function (double-float (simple-array vec-i 1)
 274                                         &key (:div-x double-float)
 275                                         (:div-y double-float)
 276                                         (:div-z double-float))
 277                           (values (simple-array (complex double-float) 3)
 278                                   &optional))
 279                 draw-scaled-spheres))
 280 ;; put points into the centers of nuclei and convolve a sphere around each
 281 (defun draw-scaled-spheres (radius centers &key (div-x 5d0)
 282                      (div-y 5d0)
 283                      (div-z 1d0))
 284   (let* ((max-x (floor (reduce #'max
 285                                (map '(simple-array fixnum 1) #'vec-i-x
 286                                     centers)) div-x))
 287          (max-y (floor (reduce #'max
 288                                (map '(simple-array fixnum 1) #'vec-i-y
 289                                     centers)) div-y))
 290          (max-z (floor (reduce #'max
 291                                (map '(simple-array fixnum 1) #'vec-i-z
 292                                     centers)) div-z))
 293          (cr (ceiling radius))
 294          (rh (floor radius 2))
 295          (x (ensure-even (+ max-x cr)))
 296          (y (ensure-even (+ max-y cr)))
 297          (z (ensure-even (+ max-z cr)))
 298          (dims (list z y x))
 299          (points (make-array dims
 300                              :element-type '(complex double-float))))
 301     (loop for c across centers do
 302          (setf (aref points
 303                      (- (floor (vec-i-z c) div-z) rh)
 304                      (- (floor (vec-i-y c) div-y) rh)
 305                      (- (floor (vec-i-x c) div-x) rh))
 306                (complex 1d0 0d0)))
 307     (convolve3-circ points (draw-sphere radius z y x))))
 308
 309 (declaim (ftype (function ((simple-array (complex double-float) 3))
 310                           (values (simple-array (unsigned-byte 8) 3)
 311                                   &optional))
 312                 convert-vol))
 313 (defun convert-vol (vol)
 314   (destructuring-bind (z y x)
 315       (array-dimensions vol)
 316    (let ((result (make-array (array-dimensions vol)
 317                              :element-type '(unsigned-byte 8))))
 318      (do-box (k j i 0 z 0 y 0 x)
 319        (setf (aref result k j i)
 320              (clamp (floor (* 255 (abs (aref vol k j i)))))))
 321      result)))
 322
 323 (declaim (ftype (function ((simple-array (complex double-float) 2)
 324                            &optional function)
 325                           (values (simple-array (unsigned-byte 8) 2)
 326                                   &optional))
 327                 convert-img))
 328 (defun convert-img (img &optional (function #'abs))
 329   (destructuring-bind (y x)
 330       (array-dimensions img)
 331    (let ((result (make-array (array-dimensions img)
 332                              :element-type '(unsigned-byte 8))))
 333      (do-rectangle (j i 0 y 0 x)
 334        (setf (aref result j i)
 335              (clamp (floor (funcall function (aref img j i))))))
 336      result)))
 337
 338 (declaim (ftype (function (double-float (simple-array vec-i 1)
 339                                         fixnum fixnum
 340                                         fixnum)
 341                           (values (simple-array (complex double-float) 3)
 342                                   &optional))
 343                 draw-spheres))
 344 ;; put points into the centers of nuclei and convolve a sphere around each
 345 (defun draw-spheres (radius centers z y x)
 346   (let* ((dims (list z y x))
 347          (points (make-array dims
 348                              :element-type '(complex double-float)))
 349          (n (length centers)))
 350     (dotimes (i n)
 351       (let ((c (aref centers i)))
 352        (setf (aref points
 353                    (* 5 (vec-i-z c))
 354                    (vec-i-y c)
 355                    (vec-i-x c))
 356              (complex 1d0 0d0))))
 357     (convolve3-circ points (fftshift3 (draw-sphere radius z y x)))))
 358
 359
 360
 361 #+nil
 362 (progn
 363  (defparameter *merge*
 364    (let ((a (make-array (array-dimensions *stack*)
 365                         :element-type '(unsigned-byte 8))))
 366      (destructuring-bind (z y x)
 367          (array-dimensions *stack*)
 368        (do-box (k j i 0 z 0 y 0 x)
 369          (setf (aref a k j i)
 370                (clamp (if (eq 0 (aref *blobs*
 371                                       k j i))
 372                           (* (aref *stack* k j i) 2)
 373                           0))))
 374        a)))
 375  (save-stack-ub "/home/martin/tmp/merge" *merge*))
 376
 377 ;; (- z k 1) (- y j 1) (- x i 1)
 378
 379
 380
 381 #+nil ;; model with isotropic pixels
 382 (save-stack-ub "/home/martin/tmp/iso/iso"
 383                (convert-vol (draw-scaled-spheres 2d0 *centers*)))
 384
 385 #+nil
 386 (save-stack-ub "/home/martin/tmp/blobs" *blobs*)
 387
 388 #+nil
 389 (save-scaled-stack 100 100 "/home/martin/tmp/sca" *blobs*)
 390
 391 #+nil ;; find maximum
 392 (reduce #'max (map 'vector #'abs (sb-ext:array-storage-vector *stack*)))
 393
 394 #+nil ;; scale to 1
 395 (destructuring-bind (z y x)
 396     (array-dimensions *stack*)
 397     (do-box (k j i 0 z 0 y 0 x)
 398       (setf (aref *stack* k j i)
 399             (/ (aref *stack* k j i) 257d0))))
 400
 401 #+nil
 402 (save-stack "/home/martin/tmp/stack" *stack*)
 403
 404
 405 #+nil ;; make a text image of the psf
 406 (let ((numerical-aperture 1.38d0))
 407   (psf:init :numerical-aperture numerical-aperture)
 408   (multiple-value-bind (u v)
 409       (psf:get-uv 0 1.5 3 :numerical-aperture numerical-aperture)
 410    (let* ((nu 31)
 411           (nv 61)
 412           (a (psf:integ-all nu nv  (/ u nu) (/ v nv)))
 413           (max (aref a 0 0))
 414           (scale (/ 99 max)))
 415      (destructuring-bind (uu vv)
 416          (array-dimensions a)
 417        (format t "~%")
 418        (dotimes (u uu)
 419          (dotimes (v vv)
 420            (format t "~2d" (truncate (* scale (aref a u v)))))
 421          (format t "~%"))))))
 422
 423
 424 #+nil
 425 (time
 426  (let ((z 18d0))
 427   (save-stack "/home/martin/tmp/psf"
 428               (fftshift3 (ft3 (sim-psf 128 128 128 z (* .5d0 z))))
 429               :function #'(lambda (x) (* .01d0 (abs x))))))
 430
 431 ;; worm stack 0.198 um in X and Y and 1 um in Z
 432
 433 #+nil
 434 (array-dimensions *blobs*)
 435
 436 #+nil ;; write ft of object
 437 (time
 438  (save-stack "/home/martin/tmp/kblobs" (fftshift3 (ft3 *blobs*))
 439              :function #'(lambda (x) (* 1d-4 (abs x)))))
 440
 441 #+nil ;; write ft of psf
 442 (time
 443  (destructuring-bind (z y x)
 444      (array-dimensions *blobs*)
 445    (save-stack "/home/martin/tmp/otf"
 446                (fftshift3 (ft3 (sim-psf z y x
 447                                         (* .198d0 z) (* .198d0 (ceiling (* (sqrt 2d0) (max y x)))))))
 448               :function #'(lambda (x) (* 1d-1 (abs x))))))
 449
 450 (declaim (ftype (function ((simple-array (complex double-float) (* * *)))
 451                           (values (simple-array (complex double-float) (* * *))
 452                                   &optional))
 453                 zshift3))
 454 (defun zshift3 (in)
 455   (let ((out (make-array (array-dimensions in)
 456                          :element-type '(complex double-float))))
 457    (destructuring-bind (w2 w1 w0)
 458        (array-dimensions in)
 459      (dotimes (k w2)
 460        (dotimes (j w1)
 461          (dotimes (i w0)
 462            (let* ((kk (if (> k (/ w2 2))
 463                           (+ w2 (/ w2 2) (- k))
 464                           (- (/ w2 2) k))))
 465              (setf (aref out k j i)
 466                    (aref in kk j i))
 467              nil)))))
 468    out))
 469
 470
 471
 472 #+nil
 473 (time
 474  (destructuring-bind (z y x)
 475      (array-dimensions *blobs*)
 476   (let ((psf (sim-psf z y x
 477                       (* .198d0 z) (* .198d0 (ceiling (* (sqrt 2d0) (max y x)))))))
 478     (save-stack "/home/martin/tmp/impsf"
 479                 (convolve3 *blobs* psf)
 480                 :function #'(lambda (x) (* 1d-8 (realpart x)))))))
 481
 482
 483 (declaim (ftype (function ((simple-array (complex double-float) 3)
 484                            &optional fixnum)
 485                           (values (simple-array (complex double-float) 2)
 486                                   &optional))
 487                 cross-section-xz))
 488 (defun cross-section-xz (a &optional (y (floor (array-dimension a 1) 2)))
 489   (destructuring-bind (z y-size x)
 490       (array-dimensions a)
 491     (unless (<= 0 y (1- y-size))
 492       (error "Y is out of bounds."))
 493     (let ((b (make-array (list z x)
 494                          :element-type `(complex double-float))))
 495       (do-rectangle (j i 0 z 0 x)
 496         (setf (aref b j i)
 497               (aref a j y i)))
 498       b)))
 499
 500 (declaim (ftype (function ((simple-array (complex double-float) 2)
 501                            &key (:function function))
 502                           (values (simple-array (unsigned-byte 8) 2)
 503                                   &optional))
 504                 normalize-img))
 505 (defun normalize-img (a &key (function #'abs))
 506   (destructuring-bind (y x)
 507       (array-dimensions a)
 508     (let ((b (make-array (list y x)
 509                          :element-type 'double-float)))
 510       (do-rectangle (j i 0 y 0 x)
 511         (setf (aref b j i)
 512               (funcall function (aref a j i))))
 513       (let* ((b1 (sb-ext:array-storage-vector b))
 514              (ma (reduce #'max b1))
 515              (mi (reduce #'min b1))
 516              (result (make-array (list y x)
 517                                  :element-type '(unsigned-byte 8)))
 518              (s (/ 255d0 (- ma mi))))
 519         (do-rectangle (j i 0 y 0 x)
 520           (setf (aref result j i)
 521                 (floor (* s (- (aref b j i) mi)))))
 522         result))))
 523
 524 (declaim (ftype (function ((simple-array (complex double-float) 3)
 525                            &key (:function function))
 526                           (values (simple-array (unsigned-byte 8) 3)
 527                                   &optional))
 528                 normalize-vol))
 529 (defun normalize-vol (a &key (function #'abs))
 530   (destructuring-bind (z y x)
 531       (array-dimensions a)
 532     (let ((b (make-array (list z y x)
 533                          :element-type 'double-float)))
 534       (do-box (k j i 0 z 0 y 0 x)
 535         (setf (aref b k j i)
 536               (funcall function (aref a k j i))))
 537       (let* ((b1 (sb-ext:array-storage-vector b))
 538              (ma (reduce #'max b1))
 539              (mi (reduce #'min b1))
 540              (result (make-array (list z y x)
 541                                  :element-type '(unsigned-byte 8)))
 542              (s (/ 255d0 (- ma mi))))
 543         (do-box (k j i 0 z 0 y 0 x)
 544           (setf (aref result k j i)
 545                 (floor (* s (- (aref b k j i) mi)))))
 546         result))))
 547
 548 #+nil ;; compare intensity and |e-field|^2
 549 (time
 550  (let ((z 128)
 551        (y 128)
 552        (x 128))
 553   (multiple-value-bind (e0 e1 e2)
 554       (psf:electric-field-psf z x y 10d0 5d0)
 555     (let ((intens (make-array (array-dimensions e0)
 556                               :element-type '(complex double-float))))
 557       (do-box (k j i 0 z 0 y 0 x)
 558         (setf (aref intens k j i) (complex (+ (psf::abs2 (aref e0 k j i))
 559                                               (psf::abs2 (aref e1 k j i))
 560                                               (psf::abs2 (aref e2 k j i))))))
 561       (let* ((k0 (fftshift3 (ft3 intens)))
 562              (k1 (fftshift3 (ft3 (psf:intensity-psf z y x 10d0 5d0))))
 563              (k- (make-array (array-dimensions k0)
 564                              :element-type '(complex double-float))))
 565         (do-box (k j i 0 z 0 y 0 x)
 566           (setf (aref k- k j i) (- (aref k0 k j i)
 567                                    (aref k1 k j i))))
 568         (save-stack "/home/martin/tmp/intens0"
 569                     k-
 570                     :function #'(lambda (x) (* 1d-1 (abs x))))
 571         (write-pgm (convert-img (cross-section-xz k-)
 572                                 #'(lambda (z) (* 1e-1 (abs z))))
 573                    "/home/martin/tmp/intens0xz.pgm"))))))
 574
 575 ;; create psf by blocking light in the back-focal plane. only a disk
 576 ;; centered on point p=(CENTER-X, CENTER-Y) with r=RADIUS is left
 577 ;; transparent. p=(0,0) and r=1 will unblock the full
 578 ;; bfp. p=(.9,0). r=.1 will illuminate with a high angle from positive
 579 ;; x.
 580
 581 ;; SIZE-X and SIZE-Z give the extend in um. the extend in the y
 582 ;; direction is set to SIZE-X.
 583 (declaim (ftype (function (double-float double-float double-float
 584                                         &key (:x fixnum) (:y fixnum) (:z fixnum)
 585                                         (:size-x double-float)
 586                                         (:size-z double-float)
 587                                         (:wavelength double-float)
 588                                         (:numerical-aperture double-float)
 589                                         (:immersion-index double-float)
 590                                         (:integrand-evaluations fixnum))
 591                           (values (simple-array (complex double-float) 3)))
 592                 angular-psf))
 593 (defun angular-psf (center-x center-y radius &key
 594                     (x 64) (y 64) (z 64) (size-x 12d0)
 595                     (size-z 12d0)
 596                     (wavelength .480d0)
 597                     (numerical-aperture 1.38d0)
 598                     (immersion-index 1.515d0)
 599                     (integrand-evaluations 31))
 600   (let* ((dx (/ x size-x)) ;; pixels/um
 601          )
 602     ;; calculate the electric field in focus
 603     (multiple-value-bind (e0 e1 e2)
 604         (psf:electric-field-psf z y x size-z size-x
 605                                 :numerical-aperture numerical-aperture
 606                                 :immersion-index immersion-index
 607                                 :wavelength wavelength
 608                                 :integrand-evaluations integrand-evaluations)
 609       (write-pgm (normalize-img (cross-section-xz e0))
 610                  "/home/martin/tmp/cut-asf-e0-no.pgm")
 611       ;; k0, k1 and k2 contain sections of the k-sphere
 612       (let* ((k0 (fftshift3 (ft3 e0)))
 613              (k1 (fftshift3 (ft3 e1)))
 614              (k2 (fftshift3 (ft3 e2)))
 615              ;; radius of the circular border of the cut k-sphere
 616              ;; sin(phi)=r/k, k=2pi/lambda, sin(phi)=NA/n -> r=k*sin(phi)
 617              (sinphi (/ numerical-aperture immersion-index))
 618              (klen (/ (* 2d0 pi) wavelength))
 619              (bfp-radius (* klen sinphi))
 620              (bfp-radius-pixels (* bfp-radius (/ dx (* 2d0 (sqrt 2d0)))))
 621              ;; wiil contain distance to center-x,center-y
 622              (rad-a (make-array (list y x)
 623                                 :element-type 'double-float)))
 624         (defparameter *efield-k* (list k0 k1 k2))
 625         (write-pgm (normalize-img (cross-section-xz k0))
 626                    "/home/martin/tmp/cut-e0-no.pgm")
 627         ;; calculate distances to center point for all points in one slice
 628         (do-rectangle (j i 0 y 0 x)
 629           (let* ((ii (- i (floor x 2) (* center-x bfp-radius-pixels)))
 630                  (jj (- j (floor y 2) (* center-y bfp-radius-pixels)))
 631                  (radius (sqrt (+ (* 1d0 ii ii) (* jj jj)))))
 632             (setf (aref rad-a j i) radius)))
 633
 634         ;; set fourier field to zero if outside of transparent circle
 635         (do-box (k j i 0 z 0 y 0 x)
 636           (when (< (* radius bfp-radius-pixels) (aref rad-a j i))
 637             (setf (aref k0 k j i) (complex 0d0)
 638                   (aref k1 k j i) (complex 0d0)
 639                   (aref k2 k j i) (complex 0d0))))
 640         (write-pgm (normalize-img (cross-section-xz k0))
 641                    "/home/martin/tmp/cut-e0.pgm")
 642
 643         ;; deteriorated electric fields
 644         (setf e0 (ift3 (fftshift3 k0))
 645               e1 (ift3 (fftshift3 k1))
 646               e2 (ift3 (fftshift3 k2)))
 647
 648         ;; calculate energy density of electric field
 649         (let ((density (make-array (array-dimensions e0)
 650                                    :element-type '(complex double-float))))
 651           (do-box (k j i 0 z 0 y 0 x)
 652             (setf (aref density k j i)
 653                   (complex (+ (psf:abs2 (aref e0 k j i))
 654                               (psf:abs2 (aref e1 k j i))
 655                               (psf:abs2 (aref e2 k j i))))))
 656           (write-pgm (normalize-img (cross-section-xz density))
 657                     "/home/martin/tmp/cut-intens.pgm")
 658           density)))))
 659
 660 #+nil ;; find centers of nuclei
 661 (time
 662  (multiple-value-bind (c ch dims)
 663        (find-centers)
 664    (defparameter *centers* c)
 665       (defparameter *center-heights* ch)
 666    (defparameter *dims* dims)
 667    (sb-ext:gc :full t)))
 668
 669 #+nil ;; draw the spheres
 670 (let ((spheres
 671        (destructuring-bind (z y x)
 672            *dims*
 673          (draw-spheres 7d0 *centers* (* 5 z) y x))))
 674   (setf *spheres* spheres)
 675   (save-stack-ub8 "/home/martin/tmp/spheres" (normalize-vol *spheres*)))
 676
 677
 678 #+nil ;; construct LCOS image
 679 (let ((coord (aref *centers* 0))
 680       (radius 7d0)
 681       (slice (make-array (array-dimensions *spheres*)
 682                          :element-type '(complex double-float))))
 683   (destructuring-bind (z y x)
 684       (array-dimensions *spheres*)
 685     ;; draw only the center sphere
 686     (let* ((xc (vec-i-x coord))
 687            (yc (vec-i-y coord))
 688            (zc (* 5 (vec-i-z coord)))
 689            (k  zc))
 690       (do-rectangle (j i 0 y 0 x)
 691        (let ((r (sqrt (+ (square (* 1d0 (- i xc)))
 692                          (square (* 1d0 (- j yc)))
 693                          (square (* 1d0 (- k zc)))))))
 694          (setf (aref slice k j i)
 695                (if (< r radius)
 696                    (complex 255d0)
 697                    (complex 0d0)))))))
 698   (defparameter *slice* slice)
 699   (save-stack-ub8 "/home/martin/tmp/slice" (convert-vol slice))
 700   #+nil (write-pgm (normalize-img (cross-section-xz slice (* 5 (vec-i-z coord))))
 701              "/home/martin/tmp/cut-intens2.pgm"))
 702
 703 (defparameter *bfp-circ-radius* .5d0)
 704 (defparameter *bfp-circ-center-x* (- .999d0 *bfp-circ-radius*))
 705
 706 #+nil
 707 (time
 708  (progn
 709   (angular-psf *bfp-circ-center-x* 0d0 *bfp-circ-radius* :x 200 :y 200 :z 200
 710                :size-x (* 200 .2d0) :size-z (* 200 .2d0)
 711                :integrand-evaluations 160)
 712   nil))
 713
 714 #+nil ;; light distribution in the specimen
 715 ;; default resolution is isotropic 12 um /64 = 187.5 nm/pixel
 716 (time
 717  (let* ((radius *bfp-circ-radius*)
 718         (x *bfp-circ-center-x*)
 719         (xx 300)
 720         (yy xx)
 721         (zz 300)
 722         (dx .2d0)
 723         (psf (angular-psf x 0d0 radius :x xx :y yy :z zz
 724                           :size-x (* xx dx) :size-z (* zz dx)
 725                           :integrand-evaluations 190))
 726         (dims (destructuring-bind (z y x)
 727                   *dims*
 728                 (list (* z 5) y x)))
 729         #+nil (psf-big (make-array dims
 730                              :element-type '(complex double-float))))
 731    #+nil   (destructuring-bind (z y x)
 732                dims
 733              (let ((ox (- (floor x 2) (floor xx 2)))
 734                    (oy (- (floor y 2) (floor yy 2)))
 735                    (oz (- (floor z 2) (floor zz 2))))
 736                (do-box (k j i 0 zz 0 yy 0 xx)
 737                  (setf (aref psf-big (+ oz k) (+ oy j) (+ ox i))
 738                        (aref psf k j i)))))
 739    #+nil (defparameter *psf-big* psf-big)
 740    #+nil (save-stack-ub8 "/home/martin/tmp/psf-big" (normalize-vol psf-big))
 741    (defparameter *psf* psf)
 742    (sb-ext:gc :full t)
 743    (defparameter *slice-x-psf* (convolve3 *slice* psf))
 744    (sb-ext:gc :full t)))
 745
 746 #+nil
 747 (save-stack-ub8 "/home/martin/tmp/psf" (normalize-vol *psf*))
 748
 749
 750 ;;                                                   /-
 751 ;;           ---                              /---
 752 ;;              \-----                    /---
 753 ;;                    \-----          /---
 754 ;;                          \---- /---
 755 ;;                            /--------
 756 ;;                        /---         \----
 757 ;;                    /---                  \-----
 758 ;;                /---                            \-----
 759 ;;              --                                      \--
 760
 761
 762
 763
 764 #+nil ;; draw lines into the light distribution in the specimen
 765 (destructuring-bind (z y x)
 766     (array-dimensions *slice-x-psf*)
 767   (let ((coord (elt *centers* 0))
 768         (vol (normalize-vol *slice-x-psf*))
 769         (dx 2.d-4))
 770     (loop for pos in (list (list (* (- (floor x 2) (- (vec-i-x coord) 7)) dx)
 771                                  (* (- (floor y 2) (vec-i-y coord)) dx))
 772                            (list (* (- (floor x 2) (+ (vec-i-x coord) 7)) dx)
 773                                  (* (- (floor y 2) (vec-i-y coord)) dx))) do
 774          (loop for angle in (list 0d0
 775                                   -.4d0
 776                                   (- (- *bfp-circ-center-x* *bfp-circ-radius*))
 777                                   (- (+ *bfp-circ-center-x* *bfp-circ-radius*))) do
 778               (draw-ray-into-vol (first pos) (second pos)
 779                                  angle 0d0
 780                                  vol
 781                                  :shift-z (- (* 5d0 (vec-i-z coord))
 782                                              (floor z 2)))))
 783     (save-stack-ub8 "/home/martin/tmp/slice-x-psf" vol)))
 784
 785
 786 #+nil ;; draw lines into the volume of the psf
 787 (destructuring-bind (z y x)
 788     (array-dimensions *psf-big*)
 789   (let ((vol (normalize-vol *psf-big*))
 790         (dx 2.d-4))
 791     (draw-ray-into-vol 0d0 0d0
 792                        (- (- 1d0 1d-4)) 0d0
 793                        vol
 794                        :center-z (floor z 2))
 795     (draw-ray-into-vol 0d0 0d0
 796                        -.7d0 0d0 vol
 797                        :center-z (floor z 2))
 798     (draw-ray-into-vol 0d0 0d0
 799                        -.4d0 0d0 vol
 800                        :center-z (floor z 2))
 801     (draw-ray-into-vol 0d0 0d0
 802                        0d0 0d0 vol
 803                        :center-z (floor z 2))
 804     (save-stack-ub8 "/home/martin/tmp/psf-big" vol)))
 805
 806
 807 #+nil ;; excited fluorophores
 808 (progn
 809   (setf *slice-x-psf-times-spheres* (.* *spheres* *slice-x-psf*))
 810   (save-stack-ub8 "/home/martin/tmp/slice-x-psf-times-spheres"
 811                  (normalize-vol *slice-x-psf-times-spheres*)))
 812
 813
 814 #+nil ;; blur with detection psf
 815 (time
 816  (let* ((radius .5d0)
 817         (x (- 1d0 radius))
 818         (xx 80)
 819         (yy xx)
 820         (zz 128)
 821         (dx .2d0)
 822         (psf (psf:intensity-psf zz yy xx (* zz dx) (* xx dx)
 823                                 :integrand-evaluations 100))
 824         (dims (destructuring-bind (z y x)
 825                   *dims*
 826                 (list (* z 5) y x)))
 827         (psf-big (make-array dims
 828                              :element-type '(complex double-float))))
 829    (setf *psf-big* psf-big)
 830    (destructuring-bind (z y x)
 831        dims
 832      (let ((ox (- (floor x 2) (floor xx 2)))
 833            (oy (- (floor y 2) (floor yy 2)))
 834            (oz (- (floor z 2) (floor zz 2))))
 835        (do-box (k j i 0 zz 0 yy 0 xx)
 836          (setf (aref psf-big (+ oz k) (+ oy j) (+ ox i))
 837                (aref psf k j i)))))
 838    (save-stack-ub8 "/home/martin/tmp/psf-detect-big" (normalize-vol psf-big))
 839    (sb-ext:gc :full t)
 840    (defparameter *camera-volume* (convolve3-circ *slice-x-psf-times-spheres*
 841                                                  (fftshift3 psf-big)))
 842    (save-stack-ub8 "/home/martin/tmp/camera-volume"
 843                    (normalize-vol *camera-volume*))
 844    (sb-ext:gc :full t)))
 845
 846
 847
 848 #+nil ;; check convolution
 849 (time
 850  (let ((a (make-array (list 64 64 64)
 851                       :element-type '(complex double-float)))
 852        (b (psf:intensity-psf 64 64 64 20d0 20d0) #+nil (make-array (list 64 64 64)
 853                       :element-type '(complex double-float))))
 854    (setf (aref a 12 12 12) (complex 255d0))
 855 #+nil   (setf (aref b 0 0 0) (complex 255d0))
 856    (save-stack-ub8 "/home/martin/tmp/conv-test" (normalize-vol (convolve3-circ a (fftshift3 b))))))
 857
 858
 859 #+nil ;; output the xz cross section centered on a sphere in the middle
 860 (let ((coord (aref *centers* 30)))
 861   (write-pgm (normalize-img (cross-section-xz *spheres* (* 5 (vec-i-z coord))))
 862              "/home/martin/tmp/cut-spheres.pgm")
 863   (write-pgm (normalize-img (cross-section-xz *psf-big*))
 864              "/home/martin/tmp/cut-psf-big.pgm")
 865   (write-pgm (normalize-img (cross-section-xz *slice-x-psf* (* 5 (vec-i-z coord))))
 866              "/home/martin/tmp/cut-slice-x-psf.pgm")
 867   (write-pgm (normalize-img (cross-section-xz (.* *spheres* *slice-x-psf*) (* 5 (vec-i-z coord))))
 868              "/home/martin/tmp/cut-exfluo.pgm"))
 869
 870 #+nil
 871 (save-stack-ub8 "/home/martin/tmp/spheres" (normalize-vol *spheres*))
 872
 873 #+nil
 874 (let ((sli (make-array (array-dimensions *spheres*)
 875                        :element-type '(complex double-float))))
 876   (destructuring-bind (z y x)
 877       (array-dimensions *spheres*)
 878     (do-box (k j i 0 z 0 y 0 x)
 879       (setf (aref sli k j i)
 880             (aref *spheres* k j i)))
 881     (let ((k (* 5 (vec-i-z (aref *centers* 30)))))
 882      (do-rectangle (j i 0 y 0 x)
 883        (setf (aref sli k j i)
 884              (* 10 (aref sli k j i)))))
 885     (defparameter *sli* sli))
 886   (save-stack-ub8 "/home/martin/tmp/spheres" (normalize-vol *sli*)))
 887
 888
 889
 890 #+nil
 891 (let ((a (sb-ext:array-storage-vector *slice*)))
 892   (reduce #'max (map 'vector #'abs a)))
 893
 894 #+nil
 895 (sb-ext:gc :full t)
 896
 897 #+nil ;; model with unscaled spheres
 898 (defparameter *blobs*
 899   (destructuring-bind (z y x)
 900       *dims*
 901     (draw-spheres 7d0 *centers* (* 5 z) y x)))
 902
 903
 904 #+nil ;; print ft of angular psf
 905 (time (let* ((radius .5d0)
 906         (x (- 1d0 radius))
 907         (psf (angular-psf x 0d0 radius)))
 908    (write-pgm (normalize-img (cross-section-xz (fftshift3 (ft3 psf))))
 909               "/home/martin/tmp/cut-intens.pgm")))
 910 #+nil
 911 (let* ((intens (psf:intensity-psf 64 64 64 10d0 5d0))
 912        (k0 intens #+nil(fftshift3 (ft3 intens))))
 913   (save-stack "/home/martin/tmp/intens1" k0
 914               :function #'(lambda (x) (* 1d-5 (abs x))))
 915   (write-pgm (convert-img
 916               (cross-section-xz k0)
 917               #'(lambda (z) (* 1d-4 (abs z))))
 918              "/home/martin/tmp/intens1xz.pgm"))
 919
 920 #+nil
 921 (time
 922  (destructuring-bind (z y x)
 923      (array-dimensions *blobs*)
 924    (save-stack "/home/martin/tmp/blobs"
 925                *blobs*
 926                )))
 927
 928
 929 #+nil ;; clean up the garbage
 930 (sb-ext:gc :full t)
 931
 932
 933 #+nil
 934 (sb-vm:memory-usage :print-spaces t :count-spaces t)
 935 #+nil
 936 (sb-vm:memory-usage)
 937
 938 #+nil
 939 (sb-vm:instance-usage :dynamic)
 940
 941 #+nil
 942 (sb-vm:list-allocated-objects :dynamic)
 943
 944 #+nil
 945 (sb-vm:print-allocated-objects :dynamic)
 946
 947 #+nil
 948 (format t "~a~%" (sb-vm::type-breakdown :dynamic))
 949
 950
 951 (declaim (ftype (function (double-float)
 952                           (values double-float &optional))
 953                 sq))
 954 (defun sq (x)
 955   (* x x))
 956
 957
 958 (declaim (ftype (function ((array double-float *))
 959                           (values double-float &optional))
 960                 rosenbrock))
 961 (defun rosenbrock (p)
 962   (let* ((x (aref p 0))
 963          (y (aref p 1))
 964          (result (+ (sq (- 1 x))
 965                     (* 100 (sq (- y (sq x)))))))
 966     (format t "~a~%" (list 'rosenbrock p result))
 967     result))
 968 #+nil
 969 (rosenbrock (make-array 2 :element-type 'double-float
 970                          :initial-contents (list 1.5d0 1.5d0)))
 971 ;; run the following code to test the downhill simplex optimizer on a
 972 ;; 2d function:
 973
 974 ;; +-----        |       |
 975 ;; |     \--     |       |
 976 ;; |        \-   |       |
 977 ;; |          \- |       |
 978 ;; |            \|       |
 979 ;; |            ||       |
 980 ;; |             |       |
 981 ;; |-------------+-------+------- <- z
 982 ;; -nf          /0       f
 983 ;; object      lens     bfp
 984
 985 #+nil
 986 (time (let ((start (make-array 2 :element-type 'double-float
 987                                :initial-contents (list 1.5d0 1.5d0))))
 988         (simplex-anneal:anneal (simplex-anneal:make-simplex start 1d0)
 989                 #'rosenbrock
 990                 :ftol 1d-5)))
 991
 992 (defun draw-ray-into-vol (x-mm y-mm bfp-ratio-x bfp-ratio-y vol
 993                           &key (dx-mm .2d-3)
 994                           (shift-z 0d0))
 995   (destructuring-bind (z y x)
 996       (array-dimensions vol)
 997    (let* ((f (lens:focal-length-from-magnification 63d0))
 998           (na 1.38d0)
 999           (ri 1.515d0)
1000           (bfp-radius (lens:back-focal-plane-radius f na))
1001           (obj (lens:make-thin-objective :normal (v 0d0 0d0 -1d0)
1002                                          :center (v)
1003                                          :focal-length f
1004                                          :radius bfp-radius
1005                                          :numerical-aperture na
1006                                          :immersion-index ri))
1007           (theta (lens:find-inverse-ray-angle x-mm y-mm obj))
1008           (phi (atan y-mm x-mm))
1009           (start (lens:v (* bfp-ratio-x bfp-radius)
1010                          (* bfp-ratio-y bfp-radius)
1011                          f))
1012           (dx dx-mm)
1013           (cz (* .5d0 z)) ;; position that is in the center of front focal plane
1014           (cy (* .5d0 y))
1015           (cx (* .5d0 x))
1016           (nf (* ri f)))
1017      (macrolet ((plane (direction position)
1018                   ;; for defining a plane that is perpendicular to an
1019                   ;; axis and crosses it at POSITION
1020                   (declare (type (member :x :y :z) direction))
1021                   (let* ((normal (ecase direction
1022                                    (:x (lens:v 1d0))
1023                                    (:y (lens:v 0d0 1d0))
1024                                    (:z (lens:v 0d0 0d0 1d0)))))
1025                     `(let* ((pos ,position)
1026                             (center (lens::v* ,normal pos))
1027                             (outer-normal (lens::normalize center)))
1028                        (declare (type double-float pos))
1029                        (lens::make-disk :normal outer-normal :center center)))))
1030        ;; define the borders of the viewing volume, distances in mm
1031        (let ((p+z (plane :z (- (* dx (- z cz))
1032                                nf)))
1033              (p-z (plane :z (- (* dx (- (- z cz)))
1034                                nf)))
1035              (p+y (plane :y (* dx (- y cy))))
1036              (p-y (plane :y (* dx (- (- y cy)))))
1037              (p+x (plane :x (* dx (- x cx))))
1038              (p-x (plane :x (* dx (- (- x cx))))))
1039          (multiple-value-bind (ro s)
1040              (lens:thin-objective-ray obj
1041                                       start
1042                                       (lens::v* (lens:v (* (cos phi) (sin theta))
1043                                                         (* (sin phi) (sin theta))
1044                                                         (cos theta))
1045                                                 -1d0))
1046            (setf s (lens::v+ s (lens:v 0d0 0d0 (* dx shift-z))))
1047            (let* ((nro (lens::normalize ro)))
1048              (macrolet ((hit (plane)
1049                           ;; (declare (type lens::disk plane))
1050                           ;; find intersection between plane and the ray
1051                           `(multiple-value-bind (dir hit-point)
1052                                (lens::plane-ray ,plane
1053                                                 ;; shift start of vector a bit
1054                                                 s
1055                                                 nro)
1056                              (declare (ignore dir))
1057                              hit-point))
1058                         (pixel (hit-expr)
1059                           ;; convert coordinates from mm into integer pixel positions
1060                           `(let ((h ,hit-expr))
1061                              (declare (type (or null lens::vec) h))
1062                              (when h
1063                                (make-vec-i
1064                                 :z (floor (+ cz (/ (+ (aref h 2) nf) dx)))
1065                                 :y (floor (+ cy (/ (aref h 1) dx)))
1066                                 :x (floor (+ cx (/ (aref h 0) dx))))))))
1067                (let* ((h+z (pixel (hit p+z)))
1068                       (h-z (pixel (hit p-z)))
1069                       (h+y (pixel (hit p+y)))
1070                       (h-y (pixel (hit p-y)))
1071                       (h+x (pixel (hit p+x)))
1072                       (h-x (pixel (hit p-x)))
1073                       ;; make a list of all the points
1074                       (hlist (list h+z h-z h+y h-y h+x h-x))
1075                       ;; throw away points that are nil or that contain
1076                       ;; coordinates outside of the array dimensions
1077                       (filtered-hlist (remove-if-not #'(lambda (v)
1078                                                          (if v
1079                                                              (and (< -1 (vec-i-x v) x)
1080                                                                   (< -1 (vec-i-y v) y)
1081                                                                   (< -1 (vec-i-z v) z))
1082                                                              nil)) hlist))
1083                       ;; sort best points by x
1084                       (choice (sort filtered-hlist #'< :key (lambda (v) (vec-i-x v)))))
1085                  (format t "~a~%" (list 'choice choice))
1086                  (scan-convert-line3
1087                   (first choice)
1088                   (second choice)
1089                   vol))))))))))
1090
1091 #+nil
1092 (let ((vol (make-array (list 128 128 128) :element-type '(unsigned-byte 8))))
1093   (loop for i in '(4.0d-3 -.2d-3) do
1094    (draw-ray-into-vol i 0d0 .99d0 .0d0 vol)
1095    (draw-ray-into-vol i 0d0 -.99d0 .0d0 vol)
1096    (draw-ray-into-vol i 0d0 0d0 .99d0 vol)
1097    (draw-ray-into-vol i 0d0 0d0 -.99d0 vol))
1098
1099   (save-stack-ub8 "/home/martin/tmp/line"
1100                   vol))
1101
1102
1103 #+nil
1104 (let ((vol (make-array (list 128 128 128) :element-type '(unsigned-byte 8))))
1105  (draw-line3 (make-vec-i :x 108 :y 112 :z  103)
1106             (make-vec-i :x 82 :y 102 :z 10)
1107             vol))
1108
1109 (declaim (ftype (function (fixnum fixnum fixnum
1110                            fixnum fixnum fixnum
1111                            (simple-array (unsigned-byte 8) 3))
1112                           (values (simple-array (unsigned-byte 8) 3) &optional))
1113                 scan-convert-line3x scan-convert-line3y scan-convert-line3z))
1114 ;; 1986 kaufman
1115 (defun scan-convert-line3x (x1 y1 z1 x2 y2 z2 vol)
1116   ;; x2 - x1 has to be the biggest difference between endpoint
1117   ;; coordinates
1118   (let* ((x x1)
1119          (delx (- x2 x1))
1120          ;; initialization for y
1121          (y y1)
1122          (ddy (- y2 y1))
1123          (dely (abs ddy))
1124          (ysign (signum ddy))
1125          (dy (- (* 2 dely) delx)) ;; decision variable along y
1126          (yinc1 (* 2 dely)) ;; increment along y for dy<0
1127          (yinc2 (* 2 (- dely delx))) ;; inc along y for dy>=0
1128          ;; same initialization for z
1129          (z z1)
1130          (ddz (- z2 z1))
1131          (delz (abs ddz))
1132          (zsign (signum ddz))
1133          (dz (- (* 2 delz) delx))
1134          (zinc1 (* 2 delz))
1135          (zinc2 (* 2 (- delz delx))))
1136     (when (<= delx 0)
1137         (error "x2 is <= x1."))
1138     (setf (aref vol z y x) 255)
1139     (loop while (< x x2) do
1140          (incf x) ;; step in x
1141          (if (< dy 0) ;; then no change in y
1142              (incf dy yinc1) ;; update dy
1143              (progn
1144                (incf dy yinc2) ;; update dy, and
1145                (incf y ysign)));; increment/decrement y
1146
1147          (if (< dz 0)
1148              (incf dz zinc1)
1149              (progn
1150                (incf dz zinc2)
1151                (incf z zsign)))
1152          (setf (aref vol z y x) 255)))
1153   vol)
1154 ;; start from scan-convert-line3x and replace x->$, y->x, $->y
1155 (defun scan-convert-line3y (x1 y1 z1 x2 y2 z2 vol)
1156   (let* ((y y1)
1157          (dely (- y2 y1))
1158          (x x1)
1159          (ddx (- x2 x1))
1160          (delx (abs ddx))
1161          (xsign (signum ddx))
1162          (dx (- (* 2 delx) dely))
1163          (xinc1 (* 2 delx))
1164          (xinc2 (* 2 (- delx dely)))
1165          (z z1)
1166          (ddz (- z2 z1))
1167          (delz (abs ddz))
1168          (zsign (signum ddz))
1169          (dz (- (* 2 delz) dely))
1170          (zinc1 (* 2 delz))
1171          (zinc2 (* 2 (- delz dely))))
1172     (when (<= dely 0)
1173         (error "y2 is <= y1."))
1174     (setf (aref vol z y x) 255)
1175     (loop while (< y y2) do
1176          (incf y)
1177          (if (< dx 0)
1178              (incf dx xinc1)
1179              (progn
1180                (incf dx xinc2)
1181                (incf x xsign)))
1182          (if (< dz 0)
1183              (incf dz zinc1)
1184              (progn
1185                (incf dz zinc2)
1186                (incf z zsign)))
1187          (setf (aref vol z y x) 255)))
1188   vol)
1189 ;; replace x->$, z->x, $->z
1190 (defun scan-convert-line3z (x1 y1 z1 x2 y2 z2 vol)
1191   (let* ((z z1)
1192          (delz (- z2 z1))
1193          (y y1)
1194          (ddy (- y2 y1))
1195          (dely (abs ddy))
1196          (ysign (signum ddy))
1197          (dy (- (* 2 dely) delz))
1198          (yinc1 (* 2 dely))
1199          (yinc2 (* 2 (- dely delz)))
1200          (x x1)
1201          (ddx (- x2 x1))
1202          (delx (abs ddx))
1203          (xsign (signum ddx))
1204          (dx (- (* 2 delx) delz))
1205          (xinc1 (* 2 delx))
1206          (xinc2 (* 2 (- delx delz))))
1207     (when (<= delz 0)
1208         (error "z2 is <= z1."))
1209     (setf (aref vol z y x) 255)
1210     (loop while (< z z2) do
1211          (incf z)
1212          (if (< dy 0)
1213              (incf dy yinc1)
1214              (progn
1215                (incf dy yinc2)
1216                (incf y ysign)))
1217
1218          (if (< dx 0)
1219              (incf dx xinc1)
1220              (progn
1221                (incf dx xinc2)
1222                (incf x xsign)))
1223          (setf (aref vol z y x) 255)))
1224   vol)
1225
1226 (declaim (ftype (function (vec-i vec-i (simple-array (unsigned-byte 8) 3))
1227                           (values (simple-array (unsigned-byte 8) 3) &optional))
1228                 scan-convert-line3))
1229 (defun scan-convert-line3 (start end vol)
1230   (let* ((diff (v--i end start))
1231          (ls (list (list (vec-i-x diff) 2)
1232                    (list (vec-i-y diff) 1)
1233                    (list (vec-i-z diff) 0)))
1234          (diffa (mapcar #'(lambda (e) (list (abs (first e))
1235                                        (second e))) ls))
1236          ;; find the direction with the biggest difference
1237          (sorted-diff-a (sort diffa #'> :key #'car))
1238          (main-direction (second (first sorted-diff-a))) ;; 2 corresponds to x, 1->y, 0->z
1239          ;; find the order in which to deliver the points
1240          (main-diff (aref diff main-direction))
1241          ;; we have to swap the points when main-diff is negative
1242          (swap-points? (< main-diff 0))
1243          ;; create the function name to dispatch to
1244          (function (intern (format nil "SCAN-CONVERT-LINE3~a" (ecase main-direction
1245                                                                 (2 'X)
1246                                                                 (1 'Y)
1247                                                                 (0 'Z))))))
1248     (when (eq 0 main-diff)
1249       (error "start and end point are the same."))
1250     (if swap-points?
1251         (funcall function
1252                  (vec-i-x end)
1253                  (vec-i-y end)
1254                  (vec-i-z end)
1255                  (vec-i-x start)
1256                  (vec-i-y start)
1257                  (vec-i-z start)
1258                  vol)
1259         (funcall function
1260                  (vec-i-x start)
1261                  (vec-i-y start)
1262                  (vec-i-z start)
1263                  (vec-i-x end)
1264                  (vec-i-y end)
1265                  (vec-i-z end)
1266                  vol))))
1267
1268
1269 #+nil
1270 (time
1271  (let ((vol (make-array (list 128 128 128) :element-type '(unsigned-byte 8))))
1272    (save-stack-ub8 "/home/martin/tmp/line"
1273                    (scan-convert-line3 (make-vec-i :x 0 :y 0 :z 0)
1274                                        (make-vec-i :x 120 :y 127 :z 127)
1275                                        vol))))
1276
1277 ;;               |
1278 ;;        -------+-------
1279 ;;     -/  h (3) |       \---   (2) q_R=NA/ri*q_max
1280 ;;    -----------+------------/------------
1281 ;;               | alpha  /---  \-
1282 ;;               |     /--        \
1283 ;;               | /---            \
1284 ;;      ---------+-----------------+-------
1285 ;;               | (0)             / (1) q_max=1/(2*pixel)
1286 ;;
1287 ;; The resolution of the image has to be big enough to fit the top
1288 ;; section of the k-sphere with radius |k|=2pi*q_max into the k space.
1289 ;; q_max (see (1)) is due to the nyquist theorem and corresponds to 1
1290 ;; over two sample widths. The radius of the backfocal plane
1291 ;; corresponds to q_R (see (2) ri is the refractive index,
1292 ;; e.g. 1.515). It is bigger for an objective with a high NA.
1293
1294 ;; A transform with uneven number of samples doesn't have a bin for
1295 ;; the nyquist sampling (draw the unit circle and divide it into n
1296 ;; equal bins starting from e^(i*0). For uneven n there will be no bin
1297 ;; on e^-i (i.e. -1, 1, -1 ...), e.g. n=3).  For even n there will be
1298 ;; n/2+1 bins ontop of the real axis (e.g. 0=1, 1=e^(-i pi/2), 2=e^-i
1299 ;; for n=4, the arguments to the exponential are (i 2 pi j/n) for the
1300 ;; j-th bin) and n/2-1 bins below (e.g 3=e^(i pi/2)).  In order to
1301 ;; simplify fftshift3 I only consider transforms with even n.
1302 ;; fftshift moves the n/2+1 bins from the front of the array to the
1303 ;; back (for n=4: [0 1 2 3] -> [3 0 1 2]).  In the shifted array the
1304 ;; highest reverse frequency (bin 3) is mapped to index 0.  The origin
1305 ;; of k-space (see (0) in the sketch) is therefor mapped to bin n/2-1
1306 ;; (bin 1 for n=4). The nyquist frequency is in the last bin n-1 (bin
1307 ;; 3 for n=4).
1308
1309 ;; We now search for the right z-sampling dz to fit the top of the
1310 ;; sphere below the nyquist bin (which corresponds to q_max=1/(2*dz)).
1311 ;; |k|=2 pi/lambda = 2 pi q_max, with wavelength lambda
1312 ;; lambda=2 dz -> dz = lambda/2.
1313
1314 ;; We could use the same sampling x and y to represent the electric
1315 ;; field. For small numerical apertures the sampling distance can be
1316 ;; increased. This time the radius q_R has to be smaller than the nyquist
1317 ;; frequency:
1318 ;; 1/(2*dx)=q_R=NA/ri * 1/lambda
1319 ;; -> dx= lambda/2 * ri/NA= dz *ri/NA=dz*1.515/1.38
1320
1321 ;; The sampling distances dz and dx that I derived above are only good
1322 ;; to represent the amplitude psf. When the intensity is to be
1323 ;; calculated the sampling distance has to be changed to accomodate
1324 ;; for the convolution in k space.
1325
1326 ;; The height of the sphere cap (h see (3) in sketch) is
1327 ;; h=q_max-q_max*cos(alpha)=q_max ( 1-cos(alpha))
1328 ;; =q_max*(1-sqrt(1-sin(alpha)^2))=q_max*(1-sqrt(1-(NA/ri)^2)) The z
1329 ;; sample distance dz2 for the intensity psf should correspond to 1/(2
1330 ;; dz2)=2 h, i.e. dz2=1/h=dz*2/(1-sqrt(1-(NA/ri)^2))>dz so the necessary z
1331 ;; sampling distance for the intensity is in general bigger than for
1332 ;; the amplitude.
1333
1334 ;; The radius of the convolved donut shape is 2 q_R. Therefor the
1335 ;; transversal sampling distance for the intensity has to be smaller:
1336 ;; dx2=dx/2.
1337
1338 ;; As we are only interested in the intensity psf we can sample the
1339 ;; amplitude psf with a sampling distance dz2. The sphere cap is
1340 ;; possibly wrapped along the k_z direction. The transversal direction
1341 ;; of the amplitude psf has to be oversampled with dx2.
1342
1343 ;; To get an angular illumination psf we multiply the values on the
1344 ;; sphere with a k_z plane containing a disk that is centered at any
1345 ;; k_x and k_y inside the back focal plane.  Later I might want to
1346 ;; replace this with a gaussian or a more elaborate window function.
1347
1348 ;; With a sampling dx2 the radius of the backfocal plane fills half of
1349 ;; the k space. The coordinate calculations below are corrected for
1350 ;; this. So setting cx to 1. and cy to 0. below would center the
1351 ;; circle on the border of the bfp.
1352
1353 #+nil
1354 (time
1355  ;; changing z,y,x without touching dx or dz leaves the area that is
1356  ;; shown in k space constant
1357  (let* ((z 64)
1358         (y 64)
1359         (x y)
1360         (na 1.38d0)
1361         (ri 1.515d0)
1362         (lambd .480d0)
1363         (dz (* .5d0 lambd))
1364         (dz2 (* dz (/ 2d0 (- 1d0 (sqrt (- 1d0
1365                                           (let ((sinphi (/ na ri)))
1366                                             (* sinphi sinphi))))))))
1367         (dx (* dz (/ ri na)))
1368         (dx2 (* .5 dx)))
1369    (format t "~a~%" (list 'aspect dx2 dz2 (/ dx2 dz2) (/ dz2 dx2)))
1370    (multiple-value-bind (e0 e1 e2)
1371        (psf:electric-field-psf z y x (* z dz2) (* x dx2)
1372                                :integrand-evaluations 140)
1373      (write-pgm (normalize-img (cross-section-xz e0))
1374                 "/home/martin/tmp/cut-psf.pgm")
1375      (let ((k0 (fftshift3 (ft3 e0)))
1376            (k1 (fftshift3 (ft3 e1)))
1377            (k2 (fftshift3 (ft3 e2))))
1378        (write-pgm (normalize-img (cross-section-xz k0))
1379                   "/home/martin/tmp/cut-psf-k.pgm")
1380        (let* ((cr .3d0)
1381               (cx (- .5d0 cr))
1382               (cy .0d0)
1383               (cr2 (* cr cr))
1384               (mul0 (make-array (array-dimensions k0)
1385                                  :element-type '(complex double-float)))
1386               (mul1 (make-array (array-dimensions k0)
1387                                  :element-type '(complex double-float)))
1388               (mul2 (make-array (array-dimensions k0)
1389                                 :element-type '(complex double-float)))
1390               (mul (make-array (list y x)
1391                               :element-type 'double-float)))
1392          (do-rectangle (j i 0 y 0 x)
1393            (let* ((xx (- (* 4d0 (- (* i (/ 1d0 x)) .5d0)) cx))
1394                   (yy (- (* 4d0 (- (* j (/ 1d0 y)) .5d0)) cy))
1395                   (r2 (+ (* xx xx) (* yy yy))))
1396              (when (< r2 cr2)
1397                (setf (aref mul j i) 1d0))))
1398          (do-box (k j i 0 z 0 y 0 x)
1399            (setf (aref mul0 k j i) (* (aref k0 k j i) (aref mul j i))
1400                  (aref mul1 k j i) (* (aref k1 k j i) (aref mul j i))
1401                  (aref mul2 k j i) (* (aref k2 k j i) (aref mul j i))))
1402          (write-pgm (normalize-img (cross-section-xz mul0))
1403                     "/home/martin/tmp/cut-psf-k-mul.pgm")
1404          (let* ((em0 (ift3 (fftshift3 mul0)))
1405                 (em1 (ift3 (fftshift3 mul1)))
1406                 (em2 (ift3 (fftshift3 mul2)))
1407                 (intens (make-array (array-dimensions e0)
1408                                    :element-type '(complex double-float))))
1409            (do-box (k j i 0 z 0 y 0 x)
1410              (setf (aref intens k j i)
1411                    (+ (* (aref em0 k j i) (conjugate (aref em0 k j i)))
1412                       (* (aref em1 k j i) (conjugate (aref em1 k j i)))
1413                       (* (aref em2 k j i) (conjugate (aref em2 k j i))))))
1414            (write-pgm (normalize-img (cross-section-xz intens))
1415                       "/home/martin/tmp/cut-psf-intens.pgm")
1416            (let ((k (fftshift3 (ft3 intens))))
1417              (write-pgm (normalize-img (cross-section-xz k))
1418                         "/home/martin/tmp/cut-psf-intk.pgm"))))))))