vol-convolve.lisp

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