run.lisp

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