vol-operators.lisp

   1 (in-package :vol)
   2
   3 (defun .*2 (vola volb)
   4   (declare ((simple-array (complex my-float) 2) vola volb)
   5            (values (simple-array (complex my-float) 2) &optional))
   6   (let ((result (make-array (array-dimensions vola)
   7                             :element-type (array-element-type vola))))
   8    (destructuring-bind (y x)
   9        (array-dimensions vola)
  10      (do-rectangle (j i 0 y 0 x)
  11        (setf (aref result j i)
  12              (* (aref vola j i)
  13                 (aref volb j i)))))
  14    result))
  15
  16 (defun .+2 (vola volb)
  17   (declare ((simple-array (complex my-float) 2) vola volb)
  18            (values (simple-array (complex my-float) 2) &optional))
  19   (let ((result (make-array (array-dimensions vola)
  20                             :element-type (array-element-type vola))))
  21    (destructuring-bind (y x)
  22        (array-dimensions vola)
  23      (do-rectangle (j i 0 y 0 x)
  24        (setf (aref result j i)
  25              (+ (aref vola j i)
  26                 (aref volb j i)))))
  27    result))
  28
  29 (defun .* (vola volb &optional volb-start)
  30   "Elementwise multiplication of VOLA and VOLB. Both volumes must have
  31 the same dimensions or VOLB must be smaller in all dimensions. In the
  32 latter case a vec-i has to be supplied in VOLB-START to define the
  33 relative position of VOLB inside VOLA."
  34   (declare ((simple-array (complex my-float) 3) vola volb)
  35            ((or null vec-i) volb-start)
  36            (values (simple-array (complex my-float) 3) &optional))
  37   (let ((result (make-array (array-dimensions volb)
  38                             :element-type '(complex my-float))))
  39    (destructuring-bind (z y x)
  40        (array-dimensions vola)
  41      (destructuring-bind (zz yy xx)
  42          (array-dimensions volb)
  43        (if volb-start
  44            ;; fill the result with volb multiplied by the
  45            ;; corresponding values from the bigger vola
  46            (let ((sx (vec-i-x volb-start))
  47                  (sy (vec-i-y volb-start))
  48                  (sz (vec-i-z volb-start)))
  49              (unless (and (<= zz (+ z sz))
  50                           (<= yy (+ y sy))
  51                           (<= xx (+ x sx)))
  52                (error "VOLB isn't contained in VOLA when shifted by VOLB-START. ~a"
  53                       (list zz (+ z sz))))
  54              (do-box (k j i 0 zz 0 yy 0 xx)
  55                (setf (aref result k j i)
  56                      (* (aref volb k j i)
  57                         (aref vola (+ k sz) (+ j sy) (+ i sx))))))
  58            (progn
  59              (unless (and (= z zz) (= y yy) (= x xx))
  60                (error "volumes don't have the same size, maybe you can supply a start vector."))
  61              (do-box (k j i 0 z 0 y 0 x)
  62                (setf (aref result k j i)
  63                      (* (aref vola k j i)
  64                         (aref volb k j i))))))))
  65    result))
  66
  67 (defun .+ (vola volb &optional (volb-start (make-vec-i)))
  68   (declare ((simple-array (complex my-float) 3) vola volb)
  69            (vec-i volb-start)
  70            (values (simple-array (complex my-float) 3) &optional))
  71   (let ((result (make-array (array-dimensions volb)
  72                             :element-type '(complex my-float))))
  73     (destructuring-bind (z y x)
  74         (array-dimensions volb)
  75       (let ((sx (vec-i-x volb-start))
  76             (sy (vec-i-y volb-start))
  77             (sz (vec-i-z volb-start)))
  78         (do-box (k j i 0 z 0 y 0 x)
  79           (setf (aref result k j i)
  80                 (+ (aref vola (+ k sz) (+ j sy) (+ i sx))
  81                    (aref volb k j i))))))
  82     result))
  83
  84 (defun .- (vola volb &optional (volb-start (make-vec-i)))
  85   (declare ((simple-array (complex my-float) 3) vola volb)
  86            (vec-i volb-start)
  87            (values (simple-array (complex my-float) 3) &optional))
  88   (let ((result (make-array (array-dimensions volb)
  89                             :element-type '(complex my-float))))
  90     (destructuring-bind (z y x)
  91         (array-dimensions volb)
  92       (let ((sx (vec-i-x volb-start))
  93             (sy (vec-i-y volb-start))
  94             (sz (vec-i-z volb-start)))
  95         (do-box (k j i 0 z 0 y 0 x)
  96           (setf (aref result k j i)
  97                 (- (aref vola (+ k sz) (+ j sy) (+ i sx))
  98                    (aref volb k j i))))))
  99     result))
 100
 101 (declaim (ftype (function (my-float (simple-array (complex my-float) 3))
 102                           (values (simple-array (complex my-float) 3) &optional))
 103                 s*))
 104 (defun s* (s vol)
 105   (let* ((a (sb-ext:array-storage-vector vol))
 106          (n (length a)))
 107     (dotimes (i n)
 108       (setf (aref a i) (* s (aref a i)))))
 109   vol)
 110
 111 (defun s*2 (s vol)
 112   (declare (my-float s)
 113            ((simple-array (complex my-float) 2) vol)
 114            (values (simple-array (complex my-float) 2) &optional))
 115   (let* ((a (sb-ext:array-storage-vector vol))
 116          (n (length a)))
 117     (dotimes (i n)
 118       (setf (aref a i) (* s (aref a i)))))
 119   vol)
 120
 121
 122 (defun mean-realpart (a)
 123   "Calculate the average value over all the samples in volume A."
 124   (declare ((simple-array (complex my-float) *) a)
 125            (values my-float &optional))
 126   (let* ((a1 (sb-ext:array-storage-vector a))
 127          (sum zero)
 128          (n (length a1)))
 129     (dotimes (i n)
 130       (incf sum (realpart (aref a1 i))))
 131     (/ sum n)))