angular-illumination.lisp

   1 #.(require :gui)
   2 (in-package :run)
   3
   4 #+nil
   5 (let* ((k 25)
   6        (nuc (first (get-visible-nuclei k)))
   7        (vol (get-lcos-volume k nuc)))
   8   (format t "~a~%" `(nuc ,nuc))
   9   (write-section "/home/martin/tmp/angular-0lcos-cut.pgm" vol)
  10   (save-stack-ub8 "/home/martin/tmp/angular-0lcos" (normalize-3-csf/ub8-realpart vol)))
  11
  12 ;; FIXME: are these coordinates in mm or relative positions for a bfp-radius of 1?
  13 ;; I think the latter, but I'm not sure.
  14 #+nil
  15 (time
  16  (let*((dx .2d0)
  17        (dz 1d0)
  18        (r 100)
  19        (z (array-dimension *spheres* 0))
  20        (psf (angular-psf :x r :y r :z (* 2 z)
  21                 :pixel-size-x dx :pixel-size-z dz
  22                 :window-radius *bfp-window-radius*
  23                 :window-x -.714d0
  24                 :window-y .16d0
  25                 :initialize t
  26                 :debug t
  27                 :integrand-evaluations 400)))
  28    (save-stack-ub8 "/home/martin/tmp/psf" (normalize3-cdf/ub8-realpart psf))
  29    nil))
  30
  31
  32 ;; calculate the excitation one nucleus
  33 #+bla
  34 (defun calc-light-field (k nucleus)
  35   (declare (fixnum k nucleus))
  36   (let* ((result nil)
  37          (lcos (get-lcos-volume k nucleus))
  38          (bfp-pos (find-optimal-bfp-window-center nucleus))
  39          (psf (let ((dx .2d0)
  40                     (dz 1d0)
  41                     (r 100)
  42                     (z (array-dimension *spheres* 0)))
  43                 (multiple-value-bind (h ddx ddz)
  44                     (angular-intensity-psf-minimal-resolution
  45                                  :x-um (* r dx) :y-um (* r dx) :z-um (* 2 z dz)
  46                                  :window-radius *bfp-window-radius*
  47                                  :window-x (vec2-x bfp-pos)
  48                                  :window-y (vec2-y bfp-pos)
  49                                  :initialize t
  50                                  :debug t
  51                                  :integrand-evaluations 1000)
  52                   (resample-3-cdf h ddx ddx ddz dx dx dz)))))
  53     (format t "~a~%" `(bfp-pos ,bfp-pos))
  54     (write-section (format nil "/home/martin/tmp/angular-1expsf-cut-~3,'0d.pgm" nucleus) psf)
  55
  56     (multiple-value-bind (conv conv-start)
  57         (convolve-nocrop lcos psf)
  58       ;; light distribution in sample
  59       (defparameter *angular-light-field* conv)
  60       (defparameter *angular-light-field-start* conv-start)
  61       (write-section (format nil "/home/martin/tmp/angular-2light-cut-~3,'0d.pgm" nucleus)
  62                      conv
  63                      (vec-i-y (aref *centers* nucleus)))
  64       (save-stack-ub8 (format nil "/home/martin/tmp/angular-2light-~3,'0d/" nucleus)
  65                       (normalize-3-cdf/ub8-realpart conv))
  66       ;; multiply fluorophore concentration with light distribution
  67       (let ((excite (.* conv *spheres* conv-start)))
  68         (defparameter *excite* excite)
  69         (write-section (format nil "/home/martin/tmp/angular-3excite-cut-~3,'0d.pgm" nucleus)
  70                        excite
  71                        (vec-i-y (aref *centers* nucleus)))
  72         (save-stack-ub8 (format nil "/home/martin/tmp/angular-3excite-~3,'0d/" nucleus)
  73                         (normalize-3-cdf/ub8-realpart excite))
  74         (destructuring-bind (z y x)
  75             (array-dimensions excite)
  76           (declare (ignorable z))
  77           (let* ((in-focus (extract-bbox-3-cdf excite
  78                                               (make-bbox :start (v 0d0 0d0 (* 1d0 k))
  79                                                          :end (v (* 1d0 (1- x))
  80                                                                  (* 1d0 (1- y))
  81                                                                  (* 1d0 k))))))
  82             (save-stack-ub8 "/home/martin/tmp/angular-4in-focus/"
  83                             (normalize-3-cdf/ub8-realpart in-focus))
  84             (let*((mplane (mean (convert-3-cdf/df-realpart in-focus)))
  85                   (mvol (mean (convert-3-cdf/df-realpart excite)))
  86                   (gamma (/ mplane mvol)))
  87               (push (list mplane mvol gamma) result)
  88               (debug-out mplane mvol gamma)
  89               (format t "plane-result ~f ~f ~f~%" mplane mvol gamma))
  90             result))))))
  91
  92 #+nil
  93 (time
  94  (progn
  95   (with-open-file (*standard-output* "/home/martin/tmp/angular-stack.log"
  96                                      :direction :output
  97                                      :if-exists :supersede
  98                                      :if-does-not-exist :create)
  99     (let ((result nil))
 100       (dotimes (k (array-dimension *spheres* 0))
 101         (let* ((nucs (get-visible-nuclei k)))
 102           (loop for nuc in nucs do
 103                (format t "~a~%" (list 'doing k nucs))
 104                (push (list k nuc (calc-light-field k nuc)) result))))
 105       (defparameter *scan-result* result)))
 106   (with-open-file (s "/home/martin/tmp/angular-stack-struc.lisp"
 107                      :direction :output
 108                      :if-exists :supersede
 109                      :if-does-not-exist :create)
 110     (write *scan-result* :stream s))))
 111
 112
 113
 114
 115
 116 #+nil
 117 (let ((vol (make-array dims :element-type '(unsigned-byte 8))))
 118   (loop for i in '(4.0d-3 -.2d-3) do
 119    (draw-ray-into-vol i 0d0 .99d0 .0d0 vol)
 120    #+nil(draw-ray-into-vol i 0d0 -.99d0 .0d0 vol)
 121    (draw-ray-into-vol i 0d0 0d0 .99d0 vol)
 122    (draw-ray-into-vol i 0d0 0d0 -.99d0 vol))
 123
 124   (save-stack-ub8 "/home/martin/tmp/line"
 125                   vol))
 126
 127 #+nil
 128 frontend::merit-function
 129
 130 #+nil
 131 (time
 132  (let* ((n 100)
 133         (a (make-array (list n n) :element-type '(unsigned-byte 8)))
 134         (nn (length *spheres-c-r*))
 135         (mosaicx (ceiling (sqrt nn)))
 136         (mosaic (make-array (list (* n mosaicx) (* n mosaicx))
 137                             :element-type '(unsigned-byte 8))))
 138    (with-open-file (*standard-output* "/dev/shm/a"
 139                                       :direction :output
 140                                       :if-exists :supersede)
 141      (dotimes (*nucleus-index* 1 nn)
 142        (dotimes (i 10)
 143          (tagbody again
 144             (multiple-value-bind (min point)
 145                 (simplex-anneal:anneal (simplex-anneal:make-simplex
 146                                         (make-vec2 :x -1d0 :y -1d0) 1d0)
 147                                        #'merit-function
 148                                        ;; set temperature bigger than the
 149                                        ;; maxima in the bfp but smaller
 150                                        ;; than border-value
 151                                        :start-temperature 2.4d0
 152                                        :eps/m .02d0
 153                                        :itmax 1000
 154                                        :ftol 1d-3
 155                                        :params )
 156               (unless (<= min 100d0)
 157                 (go again))
 158               (let* ((x (aref point 0))
 159                      (y (aref point 1))
 160                      (ix (floor (* n (+ x 1)) 2))
 161                      (iy (floor (* n (+ y 1)) 2))
 162                      (mx (mod *nucleus-index* mosaicx))
 163                      (my (floor *nucleus-index* mosaicx)))
 164                 (incf (aref mosaic (+ (* n my) iy) (+ (* n mx) ix)))
 165                 (format t "min ~a~%" (list min ix iy))))))))
 166    (write-pgm "/home/martin/tmp/scan-mosaic-max.pgm" mosaic)))
 167
 168
 169