run.lisp

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