卷积可以说是图像处理中最基本的操作。线性滤波通过不同的卷积核,可以产生很多不同的效果。假如有一个要处理的二维图像,通过二维的滤波矩阵(卷积核),对于图像的每一个像素点,计算它的领域像素和滤波器矩阵的对应元素的乘积,然后累加,作为该像素位置的值。关于 图像卷积和滤波的一些知识点 可以参考这篇博客。

下面是通过python模拟实现的图像卷积操作,模拟了 sobel 算子, prewitt 算子和拉普拉斯算子。python的np包中有 convolve 函数可以直接调用, scipy 中也有 scipy.signal.convolve 函数可以直接调用。

import matplotlib.pyplot as plt
import pylab
import cv2
import numpy as np
from PIL import Image
import os
def conv(image, kernel):
    height, width = image.shape        # 获取图像的维度
    h, w = kernel.shape                 # 卷积核的维度
    # 经过卷积操作后得到的新的图像的尺寸
    new_h = height - h + 1
    new_w = width - w + 1
    # 对新的图像矩阵进行初始化
    new_image = np.zeros((new_h, new_w), dtype=np.float)     
    # 进行卷积操作,矩阵对应元素值相乘
    for i in range(new_w):
        for j in range(new_h):
            new_image[i, j] = np.sum(image[i:i+h, j:j+w] * kernel)    # 矩阵元素相乘累加
    # 去掉矩阵乘法后的小于0的和大于255的原值,重置为0和255
    # 用clip函数处理矩阵的元素,使元素值处于(0,255)之间
    new_image = new_image.clip(0, 255)   
    # 将新图像各元素的值四舍五入,然后转成8位无符号整型
    new_image = np.rint(new_image).astype('uint8')     
    return new_image
if __name__ == "__main__":
    # 读取图像信息,并转换为numpy下的数组
    image = Image.open("图片.jpg", 'r')
    output_path = "./outputPic/"
    if not os.path.exists(output_path):
        os.mkdir(output_path)
    a = np.array(image)
    # sobel 算子
    sobel_x = np.array(([-1, 0, 1],
                        [-2, 0, 2],
                        [-1, 0, 1]))
    sobel_y = np.array(([-1, -2, -1],
                        [0, 0, 0],
                        [1, 2, 1]))
    sobel = np.array(([-1, -1, 0],
                      [-1, 0, 1],
                      [0, 1, 1]))
    # prewitt各个方向上的算子
    prewitt_x = np.array(([-1, 0, 1],
                          [-1, 0, 1],
                          [-1, 0, 1]))
    prewitt_y = np.array(([-1, -1, -1],
                          [0, 0, 0],
                          [1, 1, 1]))
    prewitt = np.array(([-2, -1, 0],
                        [-1, 0, 1],
                        [0, 1, 2]))
    # 拉普拉斯算子
    laplacian = np.array(([0, -1, 0],
                          [-1, 4, -1],
                          [0, -1, 0]))
    laplacian_2 = np.array(([-1, -1, -1],
                            [-1, 8, -1],
                            [-1, -1, -1]))
    kernel_list = ("sobel_x", "sobel_y", "sobel", "prewitt_x", "prewitt_y", "prewitt", "laplacian", "laplacian_2")
    print("Gridient detection\n")
    for w in kernel_list:
        print("starting %s....." % w)
        print("kernel:\n")
        print("R\n")
        R = conv(a[:, :, 0




    
], eval(w))
        print("G\n")
        G = conv(a[:, :, 1], eval(w))
        print("B\n")
        B = conv(a[:, :, 2], eval(w))
        I = np.stack((R, G, B), axis=2)     # 合并三个通道的结果
        Image.fromarray(I).save("%s//bigger-%s.jpg" % (output_path, w))

快速傅里叶变换(FFT)

通过上面的方法实现卷积操作之后,发现可以通过快速傅里叶变换提升卷积运算的效率。python提供了很多标准工具和封装来计算它。NumPySciPy 都有经过充分测试的封装好的FFT库,分别位于子模块 numpy.fftscipy.fftpack 。有关FFT算法的原理和推导可以参见参考链接的博客。

离散傅里叶变换

xn 到 Xk 的转化就是空域到频域的转换,转化为点值表示法。

def DFT_slow(x):
    # Compute the discrete Fourier Transform of the 1D array x
    x = np.asarray(x, dtype=float)              # 转化为ndarray
    N = x.shape[0]                              # 维度
    n = np.arange(N)                            # 0~N组成一个一维向量
    k = n.reshape((N, 1))                       # 转换为一个N维向量
    M = np.exp(-2j * np.pi * k * n / N)         # 离散傅里叶公式 -2j复数表示
    return np.dot(M, x)

快速傅里叶变换

离散傅里叶变换具有对称性,利用这种对称性,可以推出递归的快速傅里叶算法。

def FFT(x):
    # A recursive implementation of the 1D Cooley-Tukey FFT
    x = np.asarray(x, dtype=float)                 # 浅拷贝
    N = x.shape[0]
    if N % 2 > 0:
        raise ValueError("size of x must be a power of 2")
    elif N <= 32:
        return DFT_slow(x)
    else:
        X_even = FFT(x[::2])          # 从0开始,2为间隔
        X_odd = FFT(x[1::2])          # 从1开始,2为间隔
        factor = np.exp(-2j * np.pi * np.arange(N) / N)
        使用/会出现下面的错误,改为// 向下取整
        TypeError: slice indices must be integers or None or have an __index__ method
        return np.concatenate([X_even + factor[:N // 2] * X_odd,
                               X_even + factor[N // 2:] * X_odd])

向量化的FFT

使用numpy进行矩阵向量化。

def FFT_vectorized(x):
    # A vectorized, non-recurisive version of the Cooley_Tukey FFT
    x = np.asarray(x, dtype=float)
    N = x.shape[0]
    if np.log2(N) % 1 > 0:
        raise ValueError("size of x must be a power of 2")
    # N_min here is equivalent to the stopping condition above,
    # and should be a power of 2
    N_min = min(N, 32)
    # Perform an O[N^2] DFT on all length-N_min sub-problems at once
    n = np.arange(N_min)
    k = n[:, None]
    M = np.exp(-2j * np.pi * n * k / N_min)
    X = np.dot(M, x.reshape((N_min, -1)))
    # build-up each level of the recursive calculation all at once
    while X.shape[0] < N:
        X_even = X[:, :X.shape[1] // 2]
        X_odd = X[:, X.shape[1] // 2:]
        factor = np.exp(-1j * np.pi * np.arange(X.shape[0]) / X.shape[0])[:, None]
        X = np.vstack([X_even + factor * X_odd,
                       X_even - factor * X_odd])
    return X.ravel()     # 降维

用快速傅里叶变换实现卷积

根据卷积定理,时域上卷积运算对应于频域上的傅里叶变换的乘积。

def fft_convolve(a, b):
    n = len(a) + len(b) -1
    N = 2 ** (int(np.log2(n))+1)
    A = np.fft.fft(a, N)
    B = np.fft.fft(b, N)
    return np.fft.ifft(A*B)[:n]

c++实现的递归版

void FFT(Complex* a,int len){
    if(len==1) return;
    Complex* a0=new Complex[len/2];
    Complex* a1=new Complex[len/2];
    for(int i=0;i<len;i+=2){
        a0[i/2]=a[i];
        a1[i/2]=a[i+1];
    FFT(a0,len/2);FFT(a1,len/2);
    Complex wn(cos(2*Pi/len),sin(2*Pi/len));
    Complex w(1,0);
    for(int i=0;i<(len/2);i++){
        a[i]=a0[i]+w*a1[i];
        a[i+len/2]=a0[i]-w*a1[i];
        w=w*wn;
    return;

c++实现的非递归版

const double PI = acos(-1.0);
// 复数
struct complex
    double r,i;
    complex(double _r = 0.0, double _i = 0.0)
        r = _r;
        i = _i;
    complex operator +(const complex &b)
        return complex(r + b.r, i + b.i);
    complex operator -(const complex &b)
        return complex(r - b.r, i - b.i);
    complex operator *(const complex &b)
        return complex(r*b.r - i*b.i, r*b.i + i*b.r);
// 雷德算法 -- 倒位序
void Rader(complex F[], int len)
    int j = len >> 1;
    for(int i=1; i<len-1; i++)
        if(i < j)
            swap(F[i], F[j]);
        int k = len >> 1;
        while(j >= k)
            j -= k;
            k >>= 1;
        if(j < k)
            j += k;
void FFT(complex F[], int len, int on)
    Rader(F, len);
    for(int h=2; h<=len; h<<=1)    //计算长度为h的DFT
        complex wn(cos(-on*2*PI/h), sin(-on*2*PI/h));   //单位复根e^(2*PI/m),用欧拉公式展开
        for(int j=0; j<len; j+=h)
            complex w(1,0);       //旋转因子
            for(int k=j; k<j+h/2; k++)
                complex u = F[k];
                complex t = w * F[k + h/2];
                F[k] = u + t;        //蝴蝶合并操作
                F[k + h/2] = u - t;
                w = w * wn;          //更新旋转因子
    //求逆傅里叶变换
    if(on == -1)
        for(int i=0; i<len; i++)
            F[i].r /= len;
//求卷积
void Conv(complex a[], complex b[], int len)
    FFT(a, len, 1);
    FFT(b, len, 1);
    for(int i=0; i<len; i++)
        a[i] = a[i] * b[i];         //卷积定理
    FFT(a, len, -1);

参考链接:
http://blog.jobbole.com/58246/ (快速傅里叶变换)
https://blog.csdn.net/acdreamers/article/details/39005227 (快速傅里叶变换)
http://blog.miskcoo.com/2015/04/polynomial-multiplication-and-fast-fourier-transform#i-15 (快速傅里叶变换)
https://blog.csdn.net/WADuan2/article/details/79529900 (快速傅里叶变换)
https://blog.csdn.net/oliverkingli/article/details/79243731 (图像卷积)
https://blog.csdn.net/zlh_hhhh/article/details/75604333 (快速傅里叶变换)
https://www.cnblogs.com/youmuchen/p/6724780.html (图像卷积)
http://blog.sina.com.cn/s/blog_154d272c90102x98p.html (雷德算法)
https://www.cnblogs.com/Sakits/p/8405098.html (快速傅里叶变换)

卷积运算卷积可以说是图像处理中最基本的操作。线性滤波通过不同的卷积核,可以产生很多不同的效果。假如有一个要处理的二维图像,通过二维的滤波矩阵(卷积核),对于图像的每一个像素点,计算它的领域像素和滤波器矩阵的对应元素的乘积,然后累加,作为该像素位置的值。关于图像卷积和滤波的一些知识点可以参考这篇博客。下面是通过python模拟实现的图像卷积操作,模拟了sobel算子,prewitt算子和拉普... 几种加速CNN训练的方法,同时不会对模型最终的准确率有较大的影响。 全连接层是网络占用大量内存的主要原因,但是运行速度却很快,而卷积虽然减少了参数数量,但却需要消耗更多的计算资源,因此,应该从卷积操作本身入手,寻找优化的方法。 以下几种方法可以在不严重影响准确率的前提下加速卷积计算(有些... 取L为它们二者长度的最大值,不够的补零 这个在书上90页 对于线性卷积而言。它用于LTI系统,他也分为时域和频域的,线性卷积有物理意义,但是没有高效算法(这仅仅是站在之前的角度来看的),现在这一节我们将共同探讨线性卷积快速算法问题 首先请看书上149页中间的推导,这... 1. 前言 傅立叶变换具有卷积的性质,因此可以使用傅立叶变换实现数据卷积处理。图像的傅立叶变换可以参考这里:数字图像处理与Python实现-图像变换-傅立叶变换 2. 傅立叶卷积性质数学表示 如果f(x)f(x)f(x)和g(x)g(x)g(x)为两个一维时域函数;f(x,y)f(x,y)f(x,y)和g(x,y)g(x,y)g(x,y)为两个二维空域函 for (int k = ; k < ROWS; k++) { for (int l = ; l < COLS; l++) { sum += input[i + k][j + l] * kernel[k][l]; output[i][j] = sum; // 输出结果 for (int i = ; i < 3; i++) { for (int j = ; j < 3; j++) { printf("%d ", output[i][j]); printf("\n"); return ; 1. `#include <stdio.h>`:引入标准输入输出库。 2. `#define ROWS 3` 和 `#define COLS 3`:定义卷积核的行数和列数。 3. `int kernel[ROWS][COLS]`:定义卷积核。 4. `int input[5][5]`:定义输入矩阵。 5. `int output[3][3]`:定义输出矩阵。 6. `for (int i = ; i < 3; i++)` 和 `for (int j = ; j < 3; j++)`:遍历输出矩阵的每一个元素。 7. `for (int k = ; k < ROWS; k++)` 和 `for (int l = ; l < COLS; l++)`:遍历卷积核的每一个元素。 8. `sum += input[i + k][j + l] * kernel[k][l]`:计算卷积结果。 9. `output[i][j] = sum`:将卷积结果存入输出矩阵。 10. `printf("%d ", output[i][j])`:输出结果。