Opencv 离散余弦变换

JPEG 压缩——第一步为Opencv 离散余弦变换。
imori.jpg灰度化之后，先进行离散余弦变换，再进行离散余弦逆变换吧！

离散余弦变换（Discrete Cosine Transformation）是一种使用下面式子计算的频率变换：
0\leq u,\ v\leq T\\
F(u,v)=\frac{2}{T}\ C(u)\ C(v)\ \sum\limits_{y=0}^{T-1}\ \sum\limits_{x=0}^{T-1}\ I(x,y)\ \cos(\frac{(2\ x+1)\ u\ \pi}{2\ T}\ \cos(\frac{(2\ y+1)\ v\ \pi}{2\ T})\\
C(u)=
\begin{cases}
\frac{1}{\sqrt{2}}& (\text{if}\ u=0)\\
1&(\text{else})
\end{cases}
离散余弦逆变换（Inverse Discrete Cosine Transformation）是离散余弦变换的逆变换，使用下式定义。

在这里，K是决定图像复原时分辨率高低的参数。K=T时，DCT的系数全被保留，因此IDCT时分辨率最大。K=1或K=2时，图像复原时的信息量（DCT系数）减少，分辨率降低。如果适当地设定K，可以减小文件大小。
1\leq K\leq T\\
f(x,y)=\frac{2}{T}\ \sum\limits_{u=0}^{K-1}\sum\limits_{v=0}^{K-1}\ C(u)\ C(v)\ F(u,v)\ \cos(\frac{(2\ x+1)\ u\ \pi}{2\ T})\ \cos(\frac{(2\ y+1)\ v\ \pi}{2\ T})\\
C(u)=
\begin{cases}
\frac{1}{\sqrt{2}}& (\text{if}\ u=0)\\
1&(\text{else})
\end{cases}
在这里我们先将图像分割成8\times8的小块，在各个小块中使用离散余弦变换编码，使用离散余弦逆变换解码，这就是 JPEG的编码过程。现在我们也同样地，把图像分割成8\times8的小块，然后进行离散余弦变换和离散余弦逆变换。

python实现：

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

# DCT hyoer-parameter
T = 8
K = 8
channel = 3

# DCT weight
def w(x, y, u, v):
    cu = 1.
    cv = 1.
    if u == 0:
        cu /= np.sqrt(2)
    if v == 0:
        cv /= np.sqrt(2)
    theta = np.pi / (2 * T)
    return (( 2 * cu * cv / T) * np.cos((2*x+1)*u*theta) * np.cos((2*y+1)*v*theta))

# DCT
def dct(img):
    H, W, _ = img.shape

    F = np.zeros((H, W, channel), dtype=np.float32)

    for c in range(channel):
        for yi in range(0, H, T):
            for xi in range(0, W, T):
                for v in range(T):
                    for u in range(T):
                        for y in range(T):
                            for x in range(T):
                                F[v+yi, u+xi, c] += img[y+yi, x+xi, c] * w(x,y,u,v)

    return F


# IDCT
def idct(F):
    H, W, _ = F.shape

    out = np.zeros((H, W, channel), dtype=np.float32)

    for c in range(channel):
        for yi in range(0, H, T):
            for xi in range(0, W, T):
                for y in range(T):
                    for x in range(T):
                        for v in range(K):
                            for u in range(K):
                                out[y+yi, x+xi, c] += F[v+yi, u+xi, c] * w(x,y,u,v)

    out = np.clip(out, 0, 255)
    out = np.round(out).astype(np.uint8)

    return out



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

# DCT
F = dct(img)

# IDCT
out = idct(F)

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

c++实现：

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


const int height = 128, width = 128, channel = 3;

// DCT hyper-parameter
int T = 8;
int K = 8;

// DCT coefficient
struct dct_str {
  double coef[height][width][channel];
};

// Discrete Cosine transformation
dct_str dct(cv::Mat img, dct_str dct_s){
  double I;
  double F;
  double Cu, Cv;

  for (int ys = 0; ys < height; ys += T){
    for (int xs = 0; xs < width; xs += T){
      for (int c = 0; c < channel; c++){
        for (int v = 0; v < T; v ++){
          for (int u = 0; u < T; u ++){
            F = 0;

            if (u == 0){
              Cu = 1. / sqrt(2);
            } else{
              Cu = 1;
            }

            if (v == 0){
              Cv = 1. / sqrt(2);
            }else {
              Cv = 1;
            }

            for (int y = 0; y < T; y++){
              for(int x = 0; x < T; x++){
                I = (double)img.at<cv::Vec3b>(ys + y, xs + x)[c];
                F += 2. / T * Cu * Cv * I * cos((2. * x + 1) * u * M_PI / 2. / T) * cos((2. * y + 1) * v * M_PI / 2. / T);
              }
            }

            dct_s.coef[ys + v][xs + u][c] = F;
          }
        }
      }
    }
  }

  return dct_s;
}

// Inverse Discrete Cosine transformation
cv::Mat idct(cv::Mat out, dct_str dct_s){
  double f;
  double Cu, Cv;

  for(int ys = 0; ys < height; ys += T){
    for(int xs = 0; xs < width; xs += T){
      for(int c = 0; c < channel; c++){
        for(int y = 0; y < T; y++){
          for(int x = 0; x < T; x++){
            f = 0;

            for (int v = 0; v < K; v++){
              for (int u = 0; u < K; u++){
                if (u == 0){
                  Cu = 1. / sqrt(2);
                } else {
                  Cu = 1;
                }

                if (v == 0){
                  Cv = 1. / sqrt(2);
                } else { 
                  Cv = 1;
                }

                f += 2. / T * Cu * Cv * dct_s.coef[ys + v][xs + u][c] * cos((2. * x + 1) * u * M_PI / 2. / T) * cos((2. * y + 1) * v * M_PI / 2. / T);
              }
            }

            f = fmin(fmax(f, 0), 255);
            out.at<cv::Vec3b>(ys + y, xs + x)[c] = (uchar)f;
          }
        }
      }
    }
  }

  return out;
}



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

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

  // DCT coefficient
  dct_str dct_s;

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

  // DCT
  dct_s = dct(img, dct_s);

  // IDCT
  out = idct(out, dct_s);

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

  return 0;
}

输入：

输出：