run.lisp

   1 #.(require :frontend)
   2 (in-package :frontend)
   3
   4 #+nil
   5 (defparameter *model* (make-instance 'sphere-model-angular :dx .2d0 :dz 1.0d0))
   6
   7 #+nil
   8 (time
   9  (defparameter *model* (make-test-model)))
  10
  11 #+nil
  12 (write-pgm "/home/martin/tmp/model-cut.pgm"
  13            (normalize-2-csf/ub8-realpart (cross-section-xz (spheres *model*))))
  14 #+nil
  15 (times
  16  (defparameter *psf*
  17    (multiple-value-bind (conv dx dz)
  18        (angular-intensity-psf-minimal-resolution
  19         :x-um 8s0 :z-um 12s0
  20         :window-radius 1.0 :window-x 0s0 :window-y 0s0
  21         :debug t :initialize t
  22         :integrand-evaluations 100)
  23      (resample-3-csf conv dx dx dz .2 .2 1.0))))
  24 #+nil
  25 (write-pgm "/home/martin/tmp/psf-cut.pgm"
  26            (normalize-2-csf/ub8-realpart (cross-section-xz *psf*)))
  27
  28 #+nil
  29 (time (progn
  30          (defparameter *conv*
  31           (convolve (spheres *model*) *psf*))
  32         (write-pgm "/home/martin/tmp/conv-cut.pgm"
  33                    (normalize-2-csf/ub8-realpart
  34                     (cross-section-xz
  35                      #+nil (spheres *model*)
  36                      *conv*
  37                      #+nil (.- *conv* (vol::s* 1e11 (spheres *model*))))))))
  38 #+nil
  39 (write-pgm "/home/martin/tmp/conv-cut.pgm"
  40            (normalize-2-sf/ub8
  41             (cross-section-xz
  42              (.- (normalize-3-csf/sf-realpart *conv*)
  43                  (normalize-3-csf/sf-realpart (spheres *model*))))))
  44
  45 #+nil
  46 (time
  47  (progn
  48    (update-tex
  49     #+nil (normalize-3-csf/ub8-realpart *conv*)
  50     (normalize-2-csf/ub8-realpart
  51      (cross-section-xz *conv* (vec-i-y (elt (centers *model*) 0)))))
  52    nil))
  53
  54
  55 #+nil
  56 (update-scale 300 40)
  57 #+nil
  58 (update-scale 1 40)
  59
  60
  61
  62 #+nil
  63 (let* ((obj (lens:make-objective))
  64        (n (lens::immersion-index obj))
  65        (f (lens::focal-length obj))
  66        (nucleus-position (elt (centers-mm *model*) 0))
  67        (center (make-vec (vec-x nucleus-position)
  68                          (vec-y nucleus-position))))
  69   (update-view-center nucleus-position)
  70   (update-scale 2 20))
  71
  72 #+nil
  73 (defun draw-all ()
  74   (draw *model*
  75         :nucleus 0
  76         :win-x/r -0.2402d0 :win-y/r 0.0053d0
  77         ; :win-x/r -.0927d0 :win-y/r -0.2359d0
  78         :win-r/r .08d0
  79         :nr-ffp 2
  80         :nr-bfp 3
  81         :nr-theta 1))
  82
  83 #+nil
  84 (with-gui
  85   (draw-all))
  86
  87 #+nil
  88 (time
  89  (let* ((n 100)
  90         (a (make-array (list n n) :element-type '(unsigned-byte 8)))
  91         (nn (length *spheres-c-r*))
  92         (mosaicx (ceiling (sqrt nn)))
  93         (mosaic (make-array (list (* n mosaicx) (* n mosaicx))
  94                             :element-type '(unsigned-byte 8))))
  95    (with-open-file (*standard-output* "/dev/shm/a"
  96                                       :direction :output
  97                                       :if-exists :supersede)
  98      (dotimes (*nucleus-index* nn)
  99        (dotimes (i 10)
 100          (tagbody again
 101             (multiple-value-bind (min point)
 102                 (simplex-anneal:anneal (simplex-anneal:make-simplex
 103                                         (make-vec2 :x -1d0 :y -1d0) 1d0)
 104                                        #'merit-function
 105                                        ;; set temperature bigger than the
 106                                        ;; maxima in the bfp but smaller
 107                                        ;; than border-value
 108                                        :start-temperature 2.4d0
 109                                        :eps/m .02d0
 110                                        :itmax 1000
 111                                        :ftol 1d-3)
 112               (unless (<= min 100d0)
 113                 (go again))
 114               (let* ((x (aref point 0))
 115                      (y (aref point 1))
 116                      (ix (floor (* n (+ x 1)) 2))
 117                      (iy (floor (* n (+ y 1)) 2))
 118                      (mx (mod *nucleus-index* mosaicx))
 119                      (my (floor *nucleus-index* mosaicx)))
 120                 (incf (aref mosaic (+ (* n my) iy) (+ (* n mx) ix)))
 121                 (format t "min ~a~%" (list min ix iy))))))))
 122    (write-pgm "/home/martin/tmp/scan-mosaic-max.pgm" mosaic)))
 123
 124
 125 #+nil
 126 (illuminate-ray (lens:make-objective) *model* 0 :left .4d0 0d0 .2d0 :right)
 127
 128 #+nil
 129 (defmacro defstuff ()
 130   `(progn
 131      ,@(loop for i in '(*centers* *dims* *spheres* *psf* *conv-l* *conv-l-s*
 132                         *conv-plane* *conv-plane-s*
 133                         Lf L-plane*f)
 134           collect
 135             `(defparameter ,i nil))))
 136 #+nil
 137 (defstuff)
 138
 139 #+nil
 140 (write-pgm "/home/martin/tmp/fftw.pgm"
 141  (normalize-2-csf/ub8-abs
 142             (cross-section-xz
 143              (let ((a (ft (draw-sphere-csf 12.0 34 206 296)))
 144                    #+nil (b (ft (draw-sphere-csf 5.0 34 206 296))))
 145                a))))
 146 #+bla
 147 (defun init-model ()
 148   ;; find the centers of the nuclei and store into *centers*
 149   (multiple-value-bind (c ch dims)
 150       (find-centers)
 151     (declare (ignore ch))
 152     (defparameter *centers* c)
 153     (defparameter *dims* dims)
 154     (sb-ext:gc :full t))
 155
 156   ;; as a model of fluorophore concentration draw ovals around the
 157   ;; positions in *centers* and store into *spheres*
 158   (let ((spheres
 159          (destructuring-bind (z y x)
 160              *dims*
 161            (draw-ovals 12.0 *centers* z y x))))
 162     (defparameter *spheres* spheres)
 163     (write-pgm "/home/martin/tmp/comp0-spheres-cut.pgm"
 164                (normalize-2-csf/ub8-realpart
 165                (cross-section-xz *spheres*
 166                                  (vec-i-y (elt *centers* 21)))))
 167    (sb-ext:gc :full t))
 168
 169   ;; store the fluorophore concentration
 170   (save-stack-ub8 "/home/martin/tmp/spheres"
 171                   (normalize-3-csf/ub8-realpart *spheres*)))
 172
 173 #+nil
 174 (time (init-model))
 175
 176
 177 #+bla
 178 (defun init-psf ()
 179   ;; calculate intensity psf, make extend in z big enough to span the
 180   ;; full fluorophore concentration even when looking at the bottom
 181   ;; plane of it
 182   (let* ((dx .2)
 183          (dz 1.0)
 184          (psf (destructuring-bind (z y x)
 185                   *dims*
 186                 (declare (ignore y x))
 187                 (let ((r 100))
 188                  (psf:intensity-psf (* 2 z) r r (* z dz) (* r dx)
 189                                     :integrand-evaluations 40)))))
 190     (defparameter *psf* psf)
 191     (write-pgm "/home/martin/tmp/comp1-psf.pgm"
 192                (normalize-2-csf/ub8-realpart (cross-section-xz psf)))
 193     (sb-ext:gc :full t)))
 194 #+nil
 195 (time (init-psf))
 196
 197 #+bla
 198 (defun clem ()
 199   ;; Extract one specific plane from the fluorophore concentration
 200   ;; model and convolve with intensity psf. The result is the light
 201   ;; distribution in the sample.
 202   (destructuring-bind (z y x)
 203       (array-dimensions *spheres*)
 204     (declare (ignore z))
 205     (let* ((zz (vec-i-z (elt *centers* 31)))
 206            (current-slice-bbox (make-bbox :start (v 0d0 0d0 (* 1d0 zz))
 207                                           :end (v (* 1d0 (1- x))
 208                                                   (* 1d0 (1- y))
 209                                                   (* 1d0 zz))))
 210            (current-slice (extract-bbox3-cdf *spheres* current-slice-bbox)))
 211       (multiple-value-bind (conv conv-start)
 212           (convolve3-nocrop current-slice *psf*)
 213         (defparameter *conv-l* conv)
 214         (defparameter *conv-l-s* (v--i (make-vec-i :z zz)
 215                                        conv-start))
 216         (write-pgm "/home/martin/tmp/comp2-conv.pgm"
 217                    (normalize2-cdf/ub8-realpart
 218                     (cross-section-xz
 219                      conv
 220                      (vec-i-y (v+-i conv-start (elt *centers* 31)))))))
 221       (sb-ext:gc :full t)))
 222
 223   ;; store light distribution
 224   (save-stack-ub8 "/home/martin/tmp/conv-l"
 225                  (normalize3-cdf/ub8-realpart *conv-l*))
 226
 227   ;; multiply fluorophore concentration with light distribution, this
 228   ;; gives the excitation pattern in the sample
 229   (progn
 230     (defparameter Lf (.* *conv-l* *spheres* *conv-l-s*))
 231     (write-pgm "/home/martin/tmp/comp3-conv-lf.pgm"
 232               (normalize2-cdf/ub8-realpart
 233                (cross-section-xz Lf (vec-i-y (elt *centers* 31))))))
 234
 235   ;; estimate in-focus and out-of-focus excitation for this particular
 236   ;; excitation regime
 237   (destructuring-bind (z y x)
 238       (array-dimensions Lf)
 239     (let* ((zz (1- (floor z 2)))
 240            (in-focus (extract-bbox3-cdf Lf (make-bbox :start (v 0d0 0d0 (* 1d0 zz))
 241                                                       :end (v (* 1d0 (1- x))
 242                                                               (* 1d0 (1- y))
 243                                                               (* 1d0 zz))))))
 244       (save-stack-ub8 "/home/martin/tmp/Lf" (normalize3-cdf/ub8-realpart in-focus))
 245       (/ (mean-realpart in-focus)
 246          (mean-realpart Lf)))))
 247
 248 #+bla
 249 (defun widefield ()
 250   ;; convolve a plane with the psf
 251   (destructuring-bind (z y x)
 252       (array-dimensions *spheres*)
 253     (declare (ignore z))
 254     (let* ((zz (vec-i-z (elt *centers* 31)))
 255            (current-slice (make-array (list 1 y x)
 256                                       :element-type '(complex my-float)
 257                                       :initial-element (complex 1d0))))
 258       (multiple-value-bind (conv conv-start)
 259           (convolve3-nocrop current-slice *psf*)
 260         (defparameter *conv-plane* conv)
 261         (defparameter *conv-plane-s* (v--i (make-vec-i :z zz)
 262                                            conv-start))
 263         (write-pgm "/home/martin/tmp/comp4-conv-plane.pgm"
 264                    (normalize2-cdf/ub8-realpart
 265                     (cross-section-xz
 266                      conv
 267                      (vec-i-y (v+-i conv-start (elt *centers* 31)))))))
 268       (sb-ext:gc :full t)))
 269
 270   ;; multiply fluorophore concentration with light distribution, this
 271   ;; gives the excitation pattern in the sample
 272   (progn
 273     (defparameter L-plane*f (.* *conv-plane* *spheres* *conv-plane-s*))
 274     (write-pgm "/home/martin/tmp/comp5-conv-plane-lf.pgm"
 275                (normalize2-cdf/ub8-realpart
 276                 (cross-section-xz L-plane*f (vec-i-y (elt *centers* 31))))))
 277
 278   ;; estimate in-focus and out-of-focus excitation for this particular
 279   ;; excitation regime
 280   (destructuring-bind (z y x)
 281       (array-dimensions L-plane*f)
 282     (let* ((zz (1- (floor z 2)))
 283            (in-focus (extract-bbox3-cdf L-plane*f (make-bbox :start (v 0d0 0d0 (* 1d0 zz))
 284                                                              :end (v (* 1d0 (1- x))
 285                                                                      (* 1d0 (1- y))
 286                                                                      (* 1d0 zz))))))
 287       #+nil  (save-stack-ub8 "/home/martin/tmp/Lf" (normalize3-cdf/ub8-realpart in-focus))
 288       (/ (mean-realpart in-focus)
 289          (mean-realpart L-plane*f)))))
 290
 291 #+nil
 292 (time
 293  (progn
 294    (init-model)
 295    (init-psf))) ;; 7.5s
 296
 297 #+nil
 298 (time
 299  (clem)) ;; 8s, result: 6.5
 300
 301 #+nil
 302 (time
 303  (widefield)) ;; 5.7s result: 1.8
 304
 305 #||
 306 mkdir ~/tmp
 307 cp /home/martin/0519/MedianofCelegans-10-02-09-LSM700-t58.tif ~/tmp/med.tif
 308 cd ~/tmp/
 309 tiffsplit med.tif
 310 for i in *.tif ; do tifftopnm $i > `basename $i .tif`.pgm;done
 311 ||#
 312
 313
 314
 315
 316 #+bla
 317 (defun ensure-even (x)
 318   (declare (fixnum x)
 319            (values fixnum &optional))
 320   (if (eq 1 (mod x 2))
 321       (1+ x)
 322       x))
 323
 324
 325 #+nil
 326 (progn
 327  (defparameter *merge*
 328    (let ((a (make-array (array-dimensions *stack*)
 329                         :element-type '(unsigned-byte 8))))
 330      (destructuring-bind (z y x)
 331          (array-dimensions *stack*)
 332        (do-box (k j i 0 z 0 y 0 x)
 333          (setf (aref a k j i)
 334                (clamp (if (eq 0 (aref *blobs*
 335                                       k j i))
 336                           (* (aref *stack* k j i) 2)
 337                           0))))
 338        a)))
 339  (save-stack-ub "/home/martin/tmp/merge" *merge*))
 340
 341 ;; (- z k 1) (- y j 1) (- x i 1)
 342
 343
 344
 345 #+nil ;; model with isotropic pixels
 346 (save-stack-ub "/home/martin/tmp/iso/iso"
 347                (convert-vol (draw-scaled-spheres 2d0 *centers*)))
 348
 349 #+nil
 350 (save-stack-ub "/home/martin/tmp/blobs" *blobs*)
 351
 352 #+nil
 353 (save-scaled-stack 100 100 "/home/martin/tmp/sca" *blobs*)
 354
 355 #+nil ;; find maximum
 356 (reduce #'max (map 'vector #'abs (sb-ext:array-storage-vector *stack*)))
 357
 358 #+nil ;; scale to 1
 359 (destructuring-bind (z y x)
 360     (array-dimensions *stack*)
 361     (do-box (k j i 0 z 0 y 0 x)
 362       (setf (aref *stack* k j i)
 363             (/ (aref *stack* k j i) 257d0))))
 364
 365 #+nil
 366 (save-stack "/home/martin/tmp/stack" *stack*)
 367
 368
 369 #+nil ;; make a text image of the psf
 370 (let ((numerical-aperture 1.38d0))
 371   (psf:init :numerical-aperture numerical-aperture)
 372   (multiple-value-bind (u v)
 373       (psf:get-uv 0 1.5 3 :numerical-aperture numerical-aperture)
 374    (let* ((nu 31)
 375           (nv 61)
 376           (a (psf:integ-all nu nv  (/ u nu) (/ v nv)))
 377           (max (aref a 0 0))
 378           (scale (/ 99 max)))
 379      (destructuring-bind (uu vv)
 380          (array-dimensions a)
 381        (format t "~%")
 382        (dotimes (u uu)
 383          (dotimes (v vv)
 384            (format t "~2d" (truncate (* scale (aref a u v)))))
 385          (format t "~%"))))))
 386
 387
 388 #+nil
 389 (time
 390  (let ((z 18d0))
 391   (save-stack "/home/martin/tmp/psf"
 392               (fftshift3 (ft3 (sim-psf 128 128 128 z (* .5d0 z))))
 393               :function #'(lambda (x) (* .01d0 (abs x))))))
 394
 395 ;; worm stack 0.198 um in X and Y and 1 um in Z
 396
 397 #+nil
 398 (array-dimensions *blobs*)
 399
 400 #+nil ;; write ft of object
 401 (time
 402  (save-stack "/home/martin/tmp/kblobs" (fftshift3 (ft3 *blobs*))
 403              :function #'(lambda (x) (* 1d-4 (abs x)))))
 404
 405 #+nil ;; write ft of psf
 406 (time
 407  (destructuring-bind (z y x)
 408      (array-dimensions *blobs*)
 409    (save-stack "/home/martin/tmp/otf"
 410                (fftshift3 (ft3 (sim-psf z y x
 411                                         (* .198d0 z)
 412                                         (* .198d0
 413                                            (ceiling (* (sqrt 2d0) (max y x)))))))
 414               :function #'(lambda (x) (* 1d-1 (abs x))))))
 415
 416 #+nil
 417 (time
 418  (destructuring-bind (z y x)
 419      (array-dimensions *blobs*)
 420   (let ((psf (sim-psf z y x
 421                       (* .198d0 z) (* .198d0 (ceiling (* (sqrt 2d0) (max y x)))))))
 422     (save-stack "/home/martin/tmp/impsf"
 423                 (convolve3 *blobs* psf)
 424                 :function #'(lambda (x) (* 1d-8 (realpart x)))))))
 425
 426
 427
 428 #+nil ;; compare intensity and |e-field|^2
 429 (time
 430  (let ((z 128)
 431        (y 128)
 432        (x 128))
 433   (multiple-value-bind (e0 e1 e2)
 434       (psf:electric-field-psf z x y 10d0 5d0)
 435     (let ((intens (make-array (array-dimensions e0)
 436                               :element-type '(complex my-float))))
 437       (do-box (k j i 0 z 0 y 0 x)
 438         (setf (aref intens k j i) (complex (+ (psf::abs2 (aref e0 k j i))
 439                                               (psf::abs2 (aref e1 k j i))
 440                                               (psf::abs2 (aref e2 k j i))))))
 441       (let* ((k0 (fftshift3 (ft3 intens)))
 442              (k1 (fftshift3 (ft3 (psf:intensity-psf z y x 10d0 5d0))))
 443              (k- (make-array (array-dimensions k0)
 444                              :element-type '(complex my-float))))
 445         (do-box (k j i 0 z 0 y 0 x)
 446           (setf (aref k- k j i) (- (aref k0 k j i)
 447                                    (aref k1 k j i))))
 448         (save-stack "/home/martin/tmp/intens0"
 449                     k-
 450                     :function #'(lambda (x) (* 1d-1 (abs x))))
 451         (write-pgm (convert-img (cross-section-xz k-)
 452                                 #'(lambda (z) (* 1e-1 (abs z))))
 453                    "/home/martin/tmp/intens0xz.pgm"))))))
 454
 455
 456
 457 #+nil ;; find centers of nuclei 12.5s 3.1s
 458 (time
 459  (multiple-value-bind (c ch dims)
 460        (find-centers)
 461    (defparameter *centers* c)
 462       (defparameter *center-heights* ch)
 463    (defparameter *dims* dims)
 464    (sb-ext:gc :full t)))
 465
 466 #+nil ;; draw the spheres (squeezed in z) 9.7s 1.8s
 467 (time
 468  (let ((spheres
 469         (destructuring-bind (z y x)
 470             *dims*
 471           (draw-ovals 7d0 *centers* z y x))))
 472    (setf *spheres* spheres)
 473    #+nil (save-stack-ub8 "/home/martin/tmp/spheres" (normalize-vol *spheres*))
 474    (write-pgm (normalize-img (cross-section-xz *spheres*
 475                                                (vec-i-y (elt *centers* 31))))
 476               "/home/martin/tmp/spheres-cut.pgm")
 477    (sb-ext:gc :full t)))
 478
 479 #+nil ;; draw the spheres (squeezed in z) and emphasize one of them
 480 (time
 481  (let ((spheres
 482         (destructuring-bind (z y x)
 483             *dims*
 484           (let* ((dims (list z y x))
 485                  (points (make-array dims
 486                              :element-type '(complex my-float)))
 487                  (centers *centers*)
 488                  (radius 7d0)
 489                  (n (length centers)))
 490             (dotimes (i n)
 491               (let ((c (aref centers i)))
 492                 (setf (aref points
 493                             (vec-i-z c)
 494                             (vec-i-y c)
 495                             (vec-i-x c))
 496                       (complex (if (eq i 31)
 497                                    2d0
 498                                    1d0) 0d0))))
 499             (convolve3-circ points (fftshift3 (draw-oval radius z y x)))))))
 500    (save-stack-ub8 "/home/martin/tmp/spheres" (normalize-ub8 spheres))
 501    (sb-ext:gc :full t)))
 502
 503 #+nil ;; construct LCOS image
 504 (let ((coord (aref *centers* 31))
 505       (radius 7d0)
 506       (slice (make-array (array-dimensions *spheres*)
 507                          :element-type '(complex my-float))))
 508   (destructuring-bind (z y x)
 509       (array-dimensions *spheres*)
 510     ;; draw only the center sphere
 511     (let* ((xc (vec-i-x coord))
 512            (yc (vec-i-y coord))
 513            (zc (vec-i-z coord))
 514            (k  zc))
 515       (do-rectangle (j i 0 y 0 x)
 516        (let ((r (sqrt (+ (square (* 1d0 (- i xc)))
 517                          (square (* 1d0 (- j yc)))
 518                          (square (* 1d0 (- k zc)))))))
 519          (setf (aref slice k j i)
 520                (if (< r radius)
 521                    (complex 255d0)
 522                    (complex 0d0)))))))
 523   (defparameter *slice* slice)
 524   #+nil (save-stack-ub8 "/home/martin/tmp/slice" (convert-vol slice))
 525   (write-pgm (normalize-img (cross-section-xz slice (vec-i-y coord)))
 526              "/home/martin/tmp/slice-cut.pgm")
 527   (sb-ext:gc :full t))
 528
 529 #+bla (defparameter *bfp-circ-radius* .3d0)
 530 #+bla (defparameter *bfp-circ-center-x* .4d0 #+nil (- .999d0 *bfp-circ-radius*))
 531
 532 #+nil ;; 11.3s 2.6s
 533 (time
 534  (progn
 535   (angular-psf :x 80 :z 90
 536                :window-x *bfp-circ-center-x*
 537                :window-y 0d0 :window-radius *bfp-circ-radius*
 538                :numerical-aperture 1.38d0
 539                :immersion-index 1.515d0
 540                :pixel-size-x .1d0 :pixel-size-z .5d0
 541                :integrand-evaluations 160
 542                :debug t)
 543   nil))
 544
 545 #+nil ;; light distribution in the specimen
 546 ;; default resolution is isotropic 12 um /64 = 187.5 nm/pixel
 547 (time ;; 32.5s 5.4s
 548  (let* ((radius .2d0)
 549         (x .3d0)
 550         (xx 120)
 551         (zz 120)
 552         (dx .1d0)
 553         (dz .5d0)
 554         (psf (resample-half
 555               (angular-psf :window-x *bfp-circ-center-x*
 556                            :window-y 0d0
 557                            :window-radius *bfp-circ-radius*
 558                            :x (* 2 xx) :z (* 2 zz)
 559                            :pixel-size-x dx :pixel-size-z dz
 560                            :integrand-evaluations 200)))
 561         (dims (destructuring-bind (z y x)
 562                   *dims*
 563                 (list z y x))))
 564    (write-pgm (normalize-img (cross-section-xz psf))
 565               "/home/martin/tmp/small-psf-cut.pgm")
 566    (sb-ext:gc :full t)
 567    (defparameter *slice-x-psf* (convolve3 *slice* psf))
 568    (sb-ext:gc :full t)))
 569
 570 #+nil
 571 (defparameter *slice-x-psf* nil)
 572 #+nil
 573 (sb-ext:gc :full t)
 574 #+nil
 575 (write-pgm (normalize-img
 576             (cross-section-xz *slice-x-psf*
 577                               (vec-i-y (elt *centers* 31))))
 578            "/home/martin/tmp/slice-x-psf-cut.pgm")
 579
 580 #+nil
 581 (save-stack-ub8 "/home/martin/tmp/psf" (normalize-ub8 *psf*))
 582
 583
 584 #+nil ;; draw lines into the light distribution in the specimen
 585 (destructuring-bind (z y x)
 586     (array-dimensions *slice-x-psf*)
 587   (let ((coord (elt *centers* 31))
 588         (vol (normalize-ub8 *slice-x-psf*))
 589         (dx 2.d-4)
 590         (dz 1d-3))
 591     #+nil(draw-ray-into-vol (* dx (- (floor x 2) (vec-i-x coord)))
 592                        (* dx (- (floor y 2) (vec-i-y coord)))
 593                        -.6d0 0d0 vol)
 594     (loop for pos in (list (list (* (- (floor x 2) (- (vec-i-x coord) 7)) dx)
 595                                  (* (- (floor y 2) (vec-i-y coord)) dx))
 596                            (list (* (- (floor x 2) (+ (vec-i-x coord) 7)) dx)
 597                                  (* (- (floor y 2) (vec-i-y coord)) dx))) do
 598          (loop for angle in (list ;;-.010d0
 599                                   ;;-.6d0
 600                                   ;;(- (- *bfp-circ-center-x* *bfp-circ-radius*))
 601                                   ;;-.8d0
 602                                   ;;-.99d0
 603                              (- *bfp-circ-center-x*)
 604                                   ;;(- (+ *bfp-circ-center-x* *bfp-circ-radius*))
 605                                   ) do
 606               (draw-ray-into-vol (first pos) (second pos)
 607                                  angle 0d0
 608                                  vol
 609                                  :shift-z (- (vec-i-z coord)
 610                                              (floor z 2)))))
 611     nil
 612     #+nil (write-pgm (normalize-img
 613                 (cross-section-xz vol
 614                                   (vec-i-y (elt *centers* 31))))
 615                "/home/martin/tmp/slice-x-psf-lines-cut.pgm")
 616     (save-stack-ub8 "/home/martin/tmp/slice-x-psf" vol)))
 617
 618
 619 #+nil ;; excited fluorophores
 620 (progn
 621   (setf *slice-x-psf-times-spheres* (.* *spheres* *slice-x-psf*))
 622   (save-stack-ub8 "/home/martin/tmp/slice-x-psf-times-spheres"
 623                  (normalize-ub8 *slice-x-psf-times-spheres*)))
 624
 625
 626 #+nil ;; blur with detection psf
 627 (time
 628  (let* ((radius .5d0)
 629         (x (- 1d0 radius))
 630         (xx 80)
 631         (yy xx)
 632         (zz 128)
 633         (dx .2d0)
 634         (psf (psf:intensity-psf zz yy xx (* zz dx) (* xx dx)
 635                                 :integrand-evaluations 100))
 636         (dims (destructuring-bind (z y x)
 637                   *dims*
 638                 (list (* z 5) y x)))
 639         (psf-big (make-array dims
 640                              :element-type '(complex my-float))))
 641    (setf *psf-big* psf-big)
 642    (destructuring-bind (z y x)
 643        dims
 644      (let ((ox (- (floor x 2) (floor xx 2)))
 645            (oy (- (floor y 2) (floor yy 2)))
 646            (oz (- (floor z 2) (floor zz 2))))
 647        (do-box (k j i 0 zz 0 yy 0 xx)
 648          (setf (aref psf-big (+ oz k) (+ oy j) (+ ox i))
 649                (aref psf k j i)))))
 650    (save-stack-ub8 "/home/martin/tmp/psf-detect-big" (normalize-ub8 psf-big))
 651    (sb-ext:gc :full t)
 652    (defparameter *camera-volume* (convolve3-circ *slice-x-psf-times-spheres*
 653                                                  (fftshift3 psf-big)))
 654    (save-stack-ub8 "/home/martin/tmp/camera-volume"
 655                    (normalize-ub8 *camera-volume*))
 656    (sb-ext:gc :full t)))
 657
 658
 659
 660 #+nil ;; check convolution
 661 (time
 662  (let ((a (make-array (list 64 64 64)
 663                       :element-type '(complex my-float)))
 664        (b (psf:intensity-psf 64 64 64 20d0 20d0) #+nil (make-array (list 64 64 64)
 665                       :element-type '(complex my-float))))
 666    (setf (aref a 12 12 12) (complex 255d0))
 667 #+nil   (setf (aref b 0 0 0) (complex 255d0))
 668    (save-stack-ub8 "/home/martin/tmp/conv-test" (normalize-ub8 (convolve3-circ a (fftshift3 b))))))
 669
 670
 671 #+nil ;; output the xz cross section centered on a sphere in the middle
 672 (let ((coord (aref *centers* 30)))
 673   (write-pgm (normalize-img (cross-section-xz *spheres* (* 5 (vec-i-z coord))))
 674              "/home/martin/tmp/cut-spheres.pgm")
 675   (write-pgm (normalize-img (cross-section-xz *psf-big*))
 676              "/home/martin/tmp/cut-psf-big.pgm")
 677   (write-pgm (normalize-img (cross-section-xz *slice-x-psf* (* 5 (vec-i-z coord))))
 678              "/home/martin/tmp/cut-slice-x-psf.pgm")
 679   (write-pgm (normalize-img (cross-section-xz (.* *spheres* *slice-x-psf*) (* 5 (vec-i-z coord))))
 680              "/home/martin/tmp/cut-exfluo.pgm"))
 681
 682 #+nil
 683 (save-stack-ub8 "/home/martin/tmp/spheres" (normalize-ub8 *spheres*))
 684
 685 #+nil
 686 (let ((sli (make-array (array-dimensions *spheres*)
 687                        :element-type '(complex my-float))))
 688   (destructuring-bind (z y x)
 689       (array-dimensions *spheres*)
 690     (do-box (k j i 0 z 0 y 0 x)
 691       (setf (aref sli k j i)
 692             (aref *spheres* k j i)))
 693     (let ((k (* 5 (vec-i-z (aref *centers* 30)))))
 694      (do-rectangle (j i 0 y 0 x)
 695        (setf (aref sli k j i)
 696              (* 10 (aref sli k j i)))))
 697     (defparameter *sli* sli))
 698   (save-stack-ub8 "/home/martin/tmp/spheres" (normalize-ub8 *sli*)))
 699
 700
 701
 702 #+nil
 703 (let ((a (sb-ext:array-storage-vector *slice*)))
 704   (reduce #'max (map 'vector #'abs a)))
 705
 706 #+nil
 707 (sb-ext:gc :full t)
 708
 709 #+nil ;; model with unscaled spheres
 710 (defparameter *blobs*
 711   (destructuring-bind (z y x)
 712       *dims*
 713     (draw-spheres 7d0 *centers* (* 5 z) y x)))
 714
 715
 716 #+nil ;; print ft of angular psf
 717 (time (let* ((radius .5d0)
 718         (x (- 1d0 radius))
 719         (psf (angular-psf x 0d0 radius)))
 720    (write-pgm (normalize-ub8 (cross-section-xz (fftshift3 (ft3 psf))))
 721               "/home/martin/tmp/cut-intens.pgm")))
 722 #+nil
 723 (let* ((intens (psf:intensity-psf 64 64 64 10d0 5d0))
 724        (k0 intens #+nil(fftshift3 (ft3 intens))))
 725   (save-stack "/home/martin/tmp/intens1" k0
 726               :function #'(lambda (x) (* 1d-5 (abs x))))
 727   (write-pgm (convert-img
 728               (cross-section-xz k0)
 729               #'(lambda (z) (* 1d-4 (abs z))))
 730              "/home/martin/tmp/intens1xz.pgm"))
 731
 732 #+nil
 733 (time
 734  (destructuring-bind (z y x)
 735      (array-dimensions *blobs*)
 736    (save-stack "/home/martin/tmp/blobs"
 737                *blobs*
 738                )))
 739
 740
 741 #+nil ;; clean up the garbage
 742 (sb-ext:gc :full t)
 743
 744
 745 #+nil
 746 (sb-vm:memory-usage :print-spaces t :count-spaces t)
 747 #+nil
 748 (sb-vm:memory-usage)
 749
 750 #+nil
 751 (sb-vm:instance-usage :dynamic)
 752
 753 #+nil
 754 (sb-vm:list-allocated-objects :dynamic)
 755
 756 #+nil
 757 (sb-vm:print-allocated-objects :dynamic)
 758
 759 #+nil
 760 (format t "~a~%" (sb-vm::type-breakdown :dynamic))
 761
 762
 763 ;;
 764
 765
 766
 767