vol/convolve.lisp

   1 ;; 2d and 3d convolutions
   2 (in-package :vol)
   3
   4 (def-generator (convolve-circ (rank type))
   5   `(defun ,name (vola volb)
   6      (declare ((simple-array ,long-type ,rank) vola volb)
   7             (values (simple-array ,long-type ,rank) &optional))
   8      (let* ((da (array-dimensions vola))
   9             (db (array-dimensions volb))
  10             (compare-ab (map 'list #'(lambda (x y) (eq x y)) da db)))
  11        (when (some #'null compare-ab)
  12          (error
  13           ,(format
  14             nil
  15             "~a expects both input arrays to have the same dimensions." name)))
  16        ;; fta and ftb are both already divided by n, the following ift
  17        ;; will introduce will not divide by n so i have to do that
  18        (ift (s* (/ ,(coerce 1 (ecase type (csf 'single-float) (cdf 'double-float)))
  19                   (array-total-size vola))
  20                (.* (ft vola) (ft (fftshift volb))))))))
  21
  22 (defmacro def-convolve-circ-functions (ranks types)
  23   (let* ((specifics nil)
  24          (cases nil)
  25          (name (format-symbol "convolve-circ")))
  26     (loop for rank in ranks do
  27          (loop for type in types do
  28               (let ((def-name (format-symbol "def-~a-rank-type" name))
  29                     (specific-name (format-symbol "~a-~a-~a" name rank type)))
  30                 (push `(,def-name ,rank ,type) specifics)
  31                 (push `((simple-array ,(get-long-type type) ,rank)
  32                         (,specific-name a b))
  33                       cases))))
  34     (store-new-function name)
  35     `(progn ,@specifics
  36             (defun ,name (a b)
  37               (etypecase a
  38                 ,@cases
  39                 (t (error "The given type can't be handled with a generic ~a function." ',name)))))))
  40
  41 (def-convolve-circ-functions (2 3) (csf cdf))
  42
  43
  44 (defun front (i) ;; extra size needed to accommodate kernel overlap
  45                  ;; there is a difference between even and odd kernels
  46   (declare (fixnum i)
  47            (values fixnum &optional))
  48   (max 0
  49        (if (evenp i)
  50            (floor i 2)
  51            (1- (floor i 2)))))
  52
  53 ;; volb is the kernel
  54 (def-generator (convolve-nocrop (rank type))
  55   `(defun ,name (vola volb)
  56    "Convolve VOLA with VOLB. We consider VOLB as the convolution
  57 kernel. Returns (values result vec). RESULT is an arrays that is as
  58 big as necessary to accommodate the convolution and VEC contains the
  59 relative coordinates to find the original sample positions of array
  60 VOLA in RESULT."
  61    (declare ((simple-array ,long-type ,rank) vola volb)
  62             (values (simple-array ,long-type ,rank) vec-i &optional))
  63    ,(ecase
  64      rank
  65      (2 `(destructuring-bind (ya xa) (array-dimensions vola)
  66            (destructuring-bind (yb xb) (array-dimensions volb)
  67              (let* ((biga (make-array (list (+ ya yb) (+ xa xb))
  68                                       :element-type ',long-type))
  69                     (bigb (make-array (array-dimensions biga)
  70                                       :element-type ',long-type))
  71                     (fyb (front yb))
  72                     (fxb (front xb))
  73                     (fya (front ya))
  74                     (fxa (front xa))
  75                     (start (make-vec-i :x fxb :y fyb)))
  76                (do-region ((j i) (ya xa))
  77                  (setf (aref biga (+ j fyb) (+ i fxb))
  78                        (aref vola j i)))
  79                (do-region ((j i) (yb xb))
  80                  (setf (aref bigb (+ j fya) (+ i fxa))
  81                        (aref volb j i)))
  82                (values (convolve-circ biga bigb) start)))))
  83      (3 `(destructuring-bind (za ya xa) (array-dimensions vola)
  84            (destructuring-bind (zb yb xb) (array-dimensions volb)
  85              (let* ((biga (make-array (list (+ za zb) (+ ya yb) (+ xa xb))
  86                                       :element-type ',long-type))
  87                     (bigb (make-array (array-dimensions biga)
  88                                       :element-type ',long-type))
  89                     (fzb (front zb))
  90                     (fyb (front yb))
  91                     (fxb (front xb))
  92                     (fza (front za))
  93                     (fya (front ya))
  94                     (fxa (front xa))
  95                     (start (make-vec-i :x fxb :y fyb :z fzb)))
  96                (do-region ((k j i) (za ya xa))
  97                  (setf (aref biga (+ k fzb) (+ j fyb) (+ i fxb))
  98                        (aref vola k j i)))
  99                (do-region ((k j i) (zb yb xb))
 100                  (setf (aref bigb (+ k fza) (+ j fya) (+ i fxa))
 101                        (aref volb k j i)))
 102                (values (convolve-circ biga bigb) start))))))))
 103
 104 (defmacro def-convolve-nocrop-functions (ranks types)
 105   (let* ((specifics nil)
 106          (cases nil)
 107          (name (format-symbol "convolve-nocrop")))
 108     (loop for rank in ranks do
 109          (loop for type in types do
 110               (let ((def-name (format-symbol "def-~a-rank-type" name))
 111                     (specific-name (format-symbol "~a-~a-~a" name rank type)))
 112                 (push `(,def-name ,rank ,type) specifics)
 113                 (push `((simple-array ,(get-long-type type) ,rank)
 114                         (,specific-name a b))
 115                       cases))))
 116     (store-new-function name)
 117     `(progn ,@specifics
 118             (defun ,name (a b)
 119                (etypecase a
 120                  ,@cases
 121                  (t (error "The given type can't be handled with a generic ~a function." ',name)))))))
 122
 123 (def-convolve-nocrop-functions (2 3) (csf cdf))
 124
 125 (defun convolve (vola kernel)
 126   (multiple-value-bind (conv start) (convolve-nocrop vola kernel)
 127     (let* ((s (convert-1-fix/df-mul start))
 128            (d (array-dimensions vola)) ;; z y x, y x or x
 129            (dims (ecase (length d)     ;; ensure 3 entries
 130                    (1 (push 0 d) (push 0 d))
 131                    (2 (push 0 d))
 132                    (3 d)))
 133            (b (convert-1-fix/df-mul
 134                (make-array 3 :element-type 'fixnum ;; z y x -> x y z
 135                            :initial-contents (reverse dims)))))
 136       (extract-bbox conv
 137                     (make-bbox :start s
 138                                :end (v- (v+ s b) (v 1 1 1)))))))
 139
 140 #+nil
 141 (let ((a (make-array (list 20 40 30)
 142                      :element-type '(complex single-float)))
 143       (b (make-array (list 10 20 30)
 144                      :element-type '(complex single-float))))
 145   (convolve a b)
 146   nil)