run.lisp

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