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#源码解析及调试技巧