Opencv 傅立叶变换

使用离散二维傅立叶变换（Discrete Fourier Transformation），将灰度化的imori.jpg表示为频谱图。然后用二维离散傅立叶逆变换将图像复原。

二维离散傅立叶变换是傅立叶变换在图像处理上的应用方法。通常傅立叶变换用于分离模拟信号或音频等连续一维信号的频率。但是，数字图像使用[0,255]范围内的离散值表示，并且图像使用H\times W的二维矩阵表示，所以在这里使用二维离散傅立叶变换。

二维离散傅立叶变换使用下式计算，其中I表示输入图像：
G(k,l)=\frac{1}{H\ W}\ \sum\limits_{y=0}^{H-1}\ \sum\limits_{x=0}^{W-1}\ I(x,y)\ e^{-2\ \pi\ j\ (\frac{k\ x}{W}+\frac{l\ y}{H})}
在这里让图像灰度化后，再进行离散二维傅立叶变换。

频谱图为了能表示复数G，所以图上所画长度为G的绝对值。这回的图像表示时，请将频谱图缩放至[0,255]范围。

二维离散傅立叶逆变换从频率分量G按照下式复原图像：
I(x,y)=\frac{1}{H\ W}\ \sum\limits_{l=0}^{H-1}\ \sum\limits_{k=0}^{W-1}\ G(l,k)\ e^{2\ \pi\ j\ (\frac{k\ x}{W}+\frac{l\ y}{H})}

上式中\exp(j)是个复数，实际编程的时候请务必使用下式中的绝对值形态：

如果只是简单地使用for语句的话，计算量达到128^4，十分耗时。如果善用NumPy的化，则可以减少计算量（答案中已经减少到128^2）。

python实现：

import cv2
import numpy as np
import matplotlib.pyplot as plt


# DFT hyper-parameters
K, L = 128, 128
channel = 3


# DFT
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)
                #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



# Read image
img = cv2.imread("imori.jpg").astype(np.float32)

# DFT
G = dft(img)

# write poser spectal to image
ps = (np.abs(G) / np.abs(G).max() * 255).astype(np.uint8)
cv2.imwrite("out_ps.jpg", ps)

# IDFT
out = idft(G)

# Save result
cv2.imshow("result", out)
cv2.waitKey(0)
cv2.imwrite("out.jpg", out)



"""
fimg = np.fft.fft2(gray)

# 第1象限と第3象限, 第2象限と第4象限を入れ替え
fimg =  np.fft.fftshift(fimg)
print(fimg.shape)
# パワースペクトルの計算
mag = 20*np.log(np.abs(fimg))

# 入力画像とスペクトル画像をグラフ描画
plt.subplot(121)
plt.imshow(gray, cmap = 'gray')
plt.subplot(122)
plt.imshow(mag, cmap = 'gray')
plt.show()
"""

c++实现：

#include <opencv2/core.hpp>
#include <opencv2/highgui.hpp>
#include <iostream>
#include <math.h>
#include <complex>

const int height = 128, width = 128;

struct fourier_str {
  std::complex<double> coef[height][width];
};

// RGB to Gray scale
cv::Mat BGR2GRAY(cv::Mat img){
  // prepare output
  cv::Mat out = cv::Mat::zeros(height, width, CV_8UC1);

  // BGR -> Gray
  for (int y = 0; y < height; y++){
    for (int x = 0; x < width; x++){
      out.at<uchar>(y, x) = (int)((float)img.at<cv::Vec3b>(y, x)[0] * 0.0722 + \
                  (float)img.at<cv::Vec3b>(y, x)[1] * 0.7152 + \
                  (float)img.at<cv::Vec3b>(y, x)[2] * 0.2126);
    }
  }
  return out;
}

// Discrete Fourier transformation
fourier_str dft(cv::Mat img, fourier_str fourier_s){
  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 * M_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 * M_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;
}

// Main
int main(int argc, const char* argv[]){

  // read original image
  cv::Mat img = cv::imread("imori.jpg", cv::IMREAD_COLOR);

  // Fourier coefficient
  fourier_str fourier_s;

  // output image
  cv::Mat out = cv::Mat::zeros(height, width, CV_8UC1);

  // BGR -> Gray
  cv::Mat gray = BGR2GRAY(img);

  // DFT
  fourier_s = dft(gray, fourier_s);

  // IDFT
  out = idft(out, fourier_s);

  //cv::imwrite("out.jpg", out);
  cv::imshow("answer", out);
  cv::waitKey(0);
  cv::destroyAllWindows();

  return 0;
}

输入：

灰度化：

输出：

频谱图：