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