vol-convolve.lisp

   1 (defun convolve2-circ (vola volb)
   2   (declare ((simple-array (complex my-float) 2) vola volb)
   3            (values (simple-array (complex my-float) 2) &optional))
   4   (let* ((da (array-dimensions vola))
   5          (db (array-dimensions volb))
   6          (compare-ab (map 'list #'(lambda (x y) (eq x y)) da db)))
   7     (when (some #'null compare-ab)
   8       (error "convolve3-circ expects both input arrays to have the same dimensions."))
   9     (ift2 (s*2 (* one (reduce #'* da)) (.*2 (ft2 vola) (ft2 volb))))))
  10
  11 (defun convolve3-circ (vola volb)
  12   (declare ((simple-array (complex my-float) 3) vola volb)
  13            (values (simple-array (complex my-float) 3) &optional))
  14   (let* ((da (array-dimensions vola))
  15          (db (array-dimensions volb))
  16          (compare-ab (map 'list #'(lambda (x y) (eq x y)) da db)))
  17     (when (some #'null compare-ab)
  18       (error "convolve3-circ expects both input arrays to have the same dimensions.")))
  19   (ift3 (.* (ft3 vola) (ft3 volb))))
  20
  21
  22 (defun front (i) ;; extra size needed to accommodate kernel overlap
  23                  ;; there is a difference between even and odd kernels
  24   (declare (fixnum i)
  25            (values fixnum &optional))
  26   (max 0
  27        (if (evenp i)
  28            (floor i 2)
  29            (1- (floor i 2)))))
  30
  31 ;; volb is the kernel
  32 (defun convolve3-nocrop (vola volb)
  33   "Convolve VOLA with VOLB. We consider VOLB as the convolution
  34 kernel. Returns (values result vec). RESULT is an arrays that is as
  35 big as necessary to accommodate the convolution and VEC contains the
  36 relative coordinates to find the original sample positions of array
  37 VOLA in RESULT."
  38     (declare ((simple-array (complex my-float) 3) vola volb)
  39            (values (simple-array (complex my-float) 3) vec-i &optional))
  40   (destructuring-bind (za ya xa)
  41       (array-dimensions vola)
  42     (destructuring-bind (zb yb xb)
  43         (array-dimensions volb)
  44       (let* ((biga (make-array (list (+ za zb)
  45                                      (+ ya yb)
  46                                      (+ xa xb))
  47                                :element-type '(complex my-float)))
  48              (bigb (make-array (array-dimensions biga)
  49                                :element-type '(complex my-float)))
  50              (fzb (front zb))
  51              (fyb (front yb))
  52              (fxb (front xb))
  53              (fza (front za))
  54              (fya (front ya))
  55              (fxa (front xa))
  56              (start (make-vec-i :x fxb :y fyb :z fzb)))
  57         (do-box (k j i 0 za 0 ya 0 xa)
  58           (setf (aref biga (+ k fzb) (+ j fyb) (+ i fxb))
  59                 (aref vola k j i)))
  60         (do-box (k j i 0 zb 0 yb 0 xb)
  61           (setf (aref bigb (+ k fza) (+ j fya) (+ i fxa))
  62                 (aref volb k j i)))
  63         (values (convolve3-circ biga (fftshift3 bigb))
  64                 start)))))
  65
  66 (defun convolve3 (vola volb)
  67   (destructuring-bind (za ya xa)
  68       (array-dimensions vola)
  69     (multiple-value-bind (conv start)
  70         (convolve3-nocrop vola volb)
  71       (let* ((result (make-array (array-dimensions vola)
  72                                  :element-type '(complex my-float)))
  73              (oz (vec-i-z start))
  74              (oy (vec-i-y start))
  75              (ox (vec-i-x start)))
  76         (do-box (k j i 0 za 0 ya 0 xa)
  77           (setf (aref result k j i)
  78                (aref conv (+ k oz) (+ j oy) (+ i ox))))
  79         result))))
  80
  81 #+nil
  82 (let ((a (make-array (list 100 200 300)
  83                      :element-type '(complex my-float)))
  84       (b (make-array (list 10 200 30)
  85                      :element-type '(complex my-float))))
  86   (convolve3 a b)
  87   nil)