run.lisp

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