run.lisp

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