external/clem/src/typed-ops/defmatrix-convolve.lisp

   1 ;;;
   2 ;;; file: defmatrix-conolve.cl
   3 ;;; author: cyrus harmon
   4 ;;;
   5
   6 (in-package :clem)
   7
   8 (defmacro def-matrix-convolve (type-1 type-2 accumulator-element-type dest-type)
   9   (let ((element-type-1 (element-type (find-class `,type-1)))
  10         (element-type-2 (element-type (find-class `,type-2)))
  11         (dest-element-type (element-type (find-class `,dest-type))))
  12     (let ((one (coerce 1 accumulator-element-type))
  13           (zero (coerce 0 accumulator-element-type)))
  14       `(progn
  15          (defmethod %discrete-convolve ((u ,type-1)
  16                                         (v ,type-2)
  17                                         (z ,dest-type)
  18                                         &key
  19                                         (norm-v t))
  20            (declare (optimize (speed 3) (safety 0) (space 0)))
  21            ;; ur, uc, vr, vc are the number of rows and columns in u and v
  22            (let ((ur (rows u)) (uc (cols u)) (vr (rows v)) (vc (cols v)))
  23              (declare (type fixnum ur uc vr vc))
  24              ;; need a new matrix z to hold the values of the convolved matrix
  25              ;; dim z should be dim u + dim v - 1
  26              (let ((zr (+ ur vr (- 1)))
  27                    (zc (+ uc vc (- 1))))
  28                ;; now that we really have z we should get the dimensions and check
  29                ;; instead!
  30                (declare (type fixnum zr zc))
  31                (let ((uval (matrix-vals u))
  32                      (vval (matrix-vals v))
  33                      (zval (matrix-vals z))
  34                      (vsum (sum v)))
  35                  (declare (type (simple-array ,element-type-1 *) uval)
  36                           (type (simple-array ,element-type-2 *) vval)
  37                           (type (simple-array ,dest-element-type *) zval)
  38                           (type ,accumulator-element-type vsum))
  39                  (dotimes (i zr)
  40                    (declare (type fixnum i))
  41                    (let ((ustartr (max 0 (the fixnum (- i vr -1))))
  42                          (uendr (min (- ur 1) i))
  43                          (vstartr (- vr (max (the fixnum (- vr i)) 1)))
  44                          (vendr (- vr (min (- zr i) vr))))
  45                      (declare (type fixnum ustartr uendr vstartr vendr))
  46                      (dotimes (j zc)
  47                        (declare (type fixnum j))
  48                        (let ((ustartc (max 0 (the fixnum (- j vc -1))))
  49                              (uendc (min (- uc 1) j))
  50                              (vstartc (- vc (max (the fixnum (- vc j)) 1)))
  51                              (vendc (- vc (min (- zc j) vc)))
  52                              (acc ,zero))
  53                          (declare (type fixnum ustartr uendr vstartr vendr)
  54                                   (type ,accumulator-element-type acc))
  55                          (let ((normval
  56                                 (if (and norm-v
  57                                          (or (not (= vendr vendc 0))
  58                                              (< vstartr (the fixnum (1- vr)))
  59                                              (< vstartc (the fixnum (1- vc)))))
  60                                     (let ((rsum
  61                                            (%%sum-range v vendr vstartr vendc vstartc
  62                                                         ,element-type-2
  63                                                         ,(car (slot-value
  64                                                                `,(find-class `,type-2)
  65                                                                'clem::accumulator-type)))))
  66                                       (declare (type ,accumulator-element-type rsum))
  67                                       (if (not (= rsum ,zero))
  68                                           (/ vsum rsum)
  69                                           ,zero))
  70                                     ,one)))
  71                            (declare (type ,accumulator-element-type normval))
  72                            (do ((urow ustartr (the fixnum (1+ urow)))
  73                                 (vrow vstartr (the fixnum (1- vrow))))
  74                                ((> urow uendr))
  75                              (declare (type fixnum urow vrow))
  76                              (do ((ucol ustartc (the fixnum (1+ ucol)))
  77                                   (vcol vstartc (the fixnum (1- vcol))))
  78                                  ((> ucol uendc))
  79                                (declare (type fixnum ucol vcol))
  80                                (let ((uv (aref uval urow ucol))
  81                                      (vv (aref vval vrow vcol)))
  82                                  (declare (type ,element-type-1 uv)
  83                                           (type ,element-type-2 vv))
  84                                  (incf acc (* uv vv))
  85                                  )))
  86
  87                            (setf (aref zval i j)
  88                                  (maybe-truncate
  89                                   (* acc normval)
  90                                   ,accumulator-element-type
  91                                   ,dest-element-type)))))))
  92                  z))))
  93
  94          (defmethod discrete-convolve ((u ,type-1) (v ,type-2)
  95                                        &key (truncate nil) (norm-v t)
  96                                        (matrix-class ',dest-type))
  97            (declare (ignore truncate))
  98            ;; ur, uc, vr, vc are the number of rows and columns in u and v
  99            (let ((ur (rows u)) (uc (cols u)) (vr (rows v)) (vc (cols v)))
 100              (let ((zr (+ ur vr (- 1)))
 101                    (zc (+ uc vc (- 1))))
 102                (unless matrix-class
 103                  (setf matrix-class (type-of u)))
 104                (let ((z (make-instance matrix-class :rows zr :cols zc)))
 105                  (%discrete-convolve u v z :norm-v norm-v)))))
 106
 107          (defmethod %separable-discrete-convolve ((m ,type-1) (h1 ,type-2) (h2 ,type-2)
 108                                                  (z1 matrix) (z2 matrix)
 109                                                  &key (truncate nil) (norm-v nil)
 110                                                  (matrix-class (class-of m)))
 111            (declare (ignore truncate matrix-class))
 112            (let ((mr (rows m)) (mc (cols m)) (hr (rows h1)) (hc (cols h2)))
 113              (declare (ignore mr mc hr hc))
 114              (%discrete-convolve m h2 z1 :norm-v norm-v)
 115              (%discrete-convolve z1 h1 z2 :norm-v norm-v)))
 116
 117
 118          (defmethod separable-discrete-convolve ((m ,type-1) (h1 ,type-2) (h2 ,type-2)
 119                                                  &key (truncate nil) (norm-v nil)
 120                                                  (matrix-class (class-of m)))
 121            (declare (ignore truncate))
 122            (let ((mr (rows m)) (mc (cols m)) (hr (rows h1)) (hc (cols h2)))
 123              (let ((z1r mr) (z1c (+ mc hc -1))
 124                    (z2r (+ mr hr -1)) (z2c (+ mc hc -1)))
 125                (let ((z1 (make-instance matrix-class :rows z1r :cols z1c))
 126                      (z2 (make-instance matrix-class :rows z2r :cols z2c)))
 127                  (%separable-discrete-convolve m h1 h2 z1 z2 :norm-v norm-v)))))))))
 128
 129 (macrolet ((frob (type-1 type-2 type-3 type-4)
 130              `(progn
 131                 (def-matrix-convolve ,type-1 ,type-2 ,type-3 ,type-4))))
 132   (frob double-float-matrix double-float-matrix double-float double-float-matrix)
 133   (frob single-float-matrix single-float-matrix single-float single-float-matrix)
 134   (frob ub8-matrix ub8-matrix (unsigned-byte 8) ub8-matrix)
 135   (frob ub8-matrix double-float-matrix double-float ub8-matrix)
 136   (frob ub8-matrix single-float-matrix single-float ub8-matrix))
 137