Convolution theory and numpy implementation

convolution.ipynb

convolution.py

from typing import Tuple, Literal, Callable

import numpy as np
from numpy.lib.stride_tricks import as_strided

# 2D Convolution
def get_dim_ext(dim_k, dim_crop):
    l, k = dim_k
    cr, cc = dim_crop

    if cr * cc *l *k  == 0:
        raise ValueError(f"None of dimension can be zero, l:{l}, k:{k}, cc:{cc}, cr:{cr}") 
    if cr <0 or cc<0 or l<0 or k<0:
        raise ValueError(f"None of dimension can be negative, l:{l}, k:{k}, cc:{cc}, cr:{cr}") 
    if cr > l or cc >k:
        raise ValueError(f"Cropping dimension must be smaller than kernel dimension.")
    return l-cr, k-cc

def get_dim_convole(dim_i, dim_k, dim_c):
    n, m = dim_i
    l, k = dim_k

    if l> n or k>m:
        n ,l = l, n
        m, k = k, m

    cr, cc = dim_c
    
    return n+l - 2*cr+1, m+k - 2*cc +1

def _expand_matrix(data, dim_ext, edge_param):
    er, ec = dim_ext

    if edge_param[0] == "custom" and len(edge_param) >=2:
        C1, C2, R1, R2  = edge_param[1](data, (er, ec), edge_param[2:])
    else: C1, C2, R1, R2 = _edge_matrix(data, (er, ec), edge_param)

    # Column expand
    if isinstance(C1, np.ndarray):
        data = np.concatenate([C1, data, C2], axis=1)
    # Row expand
    if isinstance(R1, np.ndarray):
        data = np.concatenate([R1, data, R2], axis=0)

    return data

def _edge_matrix(data, dim_ext, edge_param):
    n, m = data.shape
    er, ec = dim_ext
    # Column Matrix: (n, ec)
    # Row Matrix: (er, m+2ec) 
    dim_column = (n, ec)
    dim_row = (er, int(m+2*ec))

    e_method_name, e_method_param = edge_param

    if e_method_name == "extend":
        
        a_col_1 = np.expand_dims(data[0:, 0], axis=0)
        a_col_2 = np.expand_dims(data[0:,-1], axis=0)
        C1 = np.repeat(a_col_1, ec, axis=0).transpose()
        C2 = np.repeat(a_col_2, ec, axis=0).transpose()

        a_row_1 = np.expand_dims(np.concatenate([C1[0],data[0], C2[0]]), axis=0)
        a_row_2 = np.expand_dims(np.concatenate([C1[-1], data[-1], C2[-1]]), axis=0)
        R1 =np.repeat(a_row_1, er, axis=0)
        R2 =np.repeat(a_row_2, er, axis=0) 

    elif e_method_name ==  "wrap":
        if er >0 and ec >0:
            Cor1 = data[-er:,   -ec:  ] 
            Cor2 = data[-er:,      :ec] 
            Cor3 = data[   :er, -ec:  ]
            Cor4 = data[   :er,    :ec]

            Row1 = data[-er:,    :]
            Row2 = data[   :er,  :]

            R1 = np.concatenate([Cor1, Row1, Cor2], axis=1)
            R2 = np.concatenate([Cor3, Row2, Cor4], axis=1)

            C1 = data[:, -ec:  ] 
            C2 = data[:,    :ec]
        elif er>0: # ec ==0
            R1 = data[-er:,    :]
            R2 = data[   :er,  :]
            C1 = C2 = None
        else:
            R1 = R2 = None
            C1 = data[:, -ec:  ] 
            C2 = data[:,    :ec]

    elif e_method_name == "reflect":
        if er * ec >0:
            C1 = np.flip(data[:,    :ec], axis=1)
            C2 = np.flip(data[:, -ec:  ], axis=1)

            Cor1 = np.flip(C1[   :er, :], axis=0) #
            Cor2 = np.flip(C2[   :er, :], axis=0) #

            Cor3 = np.flip(C1[-er:  , :], axis=0) #
            Cor4 = np.flip(C2[-er:  , :], axis=0) # 

            Row1 = np.flip(data[:er, :], axis=0)
            Row2 = np.flip(data[-er:, :], axis=0)

            R1 = np.concatenate([Cor1 , Row1, Cor2], axis =1) 
            R2 = np.concatenate([Cor3 , Row2, Cor4], axis =1)
        elif er > 0:
            R1 = np.flip(data[:er, :], axis=0)
            R2 = np.flip(data[-er:, :], axis=0)
            C1 = C2 = None
        else:
            R1 = R2 = None
            C1 = np.flip(data[:,    :ec], axis=1)
            C2 = np.flip(data[:, -ec:  ], axis=1)
    elif e_method_name == "mirror":
        if er*ec >0:
            C1 = np.flip(data[:,    1:ec+1], axis=1)
            C2 = np.flip(data[:, -(ec+1): -1  ], axis=1)

            Cor1 = np.flip(C1[   1:er+1, :], axis=0)
            Cor2 = np.flip(C2[   1:er+1, :], axis=0)
            Cor3 = np.flip(C1[-er-1: -1  , :], axis=0)
            Cor4 = np.flip(C2[-er-1: -1  , :], axis=0)

            Row1 = np.flip(data[1:er+1, :], axis=0)
            Row2 = np.flip(data[-(er+1):-1, :], axis=0)

            R1 = np.concatenate([Cor1 , Row1, Cor2], axis =1) 
            R2 = np.concatenate([Cor3 , Row2, Cor4], axis =1)
        elif er >0:
            R1 = np.flip(data[1:er+1, :], axis=0)
            R2 = np.flip(data[-(er+1):-1, :], axis=0)
            C1 = C2 = None
        else:
            R1 = R2 = None
            C1 = np.flip(data[:,    1:ec+1], axis=1)
            C2 = np.flip(data[:, -(ec+1): -1  ], axis=1)
    elif e_method_name == "constant":
        c = e_method_param
        C1 = C2 = c*np.ones(shape=dim_column) if ec>0 else None
        R1 = R2 = c*np.ones(shape=dim_row) if er>0 else None
    
    else:
        C1 = C2 = R1 = R2 = None

    return C1, C2, R1, R2

def convolve2d(  ## 1dim with 2dim convolution
    data:np.ndarray, 
    filter:np.ndarray,
    crop:Tuple[int, int]=(1, 1),
    edge_mode:Literal["extend", "wrap", "mirror", "reflect", "constant"]="constant",
    edge_params = [0],
    image_mode = False
    ):
    
    #no, mo = data.shape
    #l, k = filter.shape
    #if l*k > no*mo and not image_mode:
    #    data, filter = filter, data

    l, k = filter.shape
    filter = np.flip(filter, axis=0) # axis=1 flip occured at np.convolve routine.
    
    if image_mode:
        add_row = add_col = False
        crop = (int(np.floor((l+1)/2)), int(np.floor((k+1)/2)))
        if l%2 == 0:
            add_row = True
        if k%2 == 0:
            add_col = True
    

    er, ec = get_dim_ext((l, k), crop)

    data_ext = _expand_matrix(data, (er, ec), [edge_mode]+edge_params)

    n,m = data_ext.shape
    result = np.zeros(shape = (n-l+1, m-k+1))
    for i, h_row in enumerate(filter):
        data_i = data_ext[i:n-l+i+1]
        result += np.stack([np.convolve(h_row, r, mode="vaild") for r in data_i])
    
    if image_mode:
        if add_row:
            result = result[1:]
        if add_col:
            result = result[:, 1:]
    return result

#------------------------------------------------------
# Utils
def restore_data(data, dim_k, dim_crop):
    l, k = dim_k
    cr, cc = dim_crop
    er = l - cr
    ec = k -cc
    return data[er: -er, ec: -ec]


# To matrix system

def _vec2sub_toeplitz(vec, m): #m: column dimension of matrix
    k = len(vec)
    n = m-k +1
    rows = []
    for i in range(0, n):
        row_zeros =  np.zeros(shape=(m-k,))
        rows.append(np.insert(row_zeros, i, vec))
    return np.vstack(rows)


def convolve2toeplitz(
    data:np.ndarray, 
    filter:np.ndarray, 
    crop:Tuple[int, int]=(1, 1),
    edge_mode:Literal["extend", "wrap", "mirror", "reflect","constant"]="constant",
    edge_params = [0], 
    allow_ext=False) -> Tuple[np.ndarray, np.ndarray, Callable]:
    l, k = filter.shape
    er, ec = get_dim_ext((l, k), crop)
    data_ext = _expand_matrix(data, (er, ec), [edge_mode]+edge_params)

    n,m = data_ext.shape

    filter = np.flip(np.flip(filter, axis=1), axis=0)

    H_list = []
    for i in range(0, l):
        h_i = _vec2sub_toeplitz(filter[i], m)
        if not allow_ext:
            h_i = h_i[:, ec:-ec]
        H_list.append(h_i)
    
    topelitz_dim = H_list[0].shape

    rows =[]
    for i in range(0, n-l+1):
        i_f = i
        i_b = n- l -i_f

        row = i_f*[np.zeros(shape=topelitz_dim)] + H_list + i_b*[np.zeros(shape=topelitz_dim)]
        rows.append(row)
    
    mat = np.block(rows)

    if not allow_ext:
        c = topelitz_dim[1]
        erc = er*c
        mat = mat[:, erc:-erc]
        vec = data.flatten()
    else:
        vec = data_ext.flatten()
    return mat, vec

#-------------------------------------------------------------------------------------------------------
# Special Case
#---------------------------------------------------------------

def get_matrix_system(filter, dim_d):
    n,m = dim_d
    l, k= filter.shape

    if k != 2*m-1 or l != 2*n-1:
        raise ValueError("Invaild dimension: l, k must be 2n-1, 2m-1")
    
    rows = []
    for i in range(0, n):
        row_i = n-1-i
        row_f = 2*n-1-i

        for j in range(0, m):
            column_i = m-1 -j
            column_f = 2*m-1 -j

            t = filter[row_i : row_f, column_i:column_f]

            rows.append(t.flatten())

    return np.vstack(rows)

The above Toeplitz transform is one example of basic matrix representation.
The ouputs are same with next paper includingreshape routine.

Michal Gnacij, Krystian Lapa, Using Toeplitz matrices to obtain 2D convolution, posted: October 27th, 2022, doi:https://doi.org/10.21203/rs.3.rs-2195496/v1

See codes of their transform in their repository, 2D convolution

There is a image processing algorithm that use a toeplitz matrix

Image deconvolution via superfast inversion of a class of two-level Toeplitz matrices
10.1109/ICIP.2012.6467549

For higher dimension tensor convolution see

Nussbaumer transformation

10.1109/TASSP.1980.1163372