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