DCT變換可謂是JPEG編碼原理里面數學難度最高的一環，我也是因為DCT變換的算法才對JPEG編碼感興趣（真是不自量力）。這一章我就把我對DCT的研究心得體會分享出來，希望各位大神也不吝賜教。
1.離散余弦變換(DCT)介紹
如果想深入了解這一章，就需要從傅里葉變換開始。學過《信號與系統》或者《數學信號處理》的朋友，肯定都對傅里葉變換這一章特別有印象（mengbi），這里有一個視頻對於理解傅里葉變換有很大的幫助。
我們從離散傅里葉變換也就是DFT這里開始，公式走起：
從公式我們可以看到，如果是一個128個序列，我們可以得到65（128/2+1）個實部和65個虛部頻域的幅值。DFT的數學推導非常復雜，但是代碼及其簡單，用C#表示如下：
1 public static Complex[] Dft(Complex[] y, int len) 2 { 3 Complex[] trans = new Complex[len]; 4 for (int k = 0; k < len/2; k++) 5 { 6 for (int n = 0; n < len-1; n++) 7 { 8 trans[k].Real = trans[k].Real + y[n].Real * Math.Cos(2 * PI * k * n / len); 9 trans[k].Imag = trans[k].Imag - y[n].Real * Math.Sin(2 * PI * k * n / len); 10 } 11 } 12 return trans; 13 }
從代碼可以看出復雜度很高大概是O(n2)，所以需要進行優化，於是有了快速傅里葉變換(FFT)，它是根據離散傅氏變換的奇、偶、虛、實等特性，對離散傅立葉變換的算法進行改進獲得的，學習圖像處理肯定少不了這個算法，但是我們這一章不需要對它進行深入講解。我們再來看DCT變換的公式：
如果只是將公式翻譯成代碼，也不復雜，代碼如下：
1 public static double[,] DCT2D(double[,] input) 2 { 3 //寫死為8，實際變換可以不想等 4 int Width = 8; 5 int Height = 8; 6 double[,] coeffs = new double[Width, Height]; 7 double beta = 1d / Width + 1d / Height; 8 //整體偏移-128 9 for (int i = 0; i < Width; i++) 10 { 11 for (int j = 0; j < Height; j++) 12 { 13 input[i, j] = input[i, j] - 128f; 14 } 15 } 16 for (int u = 0; u < Width; u++) 17 { 18 for (int v = 0; v < Height; v++) 19 { 20 double sum = 0d; 21 for (int x = 0; x < Width; x++) 22 { 23 for (int y = 0; y < Height; y++) 24 { 25 double a = input[x, y]; 26 double b = Math.Cos(((2d * x + 1d) * u * Math.PI) / (2 * Width)); 27 double c = Math.Cos(((2d * y + 1d) * v * Math.PI) / (2 * Height)); 28 sum += a * b * c; 29 } 30 } 31 double alphaU = (u == 0) ? 1 / Math.Sqrt(2) : 1; 32 double alphaV = (v == 0) ? 1 / Math.Sqrt(2) : 1; 33 coeffs[u, v] = sum * beta * alphaU * alphaV; 34 } 35 } 36 return coeffs; 37 }
上面這個代碼復雜度很高，於是就有了后面的優化。
2.DCT算法優化
DCT的優化一直以來都有相關課題，我根據ImageSharp里面提到的相關文獻，記錄我所知道的兩種算法。
第一個文獻里面提到的算法非常規則和簡單，所以被廣泛使用，其流程圖如下:
上面總共有四個步驟，分別出現了三種負號分別是黑店，方框，圓圈，他們表示的數學含義如下：
看這個再對照相應源碼看會很清楚：
1 internal static void fDCT1Dllm_32f(Span<float> x, Span<float> y) 2 { 3 float t0, t1, t2, t3, t4, t5, t6, t7; 4 float c0, c1, c2, c3; 5 float[] r = new float[8]; 6 7 //下面常數是這么計算來的 8 //for(i = 0;i < 8;i++){ r[i] = (float)(cos((double)i / 16.0 * M_PI) * M_SQRT2); } 9 r[0] = 1.414214f; 10 r[1] = 1.387040f; 11 r[2] = 1.306563f; 12 r[3] = 1.175876f; 13 r[4] = 1.000000f; 14 r[5] = 0.785695f; 15 r[6] = 0.541196f; 16 r[7] = 0.275899f; 17 18 const float invsqrt2 = 0.707107f; //(float)(1.0f / M_SQRT2); 19 //const float invsqrt2h = 0.353554f; //invsqrt2*0.5f; 20 21 c1 = x[0]; 22 c2 = x[7]; 23 t0 = c1 + c2; 24 t7 = c1 - c2; 25 c1 = x[1]; 26 c2 = x[6]; 27 t1 = c1 + c2; 28 t6 = c1 - c2; 29 c1 = x[2]; 30 c2 = x[5]; 31 t2 = c1 + c2; 32 t5 = c1 - c2; 33 c1 = x[3]; 34 c2 = x[4]; 35 t3 = c1 + c2; 36 t4 = c1 - c2; 37 38 c0 = t0 + t3; 39 c3 = t0 - t3; 40 c1 = t1 + t2; 41 c2 = t1 - t2; 42 43 y[0] = c0 + c1; 44 y[4] = c0 - c1; 45 y[2] = c2 * r[6] + c3 * r[2]; 46 y[6] = c3 * r[6] - c2 * r[2]; 47 48 c3 = t4 * r[3] + t7 * r[5]; 49 c0 = t7 * r[3] - t4 * r[5]; 50 c2 = t5 * r[1] + t6 * r[7]; 51 c1 = t6 * r[1] - t5 * r[7]; 52 53 y[5] = c3 - c1; 54 y[3] = c0 - c2; 55 c0 = (c0 + c2) * invsqrt2; 56 c3 = (c3 + c1) * invsqrt2; 57 y[1] = c0 + c3; 58 y[7] = c0 - c3; 59 }
對比上面流程圖，是不是很簡單！代碼不難但是數學推導相當不簡單，我反正看了很久(mengbi)。這個代碼運算時間相對就少很多了。
另一個文獻提到的算法復雜度更低，這里我提供了一個博客地址，里面講解的還是比較淺顯易懂的，大致意思就是上面的公式如果將N確定為8，則可以經過一系列矩陣變換，將最終的算法簡化到5次乘法，29次加法。流程圖如下：
帶箭頭的連線表示減,不帶箭頭的邊表示加,a1~a5為乘法的常數。
算法代碼如下：
1 for (i = 0; i < N; i++) 2 { 3 tmp0 = output[i, 0] + output[i, 7]; 4 tmp7 = output[i, 0] - output[i, 7]; 5 tmp1 = output[i, 1] + output[i, 6]; 6 tmp6 = output[i, 1] - output[i, 6]; 7 tmp2 = output[i, 2] + output[i, 5]; 8 tmp5 = output[i, 2] - output[i, 5]; 9 tmp3 = output[i, 3] + output[i, 4]; 10 tmp4 = output[i, 3] - output[i, 4]; 11 12 // Even part 13 tmp10 = tmp0 + tmp3; 14 tmp13 = tmp0 - tmp3; 15 tmp11 = tmp1 + tmp2; 16 tmp12 = tmp1 - tmp2; 17 18 output[i, 0] = tmp10 + tmp11; 19 output[i, 4] = tmp10 - tmp11; 20 21 z1 = (tmp12 + tmp13) * (double)0.707106781; 22 output[i, 2] = tmp13 + z1; 23 output[i, 6] = tmp13 - z1; 24 25 // Odd part 26 tmp10 = tmp4 + tmp5; 27 tmp11 = tmp5 + tmp6; 28 tmp12 = tmp6 + tmp7; 29 30 // The rotator is modified from fig 4-8 to avoid extra negations. 31 z5 = (tmp10 - tmp12) * (double)0.382683433; 32 z2 = ((double)0.541196100) * tmp10 + z5; 33 z4 = ((double)1.306562965) * tmp12 + z5; 34 z3 = tmp11 * ((double)0.707106781); 35 36 z11 = tmp7 + z3; 37 z13 = tmp7 - z3; 38 39 output[i, 5] = z13 + z2; 40 output[i, 3] = z13 - z2; 41 output[i, 1] = z11 + z4; 42 output[i, 7] = z11 - z4; 43 }
3.C#中的SMID
上面通過算法的優化，已經讓DCT算法時間大大縮短了，但是我們上面的算法只是一次8個數的計算，在圖像處理中，一般都是8行8列的像素塊，就意味着調用上面算法16次，有沒有更好的法子呢，這個時候我們需要用到計算機的SIMD（Single Instruction Multiple Data）。
如果你使用C++，肯定對SSE不陌生，這個我就不多介紹了，我下面想介紹在C#中使用intrinsic functions 來操作SIMD指令，那就是Vectors。
這個庫在.NET Framework 4.6 之后引入，通過使用Vector，可以實現硬件加速功能。Vector的使用非常簡單，API設計的很容易上手，直接看我們Imagesharp的源碼，它將8*8個像素分成了16個Vector4：
1 public Vector4 V0L; 2 public Vector4 V0R; 3 4 public Vector4 V1L; 5 public Vector4 V1R; 6 7 public Vector4 V2L; 8 public Vector4 V2R; 9 10 public Vector4 V3L; 11 public Vector4 V3R; 12 13 public Vector4 V4L; 14 public Vector4 V4R; 15 16 public Vector4 V5L; 17 public Vector4 V5R; 18 19 public Vector4 V6L; 20 public Vector4 V6R; 21 22 public Vector4 V7L; 23 public Vector4 V7R;
這樣在運算時，將整個8*8像素塊的行列進行轉置，然后分成左右兩組分別進行DCT快速算法：
1 public static void FDCT8x4_RightPart(ref Block8x8F s, ref Block8x8F d) 2 { 3 Vector4 c0 = s.V0R; 4 Vector4 c1 = s.V7R; 5 Vector4 t0 = c0 + c1; 6 Vector4 t7 = c0 - c1; 7 8 c1 = s.V6R; 9 c0 = s.V1R; 10 Vector4 t1 = c0 + c1; 11 Vector4 t6 = c0 - c1; 12 13 c1 = s.V5R; 14 c0 = s.V2R; 15 Vector4 t2 = c0 + c1; 16 Vector4 t5 = c0 - c1; 17 18 c0 = s.V3R; 19 c1 = s.V4R; 20 Vector4 t3 = c0 + c1; 21 Vector4 t4 = c0 - c1; 22 23 c0 = t0 + t3; 24 Vector4 c3 = t0 - t3; 25 c1 = t1 + t2; 26 Vector4 c2 = t1 - t2; 27 28 d.V0R = c0 + c1; 29 d.V4R = c0 - c1; 30 31 float w0 = 0.541196f; 32 float w1 = 1.306563f; 33 34 d.V2R = (w0 * c2) + (w1 * c3); 35 d.V6R = (w0 * c3) - (w1 * c2); 36 37 w0 = 1.175876f; 38 w1 = 0.785695f; 39 c3 = (w0 * t4) + (w1 * t7); 40 c0 = (w0 * t7) - (w1 * t4); 41 42 w0 = 1.387040f; 43 w1 = 0.275899f; 44 c2 = (w0 * t5) + (w1 * t6); 45 c1 = (w0 * t6) - (w1 * t5); 46 47 d.V3R = c0 - c2; 48 d.V5R = c3 - c1; 49 50 c0 = (c0 + c2) * InvSqrt2; 51 c3 = (c3 + c1) * InvSqrt2; 52 53 d.V1R = c0 + c3; 54 d.V7R = c0 - c3; 55 }
然后在轉置回來，再做DCT快速算法（針對列），這樣就完成了一個8*8的DCT轉換。
public static void TransformFDCT( ref Block8x8F src, ref Block8x8F dest, bool offsetSourceByNeg128 = true) { var temp = new Block8x8F(); src.TransposeInto(ref temp); if (offsetSourceByNeg128) { temp.AddToAllInplace(new Vector4(-128)); } FDCT8x4_LeftPart(ref temp, ref dest); FDCT8x4_RightPart(ref temp, ref dest); dest.TransposeInto(ref temp); FDCT8x4_LeftPart(ref temp, ref dest); FDCT8x4_RightPart(ref temp, ref dest); dest.MultiplyInplace(C_0_125); }
4.最后的話
這一章結束了，JPEG里最難的一部分就講解完了，其中很多數學知識對於我而言只是了解的程度。不得不說將數學公式翻譯到代碼的過程雖然艱辛，但很有趣。能力一般水平有限，在表達當中難免有失妥當，還請各位多加理解。
系列目錄：
ImageSharp源碼詳解之JPEG編碼原理（1）JPEG介紹
ImageSharp源碼詳解之JPEG壓縮原理（3）DCT變換
ImageSharp源碼詳解之JPEG壓縮原理（4）量化
ImageSharp源碼詳解之JPEG壓縮原理（6）C#源碼解析及調試技巧