vol/stack-timeseries.lisp

   1 (defun split-by-char (char string)
   2     "Returns a list of substrings of string
   3 divided by ONE character CHAR each.
   4 Note: Two consecutive CHAR will be seen as
   5 if there were an empty string between them."
   6     (declare (character char)
   7              (string string))
   8     (loop for i = 0 then (1+ j)
   9        as j = (position char string :start i)
  10        collect (subseq string i j)
  11        while j))
  12
  13 #+nil
  14 (split-by-char #\x "12x124x42")
  15
  16 (defun parse-raw-filename (fn)
  17   "Parses a filename like
  18 /home/martin/d0708/stacks/c-291x354x41x91_dx200_dz1000.raw and returns
  19 291 354 41 and 91 as multiple values."
  20   (declare (string fn)
  21            (values (or null fixnum) fixnum fixnum fixnum &optional))
  22   (let* ((p- (position #\- fn :from-end t))
  23          (part (subseq fn (1+ p-)))
  24          (p_ (position #\_ part))
  25          (sizes (subseq part 0 p_))
  26          (numlist (split-by-char #\x sizes)))
  27     (unless (eq 4 (length numlist))
  28       (error "didn't read 4 dimensions as expected."))
  29     (destructuring-bind (x y z time)
  30         (mapcar #'read-from-string numlist)
  31      (values x y z time))))
  32
  33
  34 #+nil
  35 (parse-raw-filename "/home/martin/d0708/stacks/c-291x354x41x91_dx200_dz1000.raw")
  36
  37 (defun read-raw-stack-video-frame (fn time)
  38   (declare (string fn)
  39            (fixnum time)
  40            (values (simple-array (unsigned-byte 8) 3) &optional))
  41   (multiple-value-bind (x y z maxtime)
  42       (parse-raw-filename fn)
  43       (unless (< time maxtime)
  44         (error "requested time ~d is to big (must be <~d!)" time maxtime))
  45       (let* ((vol (make-array (list z y x)
  46                               :element-type '(unsigned-byte 8)))
  47              (vol1 (sb-ext:array-storage-vector vol)))
  48         (with-open-file (s fn :direction :input
  49                            :element-type '(unsigned-byte 8))
  50           (file-position s (* x y z time))
  51           (read-sequence vol1 s))
  52         vol)))
  53 #+nil
  54 (time
  55  (let* ((fn "/home/martin/d0708/stacks/c-291x354x41x91_dx200_dz1000_2.raw")
  56         (ao (decimate-xy-ub8 5
  57                              (read-raw-stack-video-frame fn 0))))
  58    (destructuring-bind (z y x)
  59        (array-dimensions ao)
  60      (let* ((timestep 20)
  61             (o (loop for radius from 1 upto 10 collect
  62                     (let* ((oval (draw-sphere-ub8 (* 1d0 radius) z y x))
  63                            (volume (count-non-zero-ub8 oval)))
  64                       (list radius volume
  65                             (ft3 (convert3-ub8/cdf-complex oval)))))))
  66        (let* ((ao (decimate-xy-ub8 5
  67                                    (read-raw-stack-video-frame fn timestep)))
  68               (a (convert3-ub8/cdf-complex ao))
  69               (ka (ft3 a)))
  70          (loop for i in o do
  71               (destructuring-bind (radius volume oval)
  72                   i
  73                 (let* ((dir (format nil "/home/martin/tmp/o~d" radius))
  74                        (conv (fftshift3 (ift3 (.* ka oval))))
  75                        (conv-df (convert3-cdf/df-realpart conv)))
  76                  (save-stack-ub8 dir
  77                                  (normalize-ub8-df/ub8-realpart conv-df))
  78                  (with-open-file (s (format nil "~a/maxima" dir)
  79                                     :if-exists :supersede
  80                                    :direction :output)
  81                    (loop for el in (find-maxima3-df conv-df) do
  82                         (destructuring-bind (height pos)
  83                             el
  84                           (format s "~f ~d ~a~%"
  85                                   (/ height volume)
  86                                   volume
  87                                   (v*-i
  88                                    (map 'vec-i #'(lambda (x) (floor x 2)) pos)
  89                                    2))))))))
  90         nil)))))