run.lisp

   1 #+nil
   2 (progn
   3    (require :vol)
   4    (require :raytrace)
   5    (require :bresenham)
   6    (require :psf))
   7 (defpackage :run
   8   (:use :cl :vector :vol :raytrace :bresenham))
   9 (in-package :run)
  10
  11 (deftype my-float-helper ()
  12   `single-float)
  13
  14 (deftype my-float (&optional (low '*) (high '*))
  15   `(single-float ,(if (eq low '*)
  16                       '*
  17                       (coerce low 'my-float-helper))
  18                  ,(if (eq high '*)
  19                       '*
  20                       (coerce high 'my-float-helper))))
  21
  22 (defconstant +one+ #.(coerce 1 'my-float))
  23
  24 #+nil
  25 (init-ft)
  26
  27 (defmacro defstuff ()
  28   `(progn
  29      ,@(loop for i in '(*centers* *dims* *spheres* *psf* *conv-l* *conv-l-s*
  30                         *conv-plane* *conv-plane-s*
  31                         Lf L-plane*f)
  32           collect
  33             `(defparameter ,i nil))))
  34
  35 (defstuff)
  36
  37
  38 (defun draw-spheres (radius centers z y x)
  39   "put points into the centers of nuclei and convolve a sphere around each"
  40   (declare (single-float radius)
  41            ((simple-array vec-i 1) centers)
  42            (fixnum z y x)
  43            (values (simple-array (complex my-float) 3) &optional))
  44   (let* ((dims (list z y x))
  45          (points (make-array dims
  46                              :element-type '(complex my-float)))
  47          (n (length centers)))
  48     (dotimes (i n)
  49       (let ((c (aref centers i)))
  50         (setf (aref points
  51                     (* 5 (vec-i-z c))
  52                     (vec-i-y c)
  53                     (vec-i-x c))
  54               (complex +one+))))
  55     (convolve-circ points
  56                    (fftshift (draw-sphere-csf radius z y x)))))
  57
  58 (defun draw-ovals (radius centers z y x)
  59   (declare (single-float radius)
  60            ((simple-array vec-i 1) centers)
  61            (fixnum z y x)
  62            (values (simple-array (complex my-float) 3) &optional))
  63   (let* ((dims (list z y x))
  64          (points (make-array dims
  65                              :element-type '(complex my-float)))
  66          (n (length centers)))
  67     (dotimes (i n)
  68       (let ((c (aref centers i)))
  69         (setf (aref points
  70                     (vec-i-z c)
  71                     (vec-i-y c)
  72                     (vec-i-x c))
  73               (complex +one+))))
  74     (convolve-circ points (fftshift (draw-oval-csf radius z y x)))))
  75
  76 (defun draw-indexed-ovals (radius centers z y x)
  77   "The first oval contains the value 1 the second 2, ..."
  78   (declare (single-float radius)
  79            ((simple-array vec-i 1) centers)
  80            (fixnum z y x)
  81            (values (simple-array (complex my-float) 3) &optional))
  82   (let* ((dims (list z y x))
  83          (points (make-array dims
  84                              :element-type '(complex my-float)))
  85          (n (length centers)))
  86     (dotimes (i n)
  87       (let ((c (aref centers i)))
  88         (setf (aref points
  89                     (vec-i-z c)
  90                     (vec-i-y c)
  91                     (vec-i-x c))
  92               (complex (+ +one+ i)))))
  93     (convolve-circ points
  94                    (fftshift (draw-oval-csf radius z y x)))))
  95
  96
  97 ;; to find centers of cells do a convolution with a sphere
  98 (defun find-centers ()
  99   (declare (values (simple-array vec-i 1)
 100                    (simple-array my-float 1)
 101                    cons &optional))
 102   (let* ((stack-byte (read-stack "/home/martin/tmp/xa*.pgm"))
 103          (dims (array-dimensions stack-byte))
 104          (stack (make-array dims :element-type '(complex my-float))))
 105     (destructuring-bind (z y x) dims
 106       (do-region ((k j i) (z y x))
 107         (setf (aref stack k j i) (complex (+ (* #.(coerce .43745 'my-float) k)
 108                                              (aref stack-byte k j i)))))
 109       ;; find centers of cells by convolving with sphere, actually an
 110       ;; oval because the z resolution is smaller than the transversal
 111       (let* ((conv (convolve-circ
 112                     stack
 113                     (fftshift
 114                      (#.(cond
 115                           ((subtypep 'my-float 'single-float) 'draw-oval-csf)
 116                           ((subtypep 'my-float 'double-float) 'draw-oval-cdf))
 117                         12.0 z y x))))
 118              (cv (convert conv 'sf 'realpart))
 119              (centers nil)
 120              (center-heights nil))
 121         (do-region ((k j i) ((- z 3) (- y 1) (- x 1)) (6 1 1))
 122           (macrolet ((c (a b c)
 123                        `(aref cv (+ k ,a) (+ j ,b) (+ i ,c))))
 124             (let ((v (c 0 0 0)))
 125               (when (and (< (c 0 0 -1) v)
 126                          (< (c 0 0 1) v)
 127                          (< (c 0 -1 0) v)
 128                          (< (c 0 1 0) v)
 129                          (< (c -1 0 0) v)
 130                          (< (c 1 0 0) v))
 131                 (push (make-vec-i :z k :y j :x i) centers)
 132                 (push v center-heights)))))
 133         (let ((c (make-array (length centers)
 134                              :element-type 'vec-i
 135                              :initial-contents centers))
 136               (ch (make-array (length center-heights)
 137                               :element-type 'my-float
 138                               :initial-contents center-heights)))
 139           (values c ch dims))))))
 140
 141 #+nil
 142 (sb-ext:gc :full t)
 143 #+nil
 144 (time
 145  (progn
 146    (format t "~a~%" (multiple-value-list (find-centers)))
 147    nil))
 148
 149 #+nil
 150 (write-pgm "/home/martin/tmp/fftw.pgm"
 151  (normalize-2-csf/ub8-abs
 152             (cross-section-xz
 153              (let ((a (ft (draw-sphere-csf 12.0 34 206 296)))
 154                    #+nil (b (ft (draw-sphere-csf 5.0 34 206 296))))
 155                a))))
 156
 157 (defun init-model ()
 158   ;; find the centers of the nuclei and store into *centers*
 159   (multiple-value-bind (c ch dims)
 160       (find-centers)
 161     (declare (ignore ch))
 162     (defparameter *centers* c)
 163     (defparameter *dims* dims)
 164     (sb-ext:gc :full t))
 165
 166   ;; as a model of fluorophore concentration draw ovals around the
 167   ;; positions in *centers* and store into *spheres*
 168   (let ((spheres
 169          (destructuring-bind (z y x)
 170              *dims*
 171            (draw-ovals 12.0 *centers* z y x))))
 172     (defparameter *spheres* spheres)
 173     (write-pgm "/home/martin/tmp/comp0-spheres-cut.pgm"
 174                (normalize-2-csf/ub8-realpart
 175                (cross-section-xz *spheres*
 176                                  (vec-i-y (elt *centers* 21)))))
 177    (sb-ext:gc :full t))
 178
 179   ;; store the fluorophore concentration
 180   (save-stack-ub8 "/home/martin/tmp/spheres"
 181                   (normalize-3-csf/ub8-realpart *spheres*)))
 182
 183 #+nil
 184 (time (init-model))
 185
 186
 187
 188 (defun init-psf ()
 189   ;; calculate intensity psf, make extend in z big enough to span the
 190   ;; full fluorophore concentration even when looking at the bottom
 191   ;; plane of it
 192   (let* ((dx .2)
 193          (dz 1.0)
 194          (psf (destructuring-bind (z y x)
 195                   *dims*
 196                 (declare (ignore y x))
 197                 (let ((r 100))
 198                  (psf:intensity-psf (* 2 z) r r (* z dz) (* r dx)
 199                                     :integrand-evaluations 40)))))
 200     (defparameter *psf* psf)
 201     (write-pgm "/home/martin/tmp/comp1-psf.pgm"
 202                (normalize-2-csf/ub8-realpart (cross-section-xz psf)))
 203     (sb-ext:gc :full t)))
 204 #+nil
 205 (time (init-psf))
 206
 207
 208 (defun clem ()
 209   ;; Extract one specific plane from the fluorophore concentration
 210   ;; model and convolve with intensity psf. The result is the light
 211   ;; distribution in the sample.
 212   (destructuring-bind (z y x)
 213       (array-dimensions *spheres*)
 214     (declare (ignore z))
 215     (let* ((zz (vec-i-z (elt *centers* 31)))
 216            (current-slice-bbox (make-bbox :start (v 0d0 0d0 (* 1d0 zz))
 217                                           :end (v (* 1d0 (1- x))
 218                                                   (* 1d0 (1- y))
 219                                                   (* 1d0 zz))))
 220            (current-slice (extract-bbox3-cdf *spheres* current-slice-bbox)))
 221       (multiple-value-bind (conv conv-start)
 222           (convolve3-nocrop current-slice *psf*)
 223         (defparameter *conv-l* conv)
 224         (defparameter *conv-l-s* (v--i (make-vec-i :z zz)
 225                                        conv-start))
 226         (write-pgm "/home/martin/tmp/comp2-conv.pgm"
 227                    (normalize2-cdf/ub8-realpart
 228                     (cross-section-xz
 229                      conv
 230                      (vec-i-y (v+-i conv-start (elt *centers* 31)))))))
 231       (sb-ext:gc :full t)))
 232
 233   ;; store light distribution
 234   (save-stack-ub8 "/home/martin/tmp/conv-l"
 235                  (normalize3-cdf/ub8-realpart *conv-l*))
 236
 237   ;; multiply fluorophore concentration with light distribution, this
 238   ;; gives the excitation pattern in the sample
 239   (progn
 240     (defparameter Lf (.* *conv-l* *spheres* *conv-l-s*))
 241     (write-pgm "/home/martin/tmp/comp3-conv-lf.pgm"
 242               (normalize2-cdf/ub8-realpart
 243                (cross-section-xz Lf (vec-i-y (elt *centers* 31))))))
 244
 245   ;; estimate in-focus and out-of-focus excitation for this particular
 246   ;; excitation regime
 247   (destructuring-bind (z y x)
 248       (array-dimensions Lf)
 249     (let* ((zz (1- (floor z 2)))
 250            (in-focus (extract-bbox3-cdf Lf (make-bbox :start (v 0d0 0d0 (* 1d0 zz))
 251                                                       :end (v (* 1d0 (1- x))
 252                                                               (* 1d0 (1- y))
 253                                                               (* 1d0 zz))))))
 254       (save-stack-ub8 "/home/martin/tmp/Lf" (normalize3-cdf/ub8-realpart in-focus))
 255       (/ (mean-realpart in-focus)
 256          (mean-realpart Lf)))))
 257
 258
 259 (defun widefield ()
 260   ;; convolve a plane with the psf
 261   (destructuring-bind (z y x)
 262       (array-dimensions *spheres*)
 263     (declare (ignore z))
 264     (let* ((zz (vec-i-z (elt *centers* 31)))
 265            (current-slice (make-array (list 1 y x)
 266                                       :element-type '(complex my-float)
 267                                       :initial-element (complex 1d0))))
 268       (multiple-value-bind (conv conv-start)
 269           (convolve3-nocrop current-slice *psf*)
 270         (defparameter *conv-plane* conv)
 271         (defparameter *conv-plane-s* (v--i (make-vec-i :z zz)
 272                                            conv-start))
 273         (write-pgm "/home/martin/tmp/comp4-conv-plane.pgm"
 274                    (normalize2-cdf/ub8-realpart
 275                     (cross-section-xz
 276                      conv
 277                      (vec-i-y (v+-i conv-start (elt *centers* 31)))))))
 278       (sb-ext:gc :full t)))
 279
 280   ;; multiply fluorophore concentration with light distribution, this
 281   ;; gives the excitation pattern in the sample
 282   (progn
 283     (defparameter L-plane*f (.* *conv-plane* *spheres* *conv-plane-s*))
 284     (write-pgm "/home/martin/tmp/comp5-conv-plane-lf.pgm"
 285                (normalize2-cdf/ub8-realpart
 286                 (cross-section-xz L-plane*f (vec-i-y (elt *centers* 31))))))
 287
 288   ;; estimate in-focus and out-of-focus excitation for this particular
 289   ;; excitation regime
 290   (destructuring-bind (z y x)
 291       (array-dimensions L-plane*f)
 292     (let* ((zz (1- (floor z 2)))
 293            (in-focus (extract-bbox3-cdf L-plane*f (make-bbox :start (v 0d0 0d0 (* 1d0 zz))
 294                                                              :end (v (* 1d0 (1- x))
 295                                                                      (* 1d0 (1- y))
 296                                                                      (* 1d0 zz))))))
 297       #+nil  (save-stack-ub8 "/home/martin/tmp/Lf" (normalize3-cdf/ub8-realpart in-focus))
 298       (/ (mean-realpart in-focus)
 299          (mean-realpart L-plane*f)))))
 300
 301 #+nil
 302 (time
 303  (progn
 304    (init-model)
 305    (init-psf))) ;; 7.5s
 306
 307 #+nil
 308 (time
 309  (clem)) ;; 8s, result: 6.5
 310
 311 #+nil
 312 (time
 313  (widefield)) ;; 5.7s result: 1.8
 314
 315 #||
 316 mkdir ~/tmp
 317 cp /home/martin/0519/MedianofCelegans-10-02-09-LSM700-t58.tif ~/tmp/med.tif
 318 cd ~/tmp/
 319 tiffsplit med.tif
 320 for i in *.tif ; do tifftopnm $i > `basename $i .tif`.pgm;done
 321 ||#
 322
 323
 324
 325
 326 (declaim (ftype (function (fixnum)
 327                           (values fixnum &optional))
 328                 ensure-even))
 329 (defun ensure-even (x)
 330   (if (eq 1 (mod x 2))
 331       (1+ x)
 332       x))
 333
 334 (declaim (ftype (function (my-float (simple-array vec-i 1)
 335                                         &key (:div-x my-float)
 336                                         (:div-y my-float)
 337                                         (:div-z my-float))
 338                           (values (simple-array (complex my-float) 3)
 339                                   &optional))
 340                 draw-scaled-spheres))
 341 ;; put points into the centers of nuclei and convolve a sphere around each
 342 (defun draw-scaled-spheres (radius centers &key (div-x 5d0)
 343                      (div-y 5d0)
 344                      (div-z 1d0))
 345   (let* ((max-x (floor (reduce #'max
 346                                (map '(simple-array fixnum 1) #'vec-i-x
 347                                     centers)) div-x))
 348          (max-y (floor (reduce #'max
 349                                (map '(simple-array fixnum 1) #'vec-i-y
 350                                     centers)) div-y))
 351          (max-z (floor (reduce #'max
 352                                (map '(simple-array fixnum 1) #'vec-i-z
 353                                     centers)) div-z))
 354          (cr (ceiling radius))
 355          (rh (floor radius 2))
 356          (x (ensure-even (+ max-x cr)))
 357          (y (ensure-even (+ max-y cr)))
 358          (z (ensure-even (+ max-z cr)))
 359          (dims (list z y x))
 360          (points (make-array dims
 361                              :element-type '(complex my-float))))
 362     (loop for c across centers do
 363          (setf (aref points
 364                      (- (floor (vec-i-z c) div-z) rh)
 365                      (- (floor (vec-i-y c) div-y) rh)
 366                      (- (floor (vec-i-x c) div-x) rh))
 367                (complex 1d0 0d0)))
 368     (convolve3-circ points
 369                     (convert3-ub8/cdf-complex (draw-sphere-ub8 radius z y x)))))
 370
 371
 372
 373
 374
 375 #+nil
 376 (progn
 377  (defparameter *merge*
 378    (let ((a (make-array (array-dimensions *stack*)
 379                         :element-type '(unsigned-byte 8))))
 380      (destructuring-bind (z y x)
 381          (array-dimensions *stack*)
 382        (do-box (k j i 0 z 0 y 0 x)
 383          (setf (aref a k j i)
 384                (clamp (if (eq 0 (aref *blobs*
 385                                       k j i))
 386                           (* (aref *stack* k j i) 2)
 387                           0))))
 388        a)))
 389  (save-stack-ub "/home/martin/tmp/merge" *merge*))
 390
 391 ;; (- z k 1) (- y j 1) (- x i 1)
 392
 393
 394
 395 #+nil ;; model with isotropic pixels
 396 (save-stack-ub "/home/martin/tmp/iso/iso"
 397                (convert-vol (draw-scaled-spheres 2d0 *centers*)))
 398
 399 #+nil
 400 (save-stack-ub "/home/martin/tmp/blobs" *blobs*)
 401
 402 #+nil
 403 (save-scaled-stack 100 100 "/home/martin/tmp/sca" *blobs*)
 404
 405 #+nil ;; find maximum
 406 (reduce #'max (map 'vector #'abs (sb-ext:array-storage-vector *stack*)))
 407
 408 #+nil ;; scale to 1
 409 (destructuring-bind (z y x)
 410     (array-dimensions *stack*)
 411     (do-box (k j i 0 z 0 y 0 x)
 412       (setf (aref *stack* k j i)
 413             (/ (aref *stack* k j i) 257d0))))
 414
 415 #+nil
 416 (save-stack "/home/martin/tmp/stack" *stack*)
 417
 418
 419 #+nil ;; make a text image of the psf
 420 (let ((numerical-aperture 1.38d0))
 421   (psf:init :numerical-aperture numerical-aperture)
 422   (multiple-value-bind (u v)
 423       (psf:get-uv 0 1.5 3 :numerical-aperture numerical-aperture)
 424    (let* ((nu 31)
 425           (nv 61)
 426           (a (psf:integ-all nu nv  (/ u nu) (/ v nv)))
 427           (max (aref a 0 0))
 428           (scale (/ 99 max)))
 429      (destructuring-bind (uu vv)
 430          (array-dimensions a)
 431        (format t "~%")
 432        (dotimes (u uu)
 433          (dotimes (v vv)
 434            (format t "~2d" (truncate (* scale (aref a u v)))))
 435          (format t "~%"))))))
 436
 437
 438 #+nil
 439 (time
 440  (let ((z 18d0))
 441   (save-stack "/home/martin/tmp/psf"
 442               (fftshift3 (ft3 (sim-psf 128 128 128 z (* .5d0 z))))
 443               :function #'(lambda (x) (* .01d0 (abs x))))))
 444
 445 ;; worm stack 0.198 um in X and Y and 1 um in Z
 446
 447 #+nil
 448 (array-dimensions *blobs*)
 449
 450 #+nil ;; write ft of object
 451 (time
 452  (save-stack "/home/martin/tmp/kblobs" (fftshift3 (ft3 *blobs*))
 453              :function #'(lambda (x) (* 1d-4 (abs x)))))
 454
 455 #+nil ;; write ft of psf
 456 (time
 457  (destructuring-bind (z y x)
 458      (array-dimensions *blobs*)
 459    (save-stack "/home/martin/tmp/otf"
 460                (fftshift3 (ft3 (sim-psf z y x
 461                                         (* .198d0 z)
 462                                         (* .198d0
 463                                            (ceiling (* (sqrt 2d0) (max y x)))))))
 464               :function #'(lambda (x) (* 1d-1 (abs x))))))
 465
 466 (declaim (ftype (function ((simple-array (complex my-float) (* * *)))
 467                           (values (simple-array (complex my-float) (* * *))
 468                                   &optional))
 469                 zshift3))
 470 (defun zshift3 (in)
 471   (let ((out (make-array (array-dimensions in)
 472                          :element-type '(complex my-float))))
 473    (destructuring-bind (w2 w1 w0)
 474        (array-dimensions in)
 475      (dotimes (k w2)
 476        (dotimes (j w1)
 477          (dotimes (i w0)
 478            (let* ((kk (if (> k (/ w2 2))
 479                           (+ w2 (/ w2 2) (- k))
 480                           (- (/ w2 2) k))))
 481              (setf (aref out k j i)
 482                    (aref in kk j i))
 483              nil)))))
 484    out))
 485
 486
 487
 488 #+nil
 489 (time
 490  (destructuring-bind (z y x)
 491      (array-dimensions *blobs*)
 492   (let ((psf (sim-psf z y x
 493                       (* .198d0 z) (* .198d0 (ceiling (* (sqrt 2d0) (max y x)))))))
 494     (save-stack "/home/martin/tmp/impsf"
 495                 (convolve3 *blobs* psf)
 496                 :function #'(lambda (x) (* 1d-8 (realpart x)))))))
 497
 498
 499
 500 #+nil ;; compare intensity and |e-field|^2
 501 (time
 502  (let ((z 128)
 503        (y 128)
 504        (x 128))
 505   (multiple-value-bind (e0 e1 e2)
 506       (psf:electric-field-psf z x y 10d0 5d0)
 507     (let ((intens (make-array (array-dimensions e0)
 508                               :element-type '(complex my-float))))
 509       (do-box (k j i 0 z 0 y 0 x)
 510         (setf (aref intens k j i) (complex (+ (psf::abs2 (aref e0 k j i))
 511                                               (psf::abs2 (aref e1 k j i))
 512                                               (psf::abs2 (aref e2 k j i))))))
 513       (let* ((k0 (fftshift3 (ft3 intens)))
 514              (k1 (fftshift3 (ft3 (psf:intensity-psf z y x 10d0 5d0))))
 515              (k- (make-array (array-dimensions k0)
 516                              :element-type '(complex my-float))))
 517         (do-box (k j i 0 z 0 y 0 x)
 518           (setf (aref k- k j i) (- (aref k0 k j i)
 519                                    (aref k1 k j i))))
 520         (save-stack "/home/martin/tmp/intens0"
 521                     k-
 522                     :function #'(lambda (x) (* 1d-1 (abs x))))
 523         (write-pgm (convert-img (cross-section-xz k-)
 524                                 #'(lambda (z) (* 1e-1 (abs z))))
 525                    "/home/martin/tmp/intens0xz.pgm"))))))
 526
 527
 528
 529 #+nil ;; find centers of nuclei 12.5s 3.1s
 530 (time
 531  (multiple-value-bind (c ch dims)
 532        (find-centers)
 533    (defparameter *centers* c)
 534       (defparameter *center-heights* ch)
 535    (defparameter *dims* dims)
 536    (sb-ext:gc :full t)))
 537
 538 #+nil ;; draw the spheres (squeezed in z) 9.7s 1.8s
 539 (time
 540  (let ((spheres
 541         (destructuring-bind (z y x)
 542             *dims*
 543           (draw-ovals 7d0 *centers* z y x))))
 544    (setf *spheres* spheres)
 545    #+nil (save-stack-ub8 "/home/martin/tmp/spheres" (normalize-vol *spheres*))
 546    (write-pgm (normalize-img (cross-section-xz *spheres*
 547                                                (vec-i-y (elt *centers* 31))))
 548               "/home/martin/tmp/spheres-cut.pgm")
 549    (sb-ext:gc :full t)))
 550
 551 #+nil ;; draw the spheres (squeezed in z) and emphasize one of them
 552 (time
 553  (let ((spheres
 554         (destructuring-bind (z y x)
 555             *dims*
 556           (let* ((dims (list z y x))
 557                  (points (make-array dims
 558                              :element-type '(complex my-float)))
 559                  (centers *centers*)
 560                  (radius 7d0)
 561                  (n (length centers)))
 562             (dotimes (i n)
 563               (let ((c (aref centers i)))
 564                 (setf (aref points
 565                             (vec-i-z c)
 566                             (vec-i-y c)
 567                             (vec-i-x c))
 568                       (complex (if (eq i 31)
 569                                    2d0
 570                                    1d0) 0d0))))
 571             (convolve3-circ points (fftshift3 (draw-oval radius z y x)))))))
 572    (save-stack-ub8 "/home/martin/tmp/spheres" (normalize-ub8 spheres))
 573    (sb-ext:gc :full t)))
 574
 575 #+nil ;; construct LCOS image
 576 (let ((coord (aref *centers* 31))
 577       (radius 7d0)
 578       (slice (make-array (array-dimensions *spheres*)
 579                          :element-type '(complex my-float))))
 580   (destructuring-bind (z y x)
 581       (array-dimensions *spheres*)
 582     ;; draw only the center sphere
 583     (let* ((xc (vec-i-x coord))
 584            (yc (vec-i-y coord))
 585            (zc (vec-i-z coord))
 586            (k  zc))
 587       (do-rectangle (j i 0 y 0 x)
 588        (let ((r (sqrt (+ (square (* 1d0 (- i xc)))
 589                          (square (* 1d0 (- j yc)))
 590                          (square (* 1d0 (- k zc)))))))
 591          (setf (aref slice k j i)
 592                (if (< r radius)
 593                    (complex 255d0)
 594                    (complex 0d0)))))))
 595   (defparameter *slice* slice)
 596   #+nil (save-stack-ub8 "/home/martin/tmp/slice" (convert-vol slice))
 597   (write-pgm (normalize-img (cross-section-xz slice (vec-i-y coord)))
 598              "/home/martin/tmp/slice-cut.pgm")
 599   (sb-ext:gc :full t))
 600
 601 (defparameter *bfp-circ-radius* .3d0)
 602 (defparameter *bfp-circ-center-x* .4d0 #+nil (- .999d0 *bfp-circ-radius*))
 603
 604 #+nil ;; 11.3s 2.6s
 605 (time
 606  (progn
 607   (angular-psf :x 80 :z 90
 608                :window-x *bfp-circ-center-x*
 609                :window-y 0d0 :window-radius *bfp-circ-radius*
 610                :numerical-aperture 1.38d0
 611                :immersion-index 1.515d0
 612                :pixel-size-x .1d0 :pixel-size-z .5d0
 613                :integrand-evaluations 160
 614                :debug t)
 615   nil))
 616
 617 #+nil ;; light distribution in the specimen
 618 ;; default resolution is isotropic 12 um /64 = 187.5 nm/pixel
 619 (time ;; 32.5s 5.4s
 620  (let* ((radius .2d0)
 621         (x .3d0)
 622         (xx 120)
 623         (zz 120)
 624         (dx .1d0)
 625         (dz .5d0)
 626         (psf (resample-half
 627               (angular-psf :window-x *bfp-circ-center-x*
 628                            :window-y 0d0
 629                            :window-radius *bfp-circ-radius*
 630                            :x (* 2 xx) :z (* 2 zz)
 631                            :pixel-size-x dx :pixel-size-z dz
 632                            :integrand-evaluations 200)))
 633         (dims (destructuring-bind (z y x)
 634                   *dims*
 635                 (list z y x))))
 636    (write-pgm (normalize-img (cross-section-xz psf))
 637               "/home/martin/tmp/small-psf-cut.pgm")
 638    (sb-ext:gc :full t)
 639    (defparameter *slice-x-psf* (convolve3 *slice* psf))
 640    (sb-ext:gc :full t)))
 641
 642 #+nil
 643 (defparameter *slice-x-psf* nil)
 644 #+nil
 645 (sb-ext:gc :full t)
 646 #+nil
 647 (write-pgm (normalize-img
 648             (cross-section-xz *slice-x-psf*
 649                               (vec-i-y (elt *centers* 31))))
 650            "/home/martin/tmp/slice-x-psf-cut.pgm")
 651
 652 #+nil
 653 (save-stack-ub8 "/home/martin/tmp/psf" (normalize-ub8 *psf*))
 654
 655
 656 #+nil ;; draw lines into the light distribution in the specimen
 657 (destructuring-bind (z y x)
 658     (array-dimensions *slice-x-psf*)
 659   (let ((coord (elt *centers* 31))
 660         (vol (normalize-ub8 *slice-x-psf*))
 661         (dx 2.d-4)
 662         (dz 1d-3))
 663     #+nil(draw-ray-into-vol (* dx (- (floor x 2) (vec-i-x coord)))
 664                        (* dx (- (floor y 2) (vec-i-y coord)))
 665                        -.6d0 0d0 vol)
 666     (loop for pos in (list (list (* (- (floor x 2) (- (vec-i-x coord) 7)) dx)
 667                                  (* (- (floor y 2) (vec-i-y coord)) dx))
 668                            (list (* (- (floor x 2) (+ (vec-i-x coord) 7)) dx)
 669                                  (* (- (floor y 2) (vec-i-y coord)) dx))) do
 670          (loop for angle in (list ;;-.010d0
 671                                   ;;-.6d0
 672                                   ;;(- (- *bfp-circ-center-x* *bfp-circ-radius*))
 673                                   ;;-.8d0
 674                                   ;;-.99d0
 675                              (- *bfp-circ-center-x*)
 676                                   ;;(- (+ *bfp-circ-center-x* *bfp-circ-radius*))
 677                                   ) do
 678               (draw-ray-into-vol (first pos) (second pos)
 679                                  angle 0d0
 680                                  vol
 681                                  :shift-z (- (vec-i-z coord)
 682                                              (floor z 2)))))
 683     nil
 684     #+nil (write-pgm (normalize-img
 685                 (cross-section-xz vol
 686                                   (vec-i-y (elt *centers* 31))))
 687                "/home/martin/tmp/slice-x-psf-lines-cut.pgm")
 688     (save-stack-ub8 "/home/martin/tmp/slice-x-psf" vol)))
 689
 690
 691 #+nil ;; excited fluorophores
 692 (progn
 693   (setf *slice-x-psf-times-spheres* (.* *spheres* *slice-x-psf*))
 694   (save-stack-ub8 "/home/martin/tmp/slice-x-psf-times-spheres"
 695                  (normalize-ub8 *slice-x-psf-times-spheres*)))
 696
 697
 698 #+nil ;; blur with detection psf
 699 (time
 700  (let* ((radius .5d0)
 701         (x (- 1d0 radius))
 702         (xx 80)
 703         (yy xx)
 704         (zz 128)
 705         (dx .2d0)
 706         (psf (psf:intensity-psf zz yy xx (* zz dx) (* xx dx)
 707                                 :integrand-evaluations 100))
 708         (dims (destructuring-bind (z y x)
 709                   *dims*
 710                 (list (* z 5) y x)))
 711         (psf-big (make-array dims
 712                              :element-type '(complex my-float))))
 713    (setf *psf-big* psf-big)
 714    (destructuring-bind (z y x)
 715        dims
 716      (let ((ox (- (floor x 2) (floor xx 2)))
 717            (oy (- (floor y 2) (floor yy 2)))
 718            (oz (- (floor z 2) (floor zz 2))))
 719        (do-box (k j i 0 zz 0 yy 0 xx)
 720          (setf (aref psf-big (+ oz k) (+ oy j) (+ ox i))
 721                (aref psf k j i)))))
 722    (save-stack-ub8 "/home/martin/tmp/psf-detect-big" (normalize-ub8 psf-big))
 723    (sb-ext:gc :full t)
 724    (defparameter *camera-volume* (convolve3-circ *slice-x-psf-times-spheres*
 725                                                  (fftshift3 psf-big)))
 726    (save-stack-ub8 "/home/martin/tmp/camera-volume"
 727                    (normalize-ub8 *camera-volume*))
 728    (sb-ext:gc :full t)))
 729
 730
 731
 732 #+nil ;; check convolution
 733 (time
 734  (let ((a (make-array (list 64 64 64)
 735                       :element-type '(complex my-float)))
 736        (b (psf:intensity-psf 64 64 64 20d0 20d0) #+nil (make-array (list 64 64 64)
 737                       :element-type '(complex my-float))))
 738    (setf (aref a 12 12 12) (complex 255d0))
 739 #+nil   (setf (aref b 0 0 0) (complex 255d0))
 740    (save-stack-ub8 "/home/martin/tmp/conv-test" (normalize-ub8 (convolve3-circ a (fftshift3 b))))))
 741
 742
 743 #+nil ;; output the xz cross section centered on a sphere in the middle
 744 (let ((coord (aref *centers* 30)))
 745   (write-pgm (normalize-img (cross-section-xz *spheres* (* 5 (vec-i-z coord))))
 746              "/home/martin/tmp/cut-spheres.pgm")
 747   (write-pgm (normalize-img (cross-section-xz *psf-big*))
 748              "/home/martin/tmp/cut-psf-big.pgm")
 749   (write-pgm (normalize-img (cross-section-xz *slice-x-psf* (* 5 (vec-i-z coord))))
 750              "/home/martin/tmp/cut-slice-x-psf.pgm")
 751   (write-pgm (normalize-img (cross-section-xz (.* *spheres* *slice-x-psf*) (* 5 (vec-i-z coord))))
 752              "/home/martin/tmp/cut-exfluo.pgm"))
 753
 754 #+nil
 755 (save-stack-ub8 "/home/martin/tmp/spheres" (normalize-ub8 *spheres*))
 756
 757 #+nil
 758 (let ((sli (make-array (array-dimensions *spheres*)
 759                        :element-type '(complex my-float))))
 760   (destructuring-bind (z y x)
 761       (array-dimensions *spheres*)
 762     (do-box (k j i 0 z 0 y 0 x)
 763       (setf (aref sli k j i)
 764             (aref *spheres* k j i)))
 765     (let ((k (* 5 (vec-i-z (aref *centers* 30)))))
 766      (do-rectangle (j i 0 y 0 x)
 767        (setf (aref sli k j i)
 768              (* 10 (aref sli k j i)))))
 769     (defparameter *sli* sli))
 770   (save-stack-ub8 "/home/martin/tmp/spheres" (normalize-ub8 *sli*)))
 771
 772
 773
 774 #+nil
 775 (let ((a (sb-ext:array-storage-vector *slice*)))
 776   (reduce #'max (map 'vector #'abs a)))
 777
 778 #+nil
 779 (sb-ext:gc :full t)
 780
 781 #+nil ;; model with unscaled spheres
 782 (defparameter *blobs*
 783   (destructuring-bind (z y x)
 784       *dims*
 785     (draw-spheres 7d0 *centers* (* 5 z) y x)))
 786
 787
 788 #+nil ;; print ft of angular psf
 789 (time (let* ((radius .5d0)
 790         (x (- 1d0 radius))
 791         (psf (angular-psf x 0d0 radius)))
 792    (write-pgm (normalize-ub8 (cross-section-xz (fftshift3 (ft3 psf))))
 793               "/home/martin/tmp/cut-intens.pgm")))
 794 #+nil
 795 (let* ((intens (psf:intensity-psf 64 64 64 10d0 5d0))
 796        (k0 intens #+nil(fftshift3 (ft3 intens))))
 797   (save-stack "/home/martin/tmp/intens1" k0
 798               :function #'(lambda (x) (* 1d-5 (abs x))))
 799   (write-pgm (convert-img
 800               (cross-section-xz k0)
 801               #'(lambda (z) (* 1d-4 (abs z))))
 802              "/home/martin/tmp/intens1xz.pgm"))
 803
 804 #+nil
 805 (time
 806  (destructuring-bind (z y x)
 807      (array-dimensions *blobs*)
 808    (save-stack "/home/martin/tmp/blobs"
 809                *blobs*
 810                )))
 811
 812
 813 #+nil ;; clean up the garbage
 814 (sb-ext:gc :full t)
 815
 816
 817 #+nil
 818 (sb-vm:memory-usage :print-spaces t :count-spaces t)
 819 #+nil
 820 (sb-vm:memory-usage)
 821
 822 #+nil
 823 (sb-vm:instance-usage :dynamic)
 824
 825 #+nil
 826 (sb-vm:list-allocated-objects :dynamic)
 827
 828 #+nil
 829 (sb-vm:print-allocated-objects :dynamic)
 830
 831 #+nil
 832 (format t "~a~%" (sb-vm::type-breakdown :dynamic))
 833
 834
 835 ;;
 836 (defun split-by-char (char string)
 837     "Returns a list of substrings of string
 838 divided by ONE character CHAR each.
 839 Note: Two consecutive CHAR will be seen as
 840 if there were an empty string between them."
 841     (declare (character char)
 842              (string string))
 843     (loop for i = 0 then (1+ j)
 844        as j = (position char string :start i)
 845        collect (subseq string i j)
 846        while j))
 847
 848 #+nil
 849 (split-by-char #\x "12x124x42")
 850
 851 (defun parse-raw-filename (fn)
 852   "Parses a filename like
 853 /home/martin/d0708/stacks/c-291x354x41x91_dx200_dz1000.raw and returns
 854 291 354 41 and 91 as multiple values."
 855   (declare (string fn)
 856            (values (or null fixnum) fixnum fixnum fixnum &optional))
 857   (let* ((p- (position #\- fn :from-end t))
 858          (part (subseq fn (1+ p-)))
 859          (p_ (position #\_ part))
 860          (sizes (subseq part 0 p_))
 861          (numlist (split-by-char #\x sizes)))
 862     (unless (eq 4 (length numlist))
 863       (error "didn't read 4 dimensions as expected."))
 864     (destructuring-bind (x y z time)
 865         (mapcar #'read-from-string numlist)
 866      (values x y z time))))
 867
 868
 869 #+nil
 870 (parse-raw-filename "/home/martin/d0708/stacks/c-291x354x41x91_dx200_dz1000.raw")
 871
 872 (defun read-raw-stack-video-frame (fn time)
 873   (declare (string fn)
 874            (fixnum time)
 875            (values (simple-array (unsigned-byte 8) 3) &optional))
 876   (multiple-value-bind (x y z maxtime)
 877       (parse-raw-filename fn)
 878       (unless (< time maxtime)
 879         (error "requested time ~d is to big (must be <~d!)" time maxtime))
 880       (let* ((vol (make-array (list z y x)
 881                               :element-type '(unsigned-byte 8)))
 882              (vol1 (sb-ext:array-storage-vector vol)))
 883         (with-open-file (s fn :direction :input
 884                            :element-type '(unsigned-byte 8))
 885           (file-position s (* x y z time))
 886           (read-sequence vol1 s))
 887         vol)))
 888 #+nil
 889 (time
 890  (let* ((fn "/home/martin/d0708/stacks/c-291x354x41x91_dx200_dz1000_2.raw")
 891         (ao (decimate-xy-ub8 5
 892                              (read-raw-stack-video-frame fn 0))))
 893    (destructuring-bind (z y x)
 894        (array-dimensions ao)
 895      (let* ((timestep 20)
 896             (o (loop for radius from 1 upto 10 collect
 897                     (let* ((oval (draw-sphere-ub8 (* 1d0 radius) z y x))
 898                            (volume (count-non-zero-ub8 oval)))
 899                       (list radius volume
 900                             (ft3 (convert3-ub8/cdf-complex oval)))))))
 901        (let* ((ao (decimate-xy-ub8 5
 902                                    (read-raw-stack-video-frame fn timestep)))
 903               (a (convert3-ub8/cdf-complex ao))
 904               (ka (ft3 a)))
 905          (loop for i in o do
 906               (destructuring-bind (radius volume oval)
 907                   i
 908                 (let* ((dir (format nil "/home/martin/tmp/o~d" radius))
 909                        (conv (fftshift3 (ift3 (.* ka oval))))
 910                        (conv-df (convert3-cdf/df-realpart conv)))
 911                  (save-stack-ub8 dir
 912                                  (normalize-ub8-df/ub8-realpart conv-df))
 913                  (with-open-file (s (format nil "~a/maxima" dir)
 914                                     :if-exists :supersede
 915                                    :direction :output)
 916                    (loop for el in (find-maxima3-df conv-df) do
 917                         (destructuring-bind (height pos)
 918                             el
 919                           (format s "~f ~d ~a~%"
 920                                   (/ height volume)
 921                                   volume
 922                                   (v*-i
 923                                    (map 'vec-i #'(lambda (x) (floor x 2)) pos)
 924                                    2))))))))
 925         nil)))))
 926 #+nil
 927 (sb-ext:gc :full t)
 928 #+nil
 929 (save-stack-ub8 "/home/martin/tmp/conv" conv)
 930
 931 #+nil
 932 (find-maxima3 conv)
 933
 934 #+nil
 935 (destructuring-bind (z y x)
 936     (array-dimensions *a*)
 937   (let ((b (draw-ovals 12d0 (find-maxima3 conv) (ensure-even z) (ensure-even y) (ensure-even x))))
 938     (write-pgm (convert-img (cross-section-xz b 42))
 939                "/home/martin/tmp/time0.pgm")))
 940
 941 (defun find-maxima3-df (conv)
 942   (declare ((simple-array my-float 3) conv)
 943            (values (simple-array * 1) &optional))
 944  (destructuring-bind (z y x)
 945      (array-dimensions conv)
 946    (let ((centers nil #+nil(make-array 0 :element-type 'vec-i
 947                               :initial-element (make-vec-i)
 948                               :adjustable t
 949                               :fill-pointer t)))
 950      (do-box (k j i 6 (- z 3) 1 (1- y) 1 (1- x))
 951        (let ((v (aref conv k j i)))
 952          (when (and (< (aref conv k j (1- i)) v)
 953                     (< (aref conv k j (1+ i)) v)
 954                     (< (aref conv k (1- j) i) v)
 955                     (< (aref conv k (1+ j) i) v)
 956                     (< (aref conv (1- k) j i) v)
 957                     (< (aref conv (1+ k) j i) v))
 958            ;; this is TOO slow
 959            #+nil(setf centers (append centers
 960                                  (list (list v (make-vec-i :z k :y j :x i)))))
 961            ;; this is faster
 962            #+nil(setf centers (nconc centers
 963                                 (list (list v (make-vec-i :z k :y j :x i)))))
 964            ;; I think push is the right thing to do
 965            (push (list v (make-vec-i :z k :y j :x i)) centers)
 966            #+nil(vector-push-extend
 967             (make-vec-i :z k :y j :x i)
 968             centers))))
 969 ;;    (nreverse centers)
 970      (coerce centers 'simple-vector) ;; I saw this in raylisps load-obj
 971    #+nil  (make-array (length centers)
 972                  :element-type 'vec-i
 973                  :initial-contents centers))))
 974