Opencv 离散余弦变换+量化|极客教程

JPEG 压缩——第二步：离散余弦变换+量化

量化离散余弦变换系数并使用离散余弦逆变换恢复。再比较变换前后图片的大小。

量化离散余弦变换系数是用于编码 JPEG 图像的技术。

量化即在对值在预定义的区间内舍入，其中floor、ceil、round等是类似的计算。

在 JPEG 图像中，根据下面所示的量化矩阵量化离散余弦变换系数。该量化矩阵取自 JPEG 软件开发联合会组织颁布的标准量化表。在量化中，将8x 8的系数除以（量化矩阵） Q 并四舍五入。之后然后再乘以 Q 。对于离散余弦逆变换，应使用所有系数。

Q = np.array(((16, 11, 10, 16, 24, 40, 51, 61),
              (12, 12, 14, 19, 26, 58, 60, 55),
              (14, 13, 16, 24, 40, 57, 69, 56),
              (14, 17, 22, 29, 51, 87, 80, 62),
              (18, 22, 37, 56, 68, 109, 103, 77),
              (24, 35, 55, 64, 81, 104, 113, 92),
              (49, 64, 78, 87, 103, 121, 120, 101),
              (72, 92, 95, 98, 112, 100, 103, 99)), dtype=np.float32)

由于量化降低了图像的大小，因此可以看出数据量已经减少。

python实现：

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

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

# DCT weight
def DCT_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] * DCT_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] * DCT_w(x,y,u,v)

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

    return out

# Quantization
def quantization(F):
    H, W, _ = F.shape

    Q = np.array(((16, 11, 10, 16, 24, 40, 51, 61),
                (12, 12, 14, 19, 26, 58, 60, 55),
                (14, 13, 16, 24, 40, 57, 69, 56),
                (14, 17, 22, 29, 51, 87, 80, 62),
                (18, 22, 37, 56, 68, 109, 103, 77),
                (24, 35, 55, 64, 81, 104, 113, 92),
                (49, 64, 78, 87, 103, 121, 120, 101),
                (72, 92, 95, 98, 112, 100, 103, 99)), dtype=np.float32)

    for ys in range(0, H, T):
        for xs in range(0, W, T):
            for c in range(channel):
                F[ys: ys + T, xs: xs + T, c] =  np.round(F[ys: ys + T, xs: xs + T, c] / Q) * Q

    return F



# MSE
def MSE(img1, img2):
    H, W, _ = img1.shape
    mse = np.sum((img1 - img2) ** 2) / (H * W * channel)
    return mse

# PSNR
def PSNR(mse, vmax=255):
    return 10 * np.log10(vmax * vmax / mse)

# bitrate
def BITRATE():
    return 1. * T * K * K / T / T


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

# DCT
F = dct(img)

# quantization
F = quantization(F)

# IDCT
out = idct(F)

# MSE
mse = MSE(img, out)

# PSNR
psnr = PSNR(mse)

# bitrate
bitrate = BITRATE()

print("MSE:", mse)
print("PSNR:", psnr)
print("bitrate:", bitrate)

# 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 = 4;

// 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;
}

// Quantization
dct_str quantization(dct_str dct_s){
  double Q[T][T] = {{16, 11, 10, 16, 24, 40, 51, 61},
                    {12, 12, 14, 19, 26, 58, 60, 55},
                    {12, 12, 14, 19, 26, 58, 60, 55},
                    {14, 17, 22, 29, 51, 87, 80, 62},
                    {18, 22, 37, 56, 68, 109, 103, 77},
                    {24, 35, 55, 64, 81, 104, 113, 92},
                    {49, 64, 78, 87, 103, 121, 120, 101},
                    {72, 92, 95, 98, 112, 100, 103, 99}
                  };

  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++){
            dct_s.coef[ys + y][xs + x][c] = round(dct_s.coef[ys + y][xs + x][c] / Q[y][x]) * Q[y][x];
          }
        }
      }
    }
  }

  return dct_s;
}


// Compute MSE
double MSE(cv::Mat img1, cv::Mat img2){
  double mse = 0;

  for(int y = 0; y < height; y++){
    for(int x = 0; x < width; x++){
      for(int c = 0; c < channel; c++){
        mse += pow(((double)img1.at<cv::Vec3b>(y, x)[c] - (double)img2.at<cv::Vec3b>(y, x)[c]), 2);
      }
    }
  }

  mse /= (height * width);
  return mse;
}

// Compute PSNR
double PSNR(double mse, double v_max){
  return 10 * log10(v_max * v_max / mse);
}

// Compute bitrate
double BITRATE(){
  return T * K * K / T * T;
}

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

  double mse;
  double psnr;
  double bitrate;

  // 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);

  // Quantization
  dct_s = quantization(dct_s);

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

  // MSE, PSNR
  mse = MSE(img, out);
  psnr = PSNR(mse, 255);
  bitrate = BITRATE();

  std::cout << "MSE: " << mse << std::endl;
  std::cout << "PSNR: " << psnr << std::endl;
  std::cout << "bitrate: " << bitrate << std::endl;

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

  return 0;
}

输入：

Opencv 离散余弦变换+量化