run.lisp

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