数学原理
使用离散二维傅立叶变换（Discrete Fourier Transformation），将灰度化的imori.jpg表示为频谱图。然后用二维离散傅立叶逆变换将图像复原。
二维离散傅立叶变换是傅立叶变换在图像处理上的应用方法。通常傅立叶变换用于分离模拟信号或音频等连续一维信号的频率。但是，数字图像使用[0,255]范围内的离散值表示，并且图像使用H×W的二维矩阵表示，所以在这里使用二维离散傅立叶变换。
二维离散傅立叶变换使用下式计算，其中I表示输入图像：
G(k,l)=\frac{1}{HW}\sum_{y=0}^{H-1}\sum_{x=0}^{W-1}I(x,y)e^{-2{\pi}j(\frac{kx}{W}+\frac{ly}{H})}
在这里让图像灰度化后，再进行离散二维傅立叶变换。
频谱图为了能表示复数G，所以图上所画长度为G的绝对值。这回的图像表示时，请将频谱图缩放至[0,255]范围。
二维离散傅立叶逆变换从频率分量G按照下式复原图像：
I(x,y)=\frac{1}{HW}\sum_{l=0}^{H-1}\sum_{k=0}^{W-1}G(l,k)e^{2{\pi}j(\frac{kx}{W}+\frac{ly}{H})}
上式中exp(j)是个复数，实际编程的时候请务必使用下式中的绝对值形态1：
如果只是简单地使用for语句的话，计算量达到1284，十分耗时。如果善用NumPy的化，则可以减少计算量（答案中已经减少到1282）。
代码实现
C++(非智能指针)
struct fourier_str
{
std::complex<double>** coef;
};
fourier_str dft(cv::Mat img, int height,int width)
{
fourier_str fourier_s;
fourier_s.coef = new std::complex<double>*[height];
for (int i = 0; i < height; ++i)
{
fourier_s.coef[i] = new std::complex<double>[width];
}
double I;
double theta;
std::complex<double> val;
for (int l = 0; l < height; l++)
{
for (int k = 0; k < width; k++)
{
val.real(0);
val.imag(0);
for (int y = 0; y < height; y++)
{
for (int x = 0; x < width; x++)
{
I = (double)img.at<uchar>(y, x);
theta = -2 * CV_PI * ((double)k * (double)x / (double)width + (double)l * (double)y / (double)height);
val += std::complex<double>(cos(theta), sin(theta)) * I;
}
}
val /= sqrt(height * width);
fourier_s.coef[l][k] = val;
}
}
return fourier_s;
}
// Inverse Discrete Fourier transformation
cv::Mat idft(cv::Mat out, fourier_str fourier_s)
{
double theta;
double g;
std::complex<double> G;
std::complex<double> val;
for (int y = 0; y < height; y++)
{
for (int x = 0; x < width; x++)
{
val.real(0);
val.imag(0);
for (int l = 0; l < height; l++)
{
for (int k = 0; k < width; k++)
{
G = fourier_s.coef[l][k];
theta = 2 * CV_PI * ((double)k * (double)x / (double)width + (double)l * (double)y / (double)height);
val += std::complex<double>(cos(theta), sin(theta)) * G;
}
}
g = std::abs(val) / sqrt(height * width);
out.at<uchar>(y, x) = (uchar)g;
}
}
return out;
}
C++(智能指针)
struct fourier_str {
std::unique_ptr<std::complex<double>[]>* coef;
};
// Discrete Fourier transformation
fourier_str dft(cv::Mat img, int height, int width)
{
fourier_str fourier_s;
fourier_s.coef = new std::unique_ptr<std::complex<double>[]>[height];
for (int i = 0; i < height; ++i) {
fourier_s.coef[i] = std::make_unique<std::complex<double>[]>(width);
}
double I;
double theta;
std::complex<double> val;
for (int l = 0; l < height; l++)
{
for (int k = 0; k < width; k++)
{
val.real(0);
val.imag(0);
for (int y = 0; y < height; y++)
{
for (int x = 0; x < width; x++)
{
I = (double)img.at<uchar>(y, x);
theta = -2 * CV_PI * ((double)k * (double)x / (double)width + (double)l * (double)y / (double)height);
val += std::complex<double>(cos(theta), sin(theta)) * I;
}
}
val /= sqrt(height * width);
fourier_s.coef[l][k] = val;
}
}
return fourier_s;
}
//Inverse Discrete Fourier transformation
cv::Mat idft(cv::Mat out, fourier_str fourier_s) {
double theta;
double g;
std::complex<double> G;
std::complex<double> val;
for (int y = 0; y < height; y++) {
for (int x = 0; x < width; x++) {
val.real(0);
val.imag(0);
for (int l = 0; l < height; l++) {
for (int k = 0; k < width; k++) {
G = fourier_s.coef[l][k];
theta = 2 * CV_PI * ((double)k * (double)x / (double)width + (double)l * (double)y / (double)height);
val += std::complex<double>(cos(theta), sin(theta)) * G;
}
}
g = std::abs(val) / sqrt(height * width);
out.at<uchar>(y, x) = (uchar)g;
}
}
return out;
}
Python
def dft(img):
H, W, _ = img.shape
# Prepare DFT coefficient
G = np.zeros((L, K, channel), dtype=np.complex)
# prepare processed index corresponding to original image positions
x = np.tile(np.arange(W), (H, 1))
y = np.arange(H).repeat(W).reshape(H, -1)
# dft
for c in range(channel):
for l in range(L):
for k in range(K):
#G[l, k, c] = np.sum(img[..., c] * np.exp(-2j * np.pi * (x * k / K + y * l / L))) / np.sqrt(K * L)
G[l, k, c] = np.sum(img[..., c] * np.exp(-2j * np.pi * (x * k / K + y * l / L))) / np.sqrt(K * L)
#for n in range(N):
# for m in range(M):
# v += gray[n, m] * np.exp(-2j * np.pi * (m * k / M + n * l / N))
#G[l, k] = v / np.sqrt(M * N)
return G
# IDFT
def idft(G):
# prepare out image
H, W, _ = G.shape
out = np.zeros((H, W, channel), dtype=np.float32)
# prepare processed index corresponding to original image positions
x = np.tile(np.arange(W), (H, 1))
y = np.arange(H).repeat(W).reshape(H, -1)
# idft
for c in range(channel):
for l in range(H):
for k in range(W):
out[l, k, c] = np.abs(np.sum(G[..., c] * np.exp(2j * np.pi * (x * k / W + y * l / H)))) / np.sqrt(W * H)
# clipping
out = np.clip(out, 0, 255)
out = out.astype(np.uint8)
return out
参考链接
https://zhuanlan.zhihu.com/p/99605178
FAQ
1,上述代码除于sqrt(height * width)是为了保证能量统一,如果idft的不使用除于sqrt(height * width),则可以使用normalize将图像重新归一化