SSE指令集優化學習：雙線性插值

對SSE的學習總算邁出了第一步，用2天時間對雙線性插值的代碼進行了優化，現將實現的過程梳理以下，算是對這段學習的一個總結。

1. 什么是SSE

說到SSE，首先要弄清楚的一個概念是SIMD（單指令多數據流，Single Instruction Multiple Data），是一種數據並行技術，能夠在一條指令中同時對多個數據執行運算操作，增加處理器的數據吞吐量。SIMD特別的適用於多媒體應用等數據密集型運算。

1.1 歷史

1996年Intel首先推出了支持MMX的Pentium處理器，極大地提高了CPU處理多媒體數據的能力，被廣泛地應用於語音合成、語音識別、音頻視頻編解碼、圖像處理和串流媒體等領域。但是MMX只支持整數運算，浮點數運算仍然要使用傳統的x87協處理器指令。由於MMX與x87的寄存器相互重疊，在MMX代碼中插入x87指令時必須先執行EMMS指令清除MMX狀態，頻繁地切換狀態將嚴重影響性能。這限制了MMX指令在需要大量浮點運算的程序，如三維幾何變換、裁剪和投影中的應用。
另一方面，由於x87古怪的堆棧式緩存器結構，使得硬件上將其流水線化和軟件上合理調度指令都很困難，這成為提高x86架構浮點性能的一個瓶頸。為了解決以上這兩個問題，AMD公司於1998年推出了包含21條指令的3DNow!指令集，並在其K6-2處理器中實現。K6-2是 第一個能執行浮點SIMD指令的x86處理器，也是第一個支持水平浮點寄存器模型的x86處理器。借助3DNow!，K6-2實現了x86處理器上最快的浮點單元，在每個時鍾周期內最多可得到4個單精度浮點數結果，是傳統x87協處理器的4倍。許多游戲廠商為3DNow!優化了程序，微軟的DirectX 7也為3DNow!做了優化，AMD處理器的游戲性能第一次超過Intel，這大大提升了AMD在消費者心目中的地位。K6-2和隨后的K6-III成為市場上的熱門貨。
1999年，隨着Athlon處理器的推出，AMD為3DNow!增加了5條新的指令，用於增強其在DSP方面的性能，它們被稱為“擴展3DNow!”（Extended 3DNow!）。
為了對抗3DNow!，Intel公司於1999年推出了SSE指令集。SSE幾乎能提供3DNow!的所有功能，而且能在一條指令中處理兩倍多的單精度浮點數；同時，SSE完全支持IEEE 754，在處理單精度浮點數時可以完全代替x87。這迅速瓦解了3DNow!的優勢。
1999年后，隨着主流操作系統和軟件都開始支持SSE並為SSE優化，AMD在其2000年發布的代號為“Thunderbird”的Athlon處理器中添加了對SSE的完全支持（“經典”的Athlon或K7只支持SSE中與MMX有關的部分，AMD稱之為“擴展MMX”即Extended MMX）。隨后，AMD致力於AMD64架構的開發；在SIMD指令集方面，AMD跟隨Intel，為自己的處理器添加SSE2和SSE3支持，而不再改進3DNow!。
2010年八月，AMD宣布將在新一代處理器中取消除了兩條數據預取指令之外3DNow!指令的支持，並鼓勵開發者將3DNow!代碼重新用SSE實現。

1.2 MMX和SSE

MMX 是Intel在Pentium MMX中引入的指令集。其缺點是占用浮點數寄存器進行運算（64位MMX寄存器實際上就是浮點數寄存器的別名）以至於MMX指令和浮點數操作不能同時工作。為了減少在MMX和浮點數模式切換之間所消耗的時間，程序員們盡可能減少模式切換的次數，也就是說，這兩種操作在應用上是互斥的。后來Intel在此基礎上發展出SSE指令集；AMD在此基礎上發展出3D Now指令集。
SSE(Streaming SIMD Extensions)是Intel在3D Now!發布一年之后，在PIII中引入的指令集，是MMX的超集。AMD后來在Athlon XP中加入了對這個指令集的支持。這個指令集增加了對8個128位寄存器XMM0-XMM7的支持，每個寄存器可以存儲4個單精度浮點數。使用這些寄存器的程序必須使用FXSAVE和FXRSTR指令來保持和恢復狀態。但是在PIII對SSE的實現中，浮點數寄存器又一次被新的指令集占用了，但是這一次切換運算模式不是必要的了，只是SSE和浮點數指令不能同時進入CPU的處理線而已。
SSE2是Intel在P4的最初版本中引入的，但是AMD后來在Opteron 和Athlon 64中也加入了對它的支持。這個指令集添加了對64位雙精度浮點數的支持，以及對整型數據的支持，也就是說這個指令集中所有的MMX指令都是多余的了，同時也避免了占用浮點數寄存器。這個指令集還增加了對CPU的緩存的控制指令。AMD對它的擴展增加了8個XMM寄存器，但是需要切換到64位模式(AMD64)才可以使用這些寄存器。Intel后來在其EM64T架構中也增加了對AMD64的支持。
SSE3是Intel在P4的Prescott版中引入的指令集，AMD在Athlon 64的第五個版本中也添加了對它的支持。這個指令集擴展的指令包含寄存器的局部位之間的運算，例如高位和低位之間的加減運算；浮點數到整數的轉換，以及對超線程技術的支持。

2 雙線性差值的優化

上面的多半是粘貼的，是前期學習SSE的資料搜集，算是對SSE的由來有一個大致的了解。下面介紹對雙線性插值的優化的學習過程。

2.1 雙線性插值

在圖像變換時，變換后圖像的像素映射到源圖像上的坐標有可能是一個浮點坐標，插值算法就是要計算出浮點坐標像素近似值。那么要如何計算浮點坐標的近似值呢。一個浮點坐標必定會被四個整數坐標所包圍，將這個四個整數坐標的像素值按照一定的比例混合就可以求出浮點坐標的像素值。混合比例為距離浮點坐標的距離，這就是雙線性插值的基本思想。關於雙線性插值的更多信息可以參見WIKI和 本人博文
要優化的雙線性插值的C++實現

	//計算縮放后的圖像大小
	dstWidth = static_cast<int>(width * fx);
	dstHeight = static_cast<int>(height * fy);

	depth /= 8;
	int dstSize = dstWidth * dstHeight * depth;
	dst = new byte[dstSize];
	memset(dst, 255, dstSize);

	byte* dstPixel = nullptr;

	double x = 0.0f; //縮放后圖像在映射到原圖像的橫坐標
	double y = 0.0f; //映射到原圖像的橫坐標

	for (int j = 0; j < dstHeight; j++)
	{
		y = j / fy;

		for (int i = 0; i < dstWidth; i++)
		{
			x = i / fx;

			dstPixel = dst + (j * dstWidth + i) * depth;

			//計算距離當前映射點(x,y)最近的4個點
			int x1, y1, x2, y2;
			x1 = static_cast<int>(x);
			x2 = x1 + 1;
			y1 = static_cast<int>(y);
			y2 = y1 + 1;

			double u = x - x1; //映射點和最左邊橫坐標的差值
			double v = y - y1; //映射點和最下邊縱坐標的差值

			byte cltr1, cltr2, cltr3, cltr4;
			//分別計算各個通道的分量
			for (int k = 0; k < depth; k++)
			{
					cltr1 = src[(y1 * width + x1) * depth + k]; // x1,y1
					cltr2 = src[(y1 * width + x2) * depth + k]; // x2,y1
					cltr3 = src[(y2 * width + x1) * depth + k]; // x1,y2
					cltr4 = src[(y2 * width + x2) * depth + k]; // x2,y2

					double f1, f2;
					f1 = u * cltr1 + (1-u) * cltr2;
					f2 = u * cltr3 + (1-u) * cltr4;
					dstPixel[k] = static_cast<byte>(v * f1 + (1-v) * f2); 
			}
		}
	}

優化的思路，使用128位的xmm寄存器，可以同時對4個32位的數據進行運算，選擇同時對一個像素的所有通道進行處理（3通道和4通道）。優化后的運行速度，因為主要是用來學習SSE指令的，沒有很做很嚴謹的對比，在debug下大致是3倍速（實驗用的圖像是24位3通道），在release下優化和優化的速度沒有明顯的提升，編譯器的優化就是牛哇...。
下面就開始介紹用SSE實現雙線性插值的過程，主要分為以下幾個步驟：

1. 數據的移動
2. 數據的組織：pack,unpck,shuffle
3. 運算

2.2 數據的移動

使用SSE的第一步就是需要將要處理的數據從內存復制到xmm寄存器中。SSE的mov指令可謂無花八門，可以mov不同長度的數據到xmm寄存器，另外有的mov指令需要內存16位對齊，有的則不需要。SSEmov指令可以將數據在xmm寄存器和內存之間進行移動，也可以用於xmm寄存器之間。下面列舉幾個常用mov指令：

1. movd 移動雙字到xmm寄存器
2. movq 移動四字到xmm寄存器
3. movapd / movaps 移動對齊的pack的雙精度浮點數/單精度浮點數
4. movupd / movups 移動未對齊的pack的雙精度浮點數/單精度浮點數
5. movdqa 移動對齊的雙四字
6. movdqu 移動未對齊的雙四字
   SSE的數據移動指令很多，這里只列舉一些，具體的可以參考Intel開發者手冊。
   在雙線性插值中，由於一個浮點坐標要使用其附近的4個像素近似得到，利用SSE指令同時計算這個4個像素的坐標，計算公式如下：

\[(y_i * width + x_i) * depth,i = 1,2 j = 1,2 \]

首先將這需要用到的數據移動到xmm寄存器中

movd xm0,width
movd xmm1,y1
movd xmm2,y2

width和y1,y2都是32位的整型，所以使用movd將其移動到xmm中。

2.3 數據的組織

數據的組織可以說是使用SSE指令的第一個難點，因為128位寄存器可以同時對4個32位數據進行運算，就需要按照該一定的順序將數據放在寄存器中。SSE中提供了兩種指令來調整數據在xmm寄存器中的順序

1. unpck 交叉組合。
2. shuffle 亂序

2.3.1 unpck

unpck是將源操作數和目的操作中字（雙字，四字具體的交叉單位依賴於不同的指令）重新組合。
unpcklps32位單精度浮點數低位交叉 (unpckhps 32位浮點數高位交叉)

unpckhpd 64位浮點數高位交叉 (unpcklpd 64位浮點數低位交叉)

低位的unpck是按照一定的長度（根據指令的不同，一般有字/雙字/四字,字長為16）取兩操作數的低位數據進行重組，重組后的數據的高位是源操作數的低位，低位是目的操作數的低位。重組的結果保存在目的操作數中。
除了上述的unpck指令外，常用的交叉指令還有

1. punpcklbw 交叉組合低位4字中的字節
2. punpcklwd 交叉組合低位雙字中字
3. punpckldq 交叉組合低位四字中的雙字
4. punpcklqdq 交叉組合低位四字
   利用上面提到的交叉指令，可以將byte變成16位的整型，16位的整型變成32為整型，32位整型變成64位整型，后面有用到再詳述。

2.3.2 shuffle

shuffle 指令功能更贊，它能夠將xmm寄存器中的位調整為我們需要的順序。這里以pshufd為例進行介紹。
pshufd xmm1,xmm2,imm8
pshufd有三個操作數，從左往右，第一個操作數是目的操作數保存結果，第二個操作數是源操作數，第三個操作數是一個8位立即數，指定以怎樣的順序將源操作數中數據保存到目的操作數。
imm8中每2位選擇源操作數的一個雙字(共4個)，00選擇第一個，01選擇第二個，10選擇第三個，11選擇第四個。利用imm8位模式選擇好源操作數的雙字后，其每2位的位段決定着這些雙字如何在目的操作數中排列。[0-1]位選擇的雙字放在目的操作數的[0-31]，[2-3]選擇的放在[32-63]，[4-5]選擇的放在[64,95]，[6-7]選擇的放在[96-127]中。
具體如下圖

例如：

_MM_ALIGN16 int a[4] = {1,2,3,4}; //要變成2,4,3,1的順序

_asm
{
	movdqa xmm0,[a];
	pshufd xmm0,xmm0,0x2D; // 0010 1101   
}

0x2d 就是0010 1101
將從源操作數選擇的雙字 從左到右，由高到底依次存放在目的操作數
00 選擇 1，放在[96-127]
10 選擇 3, 放在[64-95]
11 選擇 4, 放在[32-63]
01 選擇 2, 放在[0-31]

2.3.3 同時計算 y1 * width 和 y2 * width

上面說128位的寄存器可以同時進行4個32位數據的計算，這里要將乘法排除在外，因為32位 *32位需要64位空間才能結果不溢出。當然在保證結果不溢出的前提下，也可以使用32位存放相乘的結果。（一定要在保證結果不溢出的前提下，才能如此處理。我在寫這段程序時，開始用16位保存的相乘結果，最后產生了溢出，調試了好久才找到錯誤所在。）
要同時運算y1 * width 和 y2 * width，就需要將一個寄存器中存放width ,width，一個寄存器存放y1,y2

                    punpckldq xmm0, xmm0; // xmm0 => width,width 低64位
                    movd xmm1, y1;
                    movd xmm2, y2;
                    punpckldq xmm1, xmm2; // xmm1 =>y1,y2 低64位  

執行上面指令后，還是不能進行計算。這是因為 32位的乘法產生64位的結果，也就是
src[0-31] * dst[0-31] => src[0-64]
src[64-95] * dst[64-95] => src[64-127]
所以要調整順序 xmm0中數據的順序為 width 0 width 0 ,xmm1中順序為 y1 0 y2 0 這就要用到上面提到的shuffle指令

                    pshufd xmm1, xmm1, 0xd8; // 0x1101 1000
                    pshufd xmm0, xmm0, 0xd8;
                    pmuludq xmm0, xmm1;  

pmuludq 是乘法指令。xmm0得到的是兩個64位的數據，這里是做圖像縮放用32位存放結果也就足夠了。（因為 這里計算的是圖像中某個像素的通道的索引，即使是4k * 4k * 4 = 0x400 0000，32位足夠保存了）。
得到y1 * width 和 y2 * width后，進行交叉后可以4個一起同時和x1，x2相加

					unpcklpd xmm0, xmm0; // xmm0 y1*width y2*width y1*width y2*width
					movd xmm1, x1;
					punpckldq xmm1, xmm1; // xmm1 => x1,x1
					movd xmm2, x2;
					punpckldq xmm2, xmm2; // xmm2 => x2,x2
					unpcklpd xmm1, xmm2; // xmm1 => x1,x1,x2,x2
					paddd xmm0, xmm1; // xmm0 => y1*width + x1,y2*width + x1,y1*width+x2,y2*width+x2

下面需要計算 $ (y * width + x) * depth$ ,128位寄存器只能同時計算2個32位的乘法，所以講xmm0中的數據拆分放在兩個xmm寄存器中

					pshufd xmm1, xmm0, 0x50; // xmm1 => (x1,y1) (x1,y2)
					pshufd xmm2, xmm0, 0xfa; // xmm2 => (x2,y1) (x2,y2)

然后就和上面的類似過程，將得到的64位結果只取低32位，然后在一起放在同一個xmm寄存器中

					movd xmm3, depth;
					punpckldq xmm3, xmm3;
					unpcklpd xmm3, xmm3;

					pmuludq xmm1, xmm3;
					pmuludq xmm2, xmm3;

					pshufd xmm1, xmm1, 0xd8;
					pshufd xmm2, xmm2, 0xd8;
					punpcklqdq xmm1, xmm2;

到這里 xmm1中就存放着要求的浮點坐標周圍的4個像素的位置，依次為(x1,y1) (x1,y2) (x2,y1) (x2,y2)。將這四個值放在數組中，在后面取像素值時會用到

	//像素的地址偏移量
	__declspec(align(16)) int bitOffset[] = { 0,0,0,0 };
	movdqa[bitOffset], xmm1;//(x1,y1) (x1,y2) (x2,y1) (x2,y2)

2.3.4 位擴展

在對像素的通道進行處理時，讀取的數據是8位，而通常在處理時需要運算是32位的或者64位的（本文中和通道數據做運算的是32位浮點數），這就需要將8位的通道數據進行擴展。
在上面也提到可以利用punpcklbw punpcklwd punpckldq對8位字節進行擴展。下面的SSE指令將8位的通道數據轉換為32位的整型

                    // 取出像素依次存放在xmm0,xmm1,xmm2,xmm3
                    lea eax, bitOffset; // 取出數組地址
                    mov edx, [eax];
                    mov ecx, src;
                    movd xmm0, [ecx + edx]; //三通道像素，讀取了32位，只需要處理24位
                    //將xxm0中的低位4個字節擴展為32位整數（只處理低3個32位整數）
                    psllw xmm1, 16;  //將xxm1置為0
                    punpcklbw xmm0, xmm1;
                    punpcklwd xmm0, xmm1; // 擴展為32位整數 

punpcklbw 交叉組合低位四字中的字節，如果源操作數為0，則可以將低位的8位字節擴展為16位整型。
psllw 是以字為單位（16位）的邏輯左移指令，立即數16是左移的位數，由於是以字為單位左移16位，所以該指令可以將xmm1清0。
'punpcklwd' 交叉組合低位四字中的字，如果源操作數為0，則可以將低位的16位整型擴展為32位整型。

2.3.5 數據類型轉換

數據類型的轉換主要在雙精度64位浮點數和單精度32位浮點，64位浮點數和32位整型以及32位浮點數和32位整型之間，具體的轉換指令如下圖

指令前綴為CVTT的表示帶截斷的類型轉換。
在上面，已將三通道的24為像素數據轉換為3個32位的整型數據存放在xmm寄存器中，但是下面要做的運算是浮點型的，還需要將32位整型轉換為32位的浮點數

cvtdq2ps xmm0, xmm0; //32位整型轉換為單精度浮點數 

2.4 運算

SSE指令集能夠同時在多個數據上執行同一個操作，上面做了那么多操作，就是要將參與運算的數據以需要的格式存放到xmm寄存器中，然后對這些數據同時執行相同的操作。SSE指令集對下面三類常用的運算都有提供

1. 算術運算
   為了適應各種樣的運算環境，SSE也提供了多種樣的算術運算。
   pmuludq無符號雙字乘法，源操作數和目的操作數的第一個和第三個32位無符號整型參與運算，結果為2個64位整型，存放在目的操作數。這條指令是將乘法的結果依次給出了，但是只能同時進行2個32位的運算（結果為64位）。而有可能我們只關系乘法結果的某一部分（例如高位，判斷兩個有符號數的符號是否相同）。
   pmullw pack有符號字乘法取低位，16位的乘法運算，將32位結果的低16位保存到目的操作數。
   

pmulhw pack有符號字乘法取高位，16位的乘法運算，將32位結果的高16位保存到目的操作數。
2. 邏輯
3. 位移
4. 比較
具體的指令這里就不再贅述了，可以參考Intel開發者手冊。

3 小結

初次的學習總結就到這里了。使用SSE指令集可以總結為三步：

1. 將數據從內存移動到xmm寄存器，這里注意指令是否要求16或者32位對齊。
2. 利用各種unpck和shuffle指令，對xmm寄存器中的數據進行操作，將其組織成后面運算需要的格式。
3. 對組織好的數據執行運算，注意數據的溢出以及計算的數據類型（浮點、整型、有符號還是無符號），數據的長度（字，雙字，四字，雙四字）以及指令操作的是高位還是低位等。總之運算時要特別小心，不像在高級語言中那么的省心。
   SSE指令集提供了很豐富的操作，但是Intel的手冊查閱卻不是很方便....。另外，直接在C/C++中內嵌匯編代碼，感覺也很拖沓，下一步准備學習下編譯器提供的Intrinsics。