external/clem/src/matrixops.lisp

   1 ;;;
   2 ;;; file: matrixops.cl
   3 ;;; author: cyrus harmon
   4 ;;;
   5
   6 ;;;
   7 ;;; 2004-05-07 - This class contains matrix operations such as
   8 ;;;              gaussian-blur, gradient, etc...
   9 ;;;
  10 ;;;              Relies on the matrix package for matrix datatypes
  11 ;;;              and core functions such as discrete-convolve
  12 ;;;
  13
  14 (in-package :clem)
  15
  16
  17 ;;; discrete-convolve takes two matrices and returns
  18 ;;; a new matrix which is the convolution of the two matrices.
  19 ;;; To convolve, for each row,col of the matrix u, overlay matrix v
  20 ;;; the current cell and take the sum of the product of all of the
  21 ;;; u,v pairs. Note that v is rotated 180 degrees wrt u so that
  22 ;;; if we are calculating the value of (1,1) in the convolution of
  23 ;;; 3x3 matrices, the first number we would sum would be (0,0) x (2,2)
  24 ;;; not (0,0) x (0,0)
  25
  26 ;;; discrete-convolve takes two matrices and returns
  27 ;;; a new matrix which is the convolution of the two matrices.
  28 ;;; To convolve, for each row,col of the matrix u, overlay matrix v
  29 ;;; the current cell and take the sum of the product of all of the
  30 ;;; u,v pairs. Note that v is rotated 180 degrees wrt u so that
  31 ;;; if we are calculating the value of (1,1) in the convolution of
  32 ;;; 3x3 matrices, the first number we would sum would be (0,0) x (2,2)
  33 ;;; not (0,0) x (0,0)
  34
  35
  36 (defgeneric discrete-convolve (u v &key truncate norm-v matrix-class))
  37
  38 (defmethod discrete-convolve ((u matrix) (v matrix)
  39                               &key (truncate nil) (norm-v t)
  40                               (matrix-class nil))
  41   ;; (declare (optimize (speed 3) (safety 0) (space 0)))
  42   ;;; ur, uc, vr, vc are the number of rows and columns in u and v
  43   (print 'unspecialized-discrete-convolve!)
  44   (destructuring-bind (ur uc) (dim u)
  45     (declare (dynamic-extent ur uc) (fixnum ur uc))
  46     (destructuring-bind (vr vc) (dim v)
  47       (declare (fixnum vr vc) (dynamic-extent vr vc))
  48       ;;; need a new matrix z to hold the values of the convolved matrix
  49       ;;; dim z should be dim u + dim v - 1
  50       (let ((zr (+ ur vr (- 1)))
  51             (zc (+ uc vc (- 1))))
  52         (declare (fixnum zr zc) (dynamic-extent zr zc))
  53         (unless matrix-class
  54           (setf matrix-class (type-of u)))
  55         (let ((z (make-instance matrix-class :rows zr :cols zc))
  56               (vsum (sum v)))
  57           (dotimes (i zr)
  58             (let ((ustartr (max 0 (- i vr -1)))
  59                   (uendr (min (- ur 1) i))
  60                   (vstartr (- vr (max (- vr i) 1)))
  61                   (vendr (- vr (min (- zr i) vr))))
  62               (dotimes (j zc)
  63                 (let ((ustartc (max 0 (- j vc -1)))
  64                       (uendc (min (- uc 1) j))
  65                       (vstartc (- vc (max (- vc j) 1)))
  66                       (vendc (- vc (min (- zc j) vc)))
  67                       (acc 0))
  68                   (let ((normval (if (and norm-v (or (not (= vendr vendc 0))
  69                                                      (< vstartr (- vr 1))
  70                                                      (< vstartc (- vc 1))))
  71                                      (let ((rsum (sum-range v vendr vstartr vendc vstartc)))
  72                                        (if (not (= rsum 0))
  73                                            (/ vsum rsum)
  74                                          0))
  75                                    nil)))
  76                     (do ((urow ustartr (1+ urow))
  77                          (vrow vstartr (1- vrow)))
  78                         ((> urow uendr))
  79                       (declare (fixnum urow vrow))
  80                       (declare (dynamic-extent urow vrow))
  81                       (do ((ucol ustartc (1+ ucol))
  82                            (vcol vstartc (1- vcol)))
  83                           ((> ucol uendc))
  84                         (declare (fixnum ucol vcol))
  85                         (declare (dynamic-extent ucol vcol))
  86                         (let ((uv (val u urow ucol))
  87                               (vv (val v vrow vcol)))
  88                           (declare (dynamic-extent uv vv))
  89                           (incf acc (* uv vv))
  90                           )))
  91                     (if normval
  92                         (setf acc (fit-value (* acc normval) z)))
  93                     (if truncate
  94                         (set-val-fit z i j acc :truncate truncate)
  95                         (set-val z i j acc)))))))
  96           z)))))
  97
  98 (defun gaussian-kernel (k sigma)
  99   (let* ((d (1+ (* 2 k)))
 100          (a (make-instance 'double-float-matrix :rows d :cols d))
 101          (q (* 2 sigma sigma))
 102          (z (/ (* pi q))))
 103     (dotimes (i d)
 104       (dotimes (j d)
 105         (set-val a i j
 106                  (* (exp (- (/ (+ (* (- i k) (- i k))
 107                                   (* (- j k) (- j k)))
 108                                q)))
 109                     z))))
 110     (scalar-divide a (sum-range a 0 (- d 1) 0 (- d 1)))
 111     a))
 112
 113
 114 (defun gaussian-kernel-1d (k sigma)
 115   (let* ((d (1+ (* 2 k)))
 116          (a (make-instance 'double-float-matrix :rows d :cols 1))
 117          (q (* 2 sigma sigma))
 118          (z (/ (* pi q))))
 119     (dotimes (i d)
 120       (set-val a i 0
 121                (* (exp (- (/ (+ (* (- i k) (- i k))
 122                                 (* (- i k) (- i k)))
 123                              q)))
 124                   z)))
 125     (scalar-divide a (sum-range a 0 (- d 1) 0 0))
 126     a))
 127
 128 (defgeneric separable-discrete-convolve (m h1 h2 &key truncate norm-v))
 129
 130 (defmethod separable-discrete-convolve (m h1 h2 &key (truncate nil) (norm-v nil))
 131   (let* ((m1 (discrete-convolve m h1 :truncate truncate :norm-v norm-v))
 132          (m2 (discrete-convolve m1 h2 :truncate truncate :norm-v norm-v)))
 133     m2))
 134
 135 (defun gaussian-blur (m &key (k 2) (sigma 1) (truncate nil))
 136   (let* ((hr (gaussian-kernel-1d k sigma))
 137          (hc (transpose hr)))
 138     (separable-discrete-convolve m hr hc :truncate truncate)))
 139
 140 (defparameter *x-derivative-conv-matrix*
 141   (transpose (array->sb8-matrix #2A((1 0 -1)(1 0 -1)(1 0 -1)))))
 142
 143 (defun x-derivative (m &key (matrix-class 'double-float-matrix) (truncate t))
 144   (let ((m1 (copy-to-matrix-type m matrix-class))
 145         (m2 (copy-to-matrix-type *x-derivative-conv-matrix* matrix-class)))
 146     (mat-scale (discrete-convolve m1 m2
 147                                   :truncate truncate :matrix-class matrix-class)
 148                (/ 3))))
 149
 150 (defparameter *y-derivative-conv-matrix*
 151   (array->sb8-matrix #2A((1 0 -1)(1 0 -1)(1 0 -1))))
 152
 153 (defun y-derivative (m &key (matrix-class 'double-float-matrix) (truncate t))
 154   (mat-scale
 155    (discrete-convolve (copy-to-matrix-type m matrix-class)
 156                       (copy-to-matrix-type *y-derivative-conv-matrix* matrix-class)
 157                       :truncate truncate :matrix-class matrix-class)
 158    (/ 3)))
 159
 160 (defun gradmag (m &key (truncate nil))
 161   (let ((xd (x-derivative m :truncate truncate))
 162         (yd (y-derivative m :truncate truncate)))
 163     (mat-square! xd)
 164     (mat-square! yd)
 165     (mat-sqrt! (mat-add xd yd :in-place t))))
 166
 167 (defun graddir (m &key (truncate nil))
 168   (let ((xd (x-derivative m :truncate truncate))
 169         (yd (y-derivative m :truncate truncate)))
 170     (destructuring-bind (ur uc) (dim xd)
 171       (dotimes (i ur)
 172         (dotimes (j uc)
 173           (set-val xd i j
 174                    (if (zerop (val yd i j))
 175                        (val yd i j)
 176                        (atan (/ (val xd i j) (val yd i j))))))))
 177     xd))
 178
 179 (defparameter *laplacian-conv-matrix*
 180   (array->sb8-matrix #2A((0 1 0)(1 -4 1)(0 1 0))))
 181
 182 (defun laplacian (m &key (matrix-class 'double-float-matrix) (truncate t))
 183   (mat-scale
 184    (discrete-convolve (copy-to-matrix-type m matrix-class)
 185                       (copy-to-matrix-type *laplacian-conv-matrix* matrix-class)
 186                       :truncate truncate :matrix-class matrix-class)
 187    (/ 4)))
 188
 189 (defparameter *laplacian-conv-matrix-2*
 190   (array->sb8-matrix #2A((0 0 1 0 0)
 191                          (0 0 0 0 0)
 192                          (1 0 -4 0 1)
 193                          (0 0 0 0 0)
 194                          (0 0 1 0 0))))
 195
 196
 197 (defun laplacian-2 (m &key (matrix-class 'double-float-matrix) (truncate t))
 198   (mat-scale
 199    (discrete-convolve (copy-to-matrix-type m matrix-class)
 200                       (copy-to-matrix-type *laplacian-conv-matrix-2* matrix-class)
 201                       :truncate truncate :matrix-class matrix-class)
 202    (/ 4)
 203    :in-place t))
 204
 205
 206 (defun variance-window (a &key (k 2))
 207   (destructuring-bind (m n) (dim a)
 208     (let ((zm (1- m))
 209           (zn (1- n)))
 210       (map-matrix-copy a #'(lambda (m i j) (variance-range
 211                                             m
 212                                             (max 0 (- i k))
 213                                             (min zm (+ i k))
 214                                             (max 0 (- j k))
 215                                             (min zn (+ j k))))
 216                        :matrix-class 'double-float-matrix))))
 217
 218 (defun sample-variance-window (a &key (k 1) (truncate nil))
 219   ;;; FIXME remove truncate, I think
 220   (declare (ignore truncate))
 221   (destructuring-bind (m n) (dim a)
 222     (let ((zm (1- m))
 223           (zn (1- n)))
 224       (map-matrix-copy a #'(lambda (m i j) (sample-variance-range
 225                                             m
 226                                             (max 0 (- i k))
 227                                             (min zm (+ i k))
 228                                             (max 0 (- j k))
 229                                             (min zn (+ j k))))))))
 230
 231 (defgeneric morphological-op (u v f))
 232 (defmethod morphological-op ((u matrix) (v matrix) f)
 233   ;;; ur, uc, vr, vc are the number of rows and columns in u and v
 234   (destructuring-bind (ur uc) (dim u)
 235     (destructuring-bind (vr vc) (dim v)
 236       ;;; need a new matrix z to hold the values of the convolved matrix
 237       ;;; dim z should be dim u + dim v - 1
 238       (let ((zr (+ ur vr (- 1)))
 239             (zc (+ uc vc (- 1))))
 240         (let ((z (make-instance (class-of u) :rows zr :cols zc)))
 241           (dotimes (i zr)
 242             (let ((ustartr (max 0 (- i vr -1)))
 243                   (uendr (min (- ur 1) i))
 244                   (vstartr (- vr (max (- vr i) 1)))
 245 ;                 (vendr (- vr (min (- zr i) vr)))
 246                   )
 247               (dotimes (j zc)
 248                 (let ((ustartc (max 0 (- j vc -1)))
 249                       (uendc (min (- uc 1) j))
 250                       (vstartc (- vc (max (- vc j) 1)))
 251                       (acc '()))
 252 ;                     (vendc (- vc (min (- zc j) vc))))
 253 ;                 (print (list i j ";" ustartr uendr ";" ustartc uendc
 254 ;                              ";" vstartr vendr ";" vstartc vendc))
 255                   (do ((urow ustartr (1+ urow))
 256                        (vrow vstartr (1- vrow)))
 257                       ((> urow uendr))
 258                     (do ((ucol ustartc (1+ ucol))
 259                          (vcol vstartc (1- vcol)))
 260                         ((> ucol uendc))
 261                       (setf acc (funcall f acc (val u urow ucol) (val v vrow vcol)))))
 262 ;                     (setf acc (max acc (* (val u urow ucol) (val v vrow vcol))))))
 263                   (set-val-fit z i j acc)))))
 264           z)))))
 265
 266 (defun separable-morphological-op (m h f)
 267   (let ((rowstart (floor (/ (1- (rows h)) 2)))
 268         (rowend (floor (/ (rows h) 2)))
 269         (colstart (floor (/ (1- (cols h)) 2)))
 270         (colend (floor (/ (cols h) 2))))
 271
 272     (declare (dynamic-extent rowstart rowend colstart colend)
 273              (index-type rowstart rowend colstart colend))
 274
 275     (let ((h1 (subset-matrix h rowstart rowend 0 (1- (cols h))))
 276           (h2 (subset-matrix h 0 (1- (rows h)) colstart colend)))
 277       (let*
 278           ((m1 (morphological-op m h1 f))
 279            (m2 (morphological-op m1 h2 f)))
 280         m2))))
 281
 282 (defgeneric dilate (u v))
 283 (defmethod dilate ((u matrix) (v matrix))
 284   (separable-morphological-op u v #'(lambda (acc uval vval)
 285                                       (let ((opval (+ uval vval)))
 286                                         (cond
 287                                          ((null acc) opval)
 288                                          (t (max acc opval)))))))
 289
 290 (defgeneric erode (u v))
 291 (defmethod erode ((u matrix) (v matrix))
 292   (separable-morphological-op u v #'(lambda (acc uval vval)
 293                                       (let ((opval (- uval vval)))
 294                                         (cond
 295                                          ((null acc) opval)
 296                                          (t (min acc opval)))))))
 297
 298 #+nil
 299 (defmethod dilate-orig ((u matrix) r)
 300 ;;; ur, uc, vr, vc are the number of rows and columns in u and v
 301   (destructuring-bind (ur uc) (dim u)
 302     (let ((z (make-instance (class-of u) :rows ur :cols uc)))
 303       (dotimes (i ur)
 304         (let ((ustartr (max 0 (- i r)))
 305               (uendr (min (- ur 1) (+ i r))))
 306           (dotimes (j uc)
 307             (let ((ustartc (max 0 (- j r)))
 308                   (uendc (min (- uc 1) (+ j r)))
 309                   (max-val (val u i j)))
 310               (do ((urow ustartr (1+ urow)))
 311                   ((> urow uendr))
 312                 (do ((ucol ustartc (1+ ucol)))
 313                     ((> ucol uendc))
 314                   (setf max-val (max max-val (val u urow ucol)))))
 315               (set-val z i j max-val)))))
 316       z)))
 317
 318
 319 #+nil
 320 (defmethod erode-orig ((u matrix) r)
 321 ;;; ur, uc, vr, vc are the number of rows and columns in u and v
 322   (destructuring-bind (ur uc) (dim u)
 323     (let ((z (make-instance (class-of u) :rows ur :cols uc)))
 324       (dotimes (i ur)
 325         (let ((ustartr (max 0 (- i r)))
 326               (uendr (min (- ur 1) (+ i r))))
 327           (dotimes (j uc)
 328             (let ((ustartc (max 0 (- j r)))
 329                   (uendc (min (- uc 1) (+ j r)))
 330                   (min-val (val u i j)))
 331               (do ((urow ustartr (1+ urow)))
 332                   ((> urow uendr))
 333                 (do ((ucol ustartc (1+ ucol)))
 334                     ((> ucol uendc))
 335                   (setf min-val (min min-val (val u urow ucol)))))
 336               (set-val z i j min-val)))))
 337       z)))
 338
 339 (defgeneric threshold (u tval &key minval maxval))
 340 (defmethod threshold ((u matrix) (tval number) &key (minval 0) (maxval 255))
 341   (map-set-val-copy u #'(lambda (x) (if (> x tval)
 342                                         maxval
 343                                         minval))))
 344
 345 (defgeneric binary-threshold (u tval))
 346 (defmethod binary-threshold ((u matrix) (tval number))
 347   (destructuring-bind (rows cols)
 348       (clem:dim u)
 349     (let ((thresh (make-instance 'bit-matrix :rows (rows u) :cols (cols u))))
 350       (dotimes (i rows)
 351         (dotimes (j cols)
 352           (setf (mref thresh i j)
 353                 (if (> (mref u i j) tval) 1 0))))
 354       thresh)))
 355
 356 (defgeneric complement-matrix (u &key maxval))
 357 (defmethod complement-matrix ((u matrix) &key (maxval (max-val u)))
 358   (destructuring-bind (rows cols)
 359       (clem:dim u)
 360     (let ((comp (mat-copy-proto u)))
 361       (dotimes (i rows)
 362         (dotimes (j cols)
 363           (setf (mref comp i j) (- maxval (mref u i j)))))
 364       comp)))
 365
 366 (defun standard-deviation (&rest matrices)
 367   (flet ((square (q)
 368            (* q q)))
 369     (when matrices
 370       (let ((x (make-instance 'double-float-matrix
 371                               :rows (rows (car matrices))
 372                               :cols (cols (car matrices))))
 373             (sd (make-instance 'double-float-matrix
 374                                :rows (rows (car matrices))
 375                                :cols (cols (car matrices))))
 376             (n (length matrices)))
 377         (dotimes (i (rows sd))
 378           (dotimes (j (cols sd))
 379             (loop for m in matrices
 380                do
 381                (incf (mref x i j) (mref m i j)))))
 382         (mat-scale x (/ n) :in-place t)
 383         (dotimes (i (rows sd))
 384           (dotimes (j (cols sd))
 385             (loop for m in matrices
 386                do
 387                  (incf (mref sd i j) (square (- (mref m i j) (mref x i j)))))))
 388         (mat-scale sd (/ (1- n)) :in-place t)
 389         (mat-sqrt sd)))))
 390
 391
 392 (defun matrix-means (&rest matrices)
 393   (cond ((null matrices) nil)
 394         ((= (length matrices) 1) (car matrices))
 395         (t
 396          (when matrices
 397            (let ((x (make-instance 'double-float-matrix
 398                                    :rows (rows (car matrices))
 399                                    :cols (cols (car matrices))))
 400                  (n (length matrices)))
 401              (dotimes (i (rows x))
 402                (dotimes (j (cols x))
 403                  (loop for m in matrices
 404                     do
 405                     (incf (mref x i j) (mref m i j)))))
 406              (mat-scale x (/ n) :in-place t)
 407              x)))))
 408
 409 (defun matrix-medians (matrices)
 410   (cond ((null matrices) nil)
 411         ((= (length matrices) 1) (car matrices))
 412         (t
 413          (when matrices
 414            (let ((x (make-instance 'double-float-matrix
 415                                    :rows (rows (car matrices))
 416                                    :cols (cols (car matrices)))))
 417              (dotimes (i (rows x))
 418                (dotimes (j (cols x))
 419                  (loop for m in matrices
 420                     with l
 421                     do (push (mref m i j) l)
 422                     finally (setf (mref x i j)
 423                                   (coerce (ch-util::median l)
 424                                           'double-float)))))
 425              x)))))
 426
 427 (defgeneric matrix-l2-distance (b0 b1 &key dest))
 428
 429 ;; we could probably make this faster by passing in a double-float-matrix that we
 430 ;; could use as a temporary destination matrix
 431 (defmethod matrix-l2-distance ((b0 matrix) (b1 matrix) &key dest)
 432   (declare (optimize (speed 3)
 433                      (safety 0)))
 434   (let ((n (the fixnum (* (the fixnum (clem:rows b0)) (the fixnum (clem:cols b0))))))
 435     (declare (type fixnum n))
 436     (/
 437      ;; if we have a destination matrix, use it, otherwise we need to allocate one
 438      (the double-float
 439        (if dest
 440            (clem:sum
 441             (clem:mat-square!
 442              (progn
 443                (clem::matrix-move b0 dest)
 444                (clem:mat-subtr dest b1 :in-place t)
 445                dest)))
 446            (clem:sum
 447             (clem:mat-square
 448              (clem::copy-to-double-float-matrix (clem:m- b0 b1))))))
 449      n)))