missing quote in coerce
[woropt.git] / run.lisp
blobc3be4fd2cdcf2d9f3718d2efd05b0454101d39bd
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))
11 (defpackage :run
12 (:use :cl :vector :vol :raytrace
13 :bresenham))
14 (in-package :run)
16 #+nil
17 (init-ft)
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))))
27 (defstuff)
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))
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))
52 ;; store the fluorophore concentration
53 (save-stack-ub8 "/home/martin/tmp/spheres"
54 (normalize3-cdf/ub8-realpart *spheres*)))
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)))
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)))
98 ;; store light distribution
99 (save-stack-ub8 "/home/martin/tmp/conv-l"
100 (normalize3-cdf/ub8-realpart *conv-l*))
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))))))
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)))))
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)))
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))))))
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)))))
166 #+nil
167 (time
168 (progn
169 (init-model)
170 (init-psf))) ;; 7.5s
172 #+nil
173 (time
174 (clem)) ;; 8s, result: 6.5
176 #+nil
177 (time
178 (widefield)) ;; 5.7s result: 1.8
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
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)
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)))))))
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)
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)))))
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))))))
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))))))
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))))))
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*))
373 ;; (- z k 1) (- y j 1) (- x i 1)
377 #+nil ;; model with isotropic pixels
378 (save-stack-ub "/home/martin/tmp/iso/iso"
379 (convert-vol (draw-scaled-spheres 2d0 *centers*)))
381 #+nil
382 (save-stack-ub "/home/martin/tmp/blobs" *blobs*)
384 #+nil
385 (save-scaled-stack 100 100 "/home/martin/tmp/sca" *blobs*)
387 #+nil ;; find maximum
388 (reduce #'max (map 'vector #'abs (sb-ext:array-storage-vector *stack*)))
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))))
397 #+nil
398 (save-stack "/home/martin/tmp/stack" *stack*)
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 "~%"))))))
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))))))
427 ;; worm stack 0.198 um in X and Y and 1 um in Z
429 #+nil
430 (array-dimensions *blobs*)
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)))))
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))))))
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))
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)))))))
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"
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"))))))
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)))
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)))
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)
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)))
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))
583 (defparameter *bfp-circ-radius* .3d0)
584 (defparameter *bfp-circ-center-x* .4d0 #+nil (- .999d0 *bfp-circ-radius*))
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))
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)))
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")
634 #+nil
635 (save-stack-ub8 "/home/martin/tmp/psf" (normalize-ub8 *psf*))
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
663 :shift-z (- (vec-i-z coord)
664 (floor z 2)))))
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)))
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*)))
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)))
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))))))
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"))
736 #+nil
737 (save-stack-ub8 "/home/martin/tmp/spheres" (normalize-ub8 *spheres*))
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*)))
756 #+nil
757 (let ((a (sb-ext:array-storage-vector *slice*)))
758 (reduce #'max (map 'vector #'abs a)))
760 #+nil
761 (sb-ext:gc :full t)
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)))
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"))
786 #+nil
787 (time
788 (destructuring-bind (z y x)
789 (array-dimensions *blobs*)
790 (save-stack "/home/martin/tmp/blobs"
791 *blobs*
795 #+nil ;; clean up the garbage
796 (sb-ext:gc :full t)
799 #+nil
800 (sb-vm:memory-usage :print-spaces t :count-spaces t)
801 #+nil
802 (sb-vm:memory-usage)
804 #+nil
805 (sb-vm:instance-usage :dynamic)
807 #+nil
808 (sb-vm:list-allocated-objects :dynamic)
810 #+nil
811 (sb-vm:print-allocated-objects :dynamic)
813 #+nil
814 (format t "~a~%" (sb-vm::type-breakdown :dynamic))
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))
830 #+nil
831 (split-by-char #\x "12x124x42")
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))))
851 #+nil
852 (parse-raw-filename "/home/martin/d0708/stacks/c-291x354x41x91_dx200_dz1000.raw")
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)
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)
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)
913 #+nil
914 (find-maxima3 conv)
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")))
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))))