Faster Convolution with Separability Property in Multi-dimensional Signals

Topic: Digital image processing. 

Consider a 2D signal (or impulse response of 2D filter) with rect-shaped.

 This 2D rect signal is separable. It can be represented as dot-product of two independent 1D signals.
As well an impulse signal is also separable.



rect = 
[[0 0 0 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 0] 
 [0 0 1 1 1 1 1 0 0]
 [0 0 1 1 1 1 1 0 0] 
 [0 0 1 1 1 1 1 0 0]
 [0 0 1 1 1 1 1 0 0]
 [0 0 1 1 1 1 1 0 0]
 [0 0 0 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 0]]

as dot product of two 1D signals, 
rect =
[[0] 
 [0]
 [1]
 [1]  dot [[0 0 1 1 1 1 1 0 0]]
 [1]
 [1]
 [1]
 [0]
 [0]]

A 2D filter with impulse response h, separable into two 1D filters with impulse responses h1 and h2 respectively.
Consider a 2D image x is filtered with h. Then x is filtered with h1 dot h2.

With separability property, multi dimensional filtering can be done in multi one-dimensional filterings. Programmatically, it's simpler and faster. The larger the image, the bigger the speed advantage.

Following is a demo code in Python with OpenCV 2.
This is a case of 2D rectangular-frame image signal.

In this case, a rectangular-frame can be represented as combination of dot-product and addition of four 1D signals. 
That is, a bigger-rect minus smaller-rect. While, a rect is separable into two 1D signals. 

Note: some books consider rectangular-frame (and diamond-shaped) as non-separable, based on definition cited in the beginning of this note.

import cv2
import numpy as np
import numpy.random as npr
import time

def demo():
    # random input image
    N = 16
    img = (npr.rand(N, N) * 255.0).astype(np.float32)
    # zero outputs
    img2 = np.zeros((N, N), dtype=np.float32)
    img21 = np.zeros((N, N), dtype=np.float32)
    img22 = np.zeros((N, N), dtype=np.float32)
    img23 = np.zeros((N, N), dtype=np.float32)
    img24 = np.zeros((N, N), dtype=np.float32)

    # one 2D kernel
    kernel1 = np.array([[0, 0, 0, 0, 0, 0, 0, 0, 0],
                        [0, 1, 1, 1, 1, 1, 1, 1, 0],
                        [0, 1, 1, 1, 1, 1, 1, 1, 0],
                        [0, 1, 1, 0, 0, 0, 1, 1, 0],
                        [0, 1, 1, 0, 0, 0, 1, 1, 0],
                        [0, 1, 1, 0, 0, 0, 1, 1, 0],
                        [0, 1, 1, 1, 1, 1, 1, 1, 0],
                        [0, 1, 1, 1, 1, 1, 1, 1, 0],
                        [0, 0, 0, 0, 0, 0, 0, 0, 0]])
    
    # four 1D kernels
    kernel11 = np.array([[0], [1], [1], [1], [1], [1], [1], [1], [0]])
    kernel12 = np.array([[0, 1, 1, 1, 1, 1, 1, 1, 0]])
    kernel13 = np.array([[0], [0], [0], [1], [1], [1], [0], [0], [0]])
    kernel14 = np.array([[0, 0, 0, 1, 1, 1, 0, 0, 0]])

    # filtering with one 2D kernel
    t1 = time.time()
    cv2.filter2D(img, cv2.CV_32FC1, kernel1.astype(np.float32), img2)
    print 'time one 2D =', time.time() - t1

    # filtering with four 1D kernels
    t1 = time.time()
    cv2.filter2D(img, cv2.CV_32FC1, kernel11.astype(np.float32), img21)
    cv2.filter2D(img21, cv2.CV_32FC1, kernel12.astype(np.float32), img22)
    cv2.filter2D(img, cv2.CV_32FC1, kernel13.astype(np.float32), img23)
    cv2.filter2D(img23, cv2.CV_32FC1, kernel14.astype(np.float32), img24)
    print 'time four 1D =', time.time() - t1
    
    diff = (img22 - img24) - img2
    diff[np.abs(diff) < 0.001] = 0.0
    print 'diff = (img22 - img24) - img2 :'
    print diff

if __name__ == "__main__":
    demo()

Output:
time one 2D = 0.219000101089
time four 1D = 0.0

The same can be done for diamond-shaped 2D image,

diamond = 
[[0 0 0 0 0 0 0 0 0]
 [0 0 0 0 1 0 0 0 0] 
 [0 0 0 1 1 1 0 0 0]
 [0 0 1 1 1 1 1 0 0] 
 [0 1 1 1 1 1 1 1 0]
 [0 0 1 1 1 1 1 0 0]
 [0 0 0 1 1 1 0 0 0]
 [0 0 0 0 1 0 0 0 0]
 [0 0 0 0 0 0 0 0 0]]
This can be separated into several 1D signals, with combination of dot-product and addition operation on the convolution.

Reference:

• "Signal Processing for Communication", Chap.13, Paolo Prandoni and Martin Vetterli, Ecole Polytechnique Fédérale de Lausanne (EPFL) Press. 
• "Digital Signal Processing", a Coursera course by Ecole Polytechnique Fédérale de Lausanne (EPFL).

/end/