二維DCT變換 | Python實現

本文轉載自查看原文 2020-03-23 22:33 3687 Python

引言

最近專業課在學信息隱藏與數字水印，上到了變換域隱藏技術，提到了其中的DCT變換，遂布置了一個巨煩人的作業，讓手動給兩個\(8\times8\)的矩陣做二維DCT變換，在苦逼的算了一小時后，我決定放棄，轉而決定寫腳本來解決，\((๑•̀ㅂ•́)و✧\)，正好看網上好像只有matlab的腳本，好像沒人用Python來寫這個，遂打算搞一個（你就是純粹為了偷懶不做作業\((*￣rǒ￣)）\)

二維DCT變換原理

還是要普及一下的嘛，畢竟讓我頭疼了一下午的東西，當然也要好好給你們分享一下啦ԅ(¯﹃¯ԅ)
DCT(Discrete Cosine Transform)，又叫離散余弦變換，它的第二種類型，經常用於信號和圖像數據的壓縮。經過DCT變換后的數據能量非常集中，一般只有左上角的數值是非零的，也就是能量都集中在離散余弦變換后的直流和低頻部分
1、一維DCT變換
要弄懂二維離散余弦變換，首先我們需要先了解它在一維下的情況，具體表達式如下:
\(F(0)=\dfrac{1}{\sqrt{N}}\sum_{x=0}^{N-1} f(x)\) \(...........1\)
\(F(u)=\sqrt{\dfrac{2}{N}}\sum_{x=0}^{N-1} f(x) \cos{\dfrac{(2x+1)u\pi}{2N}}\) \(...........2\)
式中\(F(u)\)是第\(u\)個余弦變換值，\(u\)是廣義頻率變量，\(u=1,2,….,N-1;f(x)\)是時域\(N\)點序列。\(x= 1,2,….,N-1;\)
2、二維DCT變換
二維離散余弦變換可由下列表達式表示：
\(F(0,0)=\dfrac{1}{N}\sum_{x=0}^{N-1}\sum_{y=0}^{N-1} f(x,y)\) \(...........3\)
\(F(0,v)=\dfrac{\sqrt{2}}{N}\sum_{x=0}^{N-1}\sum_{y=0}^{N-1} f(x,y) \cdot \cos{\dfrac{(2y+1)v\pi}{2N}}\) \(...........4\)
\(F(u,0)=\dfrac{\sqrt{2}}{N}\sum_{x=0}^{N-1}\sum_{y=0}^{N-1} f(x,y) \cdot \cos{\dfrac{(2x+1)u\pi}{2N}}\) \(...........5\)
\(F(u,v)=\dfrac{2}{N}\sum_{x=0}^{N-1}\sum_{y=0}^{N-1} f(x,y) \cdot \cos{\dfrac{(2x+1)u\pi}{2N}}\cdot \cos{\dfrac{(2y+1)v\pi}{2N}}\) \(...........6\)
\(6\)是二維離散余弦變換的正變換公式，其中\(f(x,y)\)是空間域一個\(N*N\)的二維向量元素，即一個\(N*N\)的矩陣，\(x,y = 0,1,2，…，N-1;F(U,V)\)是經計算后得到的變換域矩陣，\(u,v = 0,1,2，….，N-1\).求和可分性是二維離散余弦變換的一個重要特征，因此我們可以用下式表示\(6\)：
\(F(u,v)=\dfrac{2}{N}\sum_{x=0}^{N-1} \cdot \cos{\dfrac{(2x+1)u\pi}{2N}}{\sum_{y=0}^{N-1} f(x,y) }\cdot \cos{\dfrac{(2y+1)v\pi}{2N}} ..........7\)
由一維和二維的離散余弦變換公式性質可以推導得到二維離散余弦變換也可以寫成矩陣相乘形式:
\(F=A[f(x,y)A^{T}]..............8\)
\(A\)為一維離散余弦變換的變換系數矩陣，\(A^{T}\)是\(A\)的轉置矩陣
對圖像進行二維離散余弦變換\(（2D-DCT）\)的步驟
1.獲得圖像的二維數據矩陣\(f(x,y)\)；
2.求離散余弦變換的系數矩陣\(A\)；
3.求系數矩陣對應的轉置矩陣\(A^{T}\)；
4.根據公式\(F=A\times[f(x,y)]\times A^{T}\)計算離散余弦變換；
注：公式的大致推導,由於\(A\)都是正交陣，所以\(A^{-1}=A^{T}\),故\(Y=A\times[f(x,y)]\times A^{T} ，f(x,y)=A^{T}\times Y \times A=A^{T}\times A\times[f(x,y)]\times A^{T}\times A\)

參考鏈接：`https://www.cnblogs.com/latencytime/p/10228938.html`

參考鏈接：`https://blog.csdn.net/allen_sdz/article/details/83279210`

Python編程實現

大家注意上述的第\(8\)個式子，將變換的兩個\(\sum_{i=0}^n\)轉變成了變換矩陣和轉置矩陣以及代轉換矩陣之間乘積的問題（注意這里的乘是指的叉乘）
遂可以發現\(DCT\)后的矩陣應等於\(A\times X \times A^{T}\)
此外，對於\(x,y\)同時為0的情況，其參數要單獨考慮
大致的思路就是，由於直接計算變換比較繁瑣，所以我們就先對於一個單位陣進行操作運算，將其變成一個DCT的變換矩陣，而后在與代算矩陣和轉置矩陣叉乘即可。
至此我們就可以着手寫腳本了，這里主要是用了兩個庫numpy和math

# -*- coding:utf-8 -*-
# Author:Konmu
# DCT二維變換

from numpy import array as matrix, arange,zeros,transpose,matmul,ones
from math import sqrt,cos,pi

'''
作業代轉化矩陣1
a=matrix([[0,255,0,255,0,255,0,255],
          [255,0,255,0,255,0,255,0],
          [0,255,0,255,0,255,0,255],
          [255,0,255,0,255,0,255,0],
          [0,255,0,255,0,255,0,255],
          [255,0,255,0,255,0,255,0],
          [0,255,0,255,0,255,0,255],
          [255,0,255,0,255,0,255,0]])
'''
a=ones((8,8))#生成單位陣
for i in range(8):
    a[i]=matrix([128]*8)
# 生成全是128的矩陣（作業代轉化矩陣2）
'''
測試數據
a=matrix([[61,19,50,20],
          [82,26,61,45],
          [89,90,82,43],
          [93,59,53,97]])
'''
A=zeros((8,8))#生成0矩陣
shape=a.shape[1]#獲取維數
for i in range(8):
    for j in range(8):
        if(i == 0):
            x=sqrt(1/shape)
        else:
            x=sqrt(2/shape)
        A[i][j]=x*cos(pi*(j+0.5)*i/shape)#與維數相關

A_T=A.transpose()#矩陣轉置
Y1=matmul(A,a)#矩陣叉乘
Y=matmul(Y1,A_T)
print(Y)
'''
想要近似值可以嘗試這樣輸出
for i in range(shape):
    for j in range(shape):
        print('{:^8.4f}'.format(Y[i][j]),end='\n')
print()
'''