卷积可以说是图像处理中最基本的操作。线性滤波通过不同的卷积核,可以产生很多不同的效果。假如有一个要处理的二维图像,通过二维的滤波矩阵(卷积核),对于图像的每一个像素点,计算它的领域像素和滤波器矩阵的对应元素的乘积,然后累加,作为该像素位置的值。关于
图像卷积和滤波的一些知识点
可以参考这篇博客。
下面是通过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)
new_image = new_image.clip(0, 255)
new_image = np.rint(new_image).astype('uint8')
return new_image
if __name__ == "__main__":
image = Image.open("图片.jpg", 'r')
output_path = "./outputPic/"
if not os.path.exists(output_path):
os.mkdir(output_path)
a = np.array(image)
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_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))
通过上面的方法实现卷积操作之后,发现可以通过快速傅里叶变换提升卷积运算的效率。python提供了很多标准工具和封装来计算它。NumPy
和 SciPy
都有经过充分测试的封装好的FFT库,分别位于子模块 numpy.fft
和 scipy.fftpack
。有关FFT算法的原理和推导可以参见参考链接的博客。
xn 到 Xk 的转化就是空域到频域的转换,转化为点值表示法。
def DFT_slow(x):
x = np.asarray(x, dtype=float)
N = x.shape[0]
n = np.arange(N)
k = n.reshape((N, 1))
M = np.exp(-2j * np.pi * k * n / N)
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])
使用numpy进行矩阵向量化。
def FFT_vectorized(x):
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 = min(N, 32)
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)))
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])`:输出结果。