SSE圖像算法優化系列二十四: 基於形態學的圖像後期抗鋸齒算法--MLAA優化研究。
偶爾看到這樣的一個算法，覺得還是蠻有意思的，花了將近10天多的時間研究了下相關代碼。
以下為百度的結果：MLAA全稱Morphological Antialiasing，意為形態抗鋸齒是AMD推出的完全基於CPU處理的抗鋸齒解決方案。對於遊戲廠商使用的MSAA抗鋸齒技術不同，Intel最新推出的MLAA將跨越邊緣像素的前景和背景色進行混合，用第2種顏色來填充該像素，從而更有效地改進圖像邊緣的變現效果，這就是MLAA技術。
其實就是這個是由Intel的工程師先於2009年提出的技術，但是由AMD將其發發揚光大。
整個算法的渲染工作全部是交給CPU來完成，在這裏GPU的作用只是將最終渲染出來的畫面傳給顯示器。所以這項技術最大的優勢是可以讓GPU不再承擔抗鋸齒的工作，大大降低GPU在運行3D遊戲時的壓力。相對於以前的抗鋸齒技術，MLAA采用Post-filtering（後濾波）機制，好處就在於可以按照顏色是否連續來驅動抗鋸齒，而以前只能在初始邊緣來抗鋸齒。
也就是說這項技術可以在後期來修補那些由鋸齒的圖，因此我們可以想到其另外一些用處，後續會對這方面進行一個簡單的擴展。
如上面兩圖，左側圖中樹葉的邊緣有明顯的鋸齒狀圖像，而右側為經過MLAA算法處理後的圖，邊緣光滑了許多，而且其他部位未受任何的畫質影響。
關於這方面的論文和資料主要有Morphological Antialiasing.pdf ，Intel官方還對改算法進行了一些額外的說明，詳見：https://software.intel.com/sites/default/files/m/d/4/1/d/8/MLAA.pdf，並且附帶了相關源代碼，代碼可從https://software.intel.com/zh-cn/articles/morphological-antialiasing-mlaa-sample處下載。我們重點來研讀和改進下這部分代碼。
下載代碼後，找到\SAA-samples\SAA文件夾下的 MLAAPostProcess.h和MLAAPostProcess.cpp文件，這就是我們最關心的CPU實現的代碼。
根據論文的描述，MLAA算法共包含3個步驟：
1、尋找在特定的圖像像素之間的不連續性，在有些圖像中梯度幅值較大的並不是邊緣點。
2、確定預定模式，確定渲染的圖像。
3、在預定模式中進行領裏邊緣色彩混合處理。確定模式中的相應模板。
下面的文章如果你沒有看過代碼或者沒有看過論文，你根本不知道我在說什麽。
第一步的計算連續性，在實現上實際是計算一個像素點和其右側及下方一個像素的顏色差異絕對值的大小，如果某個點和其下方像素的 AbsDiff大於某個指定的閾值，則設置這個點為水平方向邊緣的標記（EdgeFlagH)，如果和其右側像素的顏色差異大於閾值，則設置這個點為為垂直方向邊緣的標記（EdgeFlagV）。一個點可以只是水平邊緣或垂直邊緣，也可以只是其中一個，或者兩者都不是。
在獲得連續性的基礎上，第二步是沿著圖像的某一個方向，比如寬度方向分析邊緣的形狀，這裏可以有Z型，U形和L型，針對不同的形狀我們有不同的處理方式。
最後，就是在獲得形狀後，按照一定的規則對這些比較硬性的拐角處(L和Z型都是直角的彎），進行融合和柔化。
由於MLAA算法的初衷是處理GPU初處理的圖，因此其主要針對的圖像必然是32位圖像，而且圖像的寬度和高度一般來說都是4的倍數，因此，在Intel給出的代碼中我們可以看到都是處理BGRA格式的圖像的。
具體到代碼實時上，因為是和顯示打交道的，因此 ，算法的實時性必須有可靠的保證，不然這個算法粗在的意義就會大打折扣，網絡上說人有提供了改算法的GPU實現，但是同時又提到那個GPU代碼還不如沒有。我想也是，Intel作為改算法的提出者，所分享的代碼確實還是相當又學習意義的，我這裏來稍微分析下。
首先是不連續性的相關代碼，對於32位圖形，Intel只計算其BGR三通道和周邊像素的差異，當三通道其中一個通道的差異絕對值大於16時，我們就認為這個點是某個方位的邊緣處。另外，在Intel的代碼中，其將這些標誌信息直接隱藏到了BGRA的A通道中，這個可能於DDS格式有關，DDS格式有一些並沒有用到所有的Alpha的8位信息。這樣做的一個好處是不需要額外的內存用來保存邊緣標誌。
為了速度起見，同時考慮在CPU端運行，這部分可充分利用SSE進行優化。我們先貼出Intel的相關代碼進行分析。
//-------------------------------------------------------------------------------------------------------------------
// This task analyzes the color buffer for discontinuities between pixels to set the edge flags in the work buffer.
//-------------------------------------------------------------------------------------------------------------------
void MLAAFindDiscontinuitiesTask(VOID* pvInfo, INT, UINT uTaskId, UINT)
{
MLAAData* pMLAATaskData = static_cast<MLAAData*>(pvInfo);
UINT* ColorBuffer = pMLAATaskData->pColorBuffer;
UINT BufferWidth = pMLAATaskData->uWidth;
UINT BufferHeight = pMLAATaskData->uHeight;
UINT BufferPitch = pMLAATaskData->uRowPitch; // in pixels, not bytes.
UINT TaskFirstRow = uTaskId * RowsOrColsPerTask;
UINT TaskLastRow = TaskFirstRow + RowsOrColsPerTask;
for(UINT iRow = TaskFirstRow; iRow < TaskLastRow; iRow += 4)
{
UINT BlockFirstRowOffset = iRow * BufferPitch;
for(UINT iCol = 0; iCol < BufferWidth; iCol += 4)
{
__m128i *PixelData = reinterpret_cast<__m128i*>(ColorBuffer + BlockFirstRowOffset + iCol);
// BufferPitch is in pixels i.e. 32-bit units, and PixelData points to a 128-bit data type...
UINT BufferPitch128 = BufferPitch/4;
// Load pixel block from color buffer. vPixels0 contains the 4 pixels indexed by iRow and iCol,
// vPixels1, the 4 pixels just below the pixels in vPixels0, and so on and so forth.
__m128i vPixels0 = _mm_load_si128(PixelData);
__m128i vPixels1 = _mm_load_si128(PixelData + BufferPitch128);
__m128i vPixels2 = _mm_load_si128(PixelData + 2*BufferPitch128);
__m128i vPixels3 = _mm_load_si128(PixelData + 3*BufferPitch128);
__m128i vPixels4 = (iRow == BufferHeight - 4) ?
vPixels3 // For the last block (vertically), add a virtual row by duplicating the last real block row.
: _mm_load_si128(PixelData + BufferPitch);
// zero alpha, we are using it to store discontinuity flags.
__m128i vZeroAlpha = _mm_set1_epi32(0x00FFFFFF);
vPixels0 = _mm_and_si128(vPixels0, vZeroAlpha);
vPixels1 = _mm_and_si128(vPixels1, vZeroAlpha);
vPixels2 = _mm_and_si128(vPixels2, vZeroAlpha);
vPixels3 = _mm_and_si128(vPixels3, vZeroAlpha);
// Check for horizontal pixel discontinuities, one row of 4 pixels checked per call.
// (we compare a row of 4 pixels with its bottom neighbor)
ComparePixelsSSE(vPixels0, vPixels1, EdgeFlagH);
ComparePixelsSSE(vPixels1, vPixels2, EdgeFlagH);
ComparePixelsSSE(vPixels2, vPixels3, EdgeFlagH);
ComparePixelsSSE(vPixels3, vPixels4, EdgeFlagH);
// Transpose pixel block so we can use ComparePixelsSSE to check for vertical discontinuities.
_MM_TRANSPOSE4_PS(
reinterpret_cast<__m128&>(vPixels0),
reinterpret_cast<__m128&>(vPixels1),
reinterpret_cast<__m128&>(vPixels2),
reinterpret_cast<__m128&>(vPixels3) );
vPixels4 = (iCol == BufferWidth - 4) ?
vPixels3 // For the last block (horizontally), add a virtual column by duplicating the last real column.
: _mm_setr_epi32(
* reinterpret_cast<UINT*>(PixelData + 1),
* reinterpret_cast<UINT*>(PixelData + BufferPitch128 + 1),
* reinterpret_cast<UINT*>(PixelData + 2*BufferPitch128 + 1),
* reinterpret_cast<UINT*>(PixelData + 3*BufferPitch128 + 1) );
// Now vPixels0..4 contains the block of pixel data in column order, e.g.
// vPixels0 now stores the leftmost column of 4 pixels, vPixels1 the 2nd leftmost one, etc.
// Now check for vertical pixel discontinuities, one column of 4 pixels checked per call.
// (we compare a column of 4 pixels with its neighbor on the right)
ComparePixelsSSE(vPixels0, vPixels1, EdgeFlagV);
ComparePixelsSSE(vPixels1, vPixels2, EdgeFlagV);
ComparePixelsSSE(vPixels2, vPixels3, EdgeFlagV);
ComparePixelsSSE(vPixels3, vPixels4, EdgeFlagV);
// Transpose back and store in color buffer the pixel data with added discontinuities flags.
_MM_TRANSPOSE4_PS(
reinterpret_cast<__m128&>(vPixels0),
reinterpret_cast<__m128&>(vPixels1),
reinterpret_cast<__m128&>(vPixels2),
reinterpret_cast<__m128&>(vPixels3) );
_mm_store_si128(PixelData, vPixels0);
_mm_store_si128(PixelData + BufferPitch128, vPixels1);
_mm_store_si128(PixelData + 2*BufferPitch128, vPixels2);
_mm_store_si128(PixelData + 3*BufferPitch128, vPixels3);
}
}
}
第一感覺就是代碼的註釋很豐富，不愧是大家制作，在解釋這段代碼之前，我們先來看看ComparePixelsSSE函數。
//--------------------------------------------------------------------------------------------------------------
// Given a row of 4 pixels, check for color discontinuities between a pixel and its neighbor.
// SSE implementation, so we process the 4 consecutive pixels at a time.
// If a discontinuity is detected, the flag passed as 3rd arg. is set in the alpha component of the pixel.
// vPixels0 is the vector of 4 pixels to be flagged if discontinuities are detected.
// vPixels1 is the vector of 4 neighbor pixels.
//--------------------------------------------------------------------------------------------------------------
inline void ComparePixelsSSE(__m128i& vPixels0, __m128i vPixels1, const INT EdgeFlag)
{
// Each byte of vDiff is the absolute difference between a pixel‘s color channel value, and the one from its neighbor.
__m128i vDiff = _mm_sub_epi8(_mm_max_epu8(vPixels0, vPixels1), _mm_min_epu8(vPixels0, vPixels1));
// We are only interested if the difference is greater than 16, and not interested in alpha differences.
const INT Threshold = 0x00F0F0F0;
__m128i vThresholdMask = _mm_set1_epi32(Threshold);
vDiff = _mm_and_si128(vDiff, vThresholdMask);
// Each of the four lanes of the selector is 0 if no discontinuity detected RGB, 0xFFFFFFFF else.
__m128i vSelector = _mm_cmpeq_epi32(vDiff, _mm_setzero_si128());
__m128i vEdgeFlag = _mm_set1_epi32(EdgeFlag);
vEdgeFlag = _mm_andnot_si128(vSelector, vEdgeFlag);
// vEdgeFlag now contains EdgeFlag for the pixels where a discontinuity was detected, 0 otherwise.
vPixels0 = _mm_or_si128(vPixels0, vEdgeFlag);
}
ComparePixelsSSE的作用就是比較8個BGRA像素的差異，如果像素的BGR差異值有一個大於閾值，則設置A的某個位為Flag。
第一行vDiff的計算就是計算差異值的絕對值，他一次性可以計算16個字節值的差異，因為SSE沒有直接提供這樣的函數，因此，這裏使用max和min函數結合實現，也是非常的巧妙，其實還有另外一個方式可以實現這個功能，如下所示：
// 返回8位字節數數據的差的絕對值
inline __m128i _mm_absdiff_epu8(__m128i a, __m128i b)
{
return _mm_or_si128(_mm_subs_epu8(a, b), _mm_subs_epu8(b, a));
}
利用了飽和計算，同樣也是十分的巧妙。後續測試表面上述方法速度似乎還能有一定的優勢。
vDiff = _mm_and_si128(vDiff, vThresholdMask); 這一句結合__m128i vSelector = _mm_cmpeq_epi32(vDiff, _mm_setzero_si128());是這個函數的靈魂所在，我們註意到vThresholdMask對應的字節範圍內數據的高4位都是1，低4位都為0，因此如果vDiff中某個值大於16（高4位有值不為1），則進行and運算後返回值必然不為0，4個字節中只要有一個不為0，組成的32位數也必然不為0，這樣一個BGRA像素只要有一個通道有AbsDiff大於0的值，就可以通過_mm_cmpeq_epi32的比較（和0比較）返回值得以體現。
後面使用_mm_andnot_si128是因為我們_mm_cmpeq_epi32的函數在等於0時返回0xFFFFFF，而我們實際上在此時想讓他為0x00000000，因為系統不提供_mm_cmpneq_epi32這樣的函數，所以要先取反下(not)，然後在和Flag進行And運算，特別需要註意的是，不想大多數的SSE函數，_mm_andnot_si128這個SSE函數對參數的順序是銘感的，如果放錯了位置，則無法得到正確的結果。 另外，這裏const INT EdgeFlagH = (1 << 31)以及const INT EdgeFlagV = (1 << 30)，Intel也是考慮的很好，第一，和他們進行and操作不會影響BGR的值，第二，後面還有一個技巧也和這個值有關。
最後的_mm_or_si128主要是考慮水平和垂直方向的邊緣的綜合識別。
雖然我們知道對32位數，SSE有_mm_cmpge_epi32這個比較函數，但是在這裏確實無法直接使用他。
我們再來回頭看MLAAFindDiscontinuitiesTask這個函數，他的核心是一次性讀取4行4列共16個BGRA像素，占用4個XMM寄存器得大小，然後在函數內部一次性的計算出水平邊緣和垂直邊緣的標記，這樣左一個核心的好處是能減少讀取內存的次數。那麽這裏最靈活的運用就是_MM_TRANSPOSE4_PS這個宏的使用。我們看他的名字，就知道他是針對浮點進行轉置使用，而我們的32位圖像加載後是__m128i類型的，但是其實在計算機內部，不管你表面是浮點的還是整形的，都是把數據保存在XMM寄存上，因此，我們用reinterpret_cast或者_mm_castsi128_ps這樣的一些語法糖來強制轉換，讓編譯器能編譯通過，_MM_TRANSPOSE4_PS這個照樣可以處理整形的。
比如，_mm_castsi128_ps這個函數intel給出的解釋是Cast vector of type __m128i to type __m128. This intrinsic is only used for compilation and does not generate any instructions, thus it has zero latency. 也就是說他並沒有產生任何額外的指令，只是個語法糖。
當然，我們不用_MM_TRANSPOSE4_PS也可以用整形的unpack系列函數實現32位整形的轉置，實測他們效率基本無區別。
話題扯得遠了點，MLAAFindDiscontinuitiesTask中，首先一次性加載四行各4各像素後，一次進行水平方向的比較，水平比較完後，這樣對這些數據進行轉置，有可以繼續進行垂直方向的比較，比較完成後，再轉置回來，就一次性完成了水平和垂直方向的所有比較。因此，速度就能得到大幅的提高。
這一步的優化完成，第二步中我們要尋找不同的模式，其中一個關鍵的步驟就是尋找具有同一邊緣標記的連續線段。這部分再函數裏實現，如下所示：
//-----------------------------------------------------------------------------------------------------------------------------------------
// Given a range of pixels offsets [OffsetCurrentPixel, OffsetEndRow], walk this range to find a discontinuity line.
// We call a "separation line" or "discontinuity line" a sequence of consecutive pixels that all have the same edge flag set.
// The 3 return values are: the length of the separation line, and the offsets in the buffer of the first and last pixel of the line.
//-----------------------------------------------------------------------------------------------------------------------------------------
inline UINT FindSeparationLine(int& OutOffsetLineStart, int& OutOffsetLineEnd, UINT* WorkBuffer, UINT OffsetCurrentPixel, UINT OffsetEndRow, const INT EdgeFlag)
{
if (OffsetCurrentPixel >= OffsetEndRow)
{ // We are done scanning this row/column; no separation line left to find.
return 0;
}
// Find first extremity of the line... 找尋線的端頭
OutOffsetLineStart = -1;
int Shift = (EdgeFlag == EdgeFlagV) ? 1 : 0;
for (;;)
{
__m128i PixelFlags = _mm_loadu_si128((__m128i*)(WorkBuffer + OffsetCurrentPixel));
PixelFlags = _mm_slli_epi32(PixelFlags, Shift);
// Creates a 4-bit mask from the edge flag of each of the 4 pixels
int HFlags = _mm_movemask_ps(_mm_castsi128_ps(PixelFlags));
if (HFlags)
{
unsigned long Index;
_BitScanForward(&Index, HFlags);
OffsetCurrentPixel += Index;
OutOffsetLineStart = OffsetCurrentPixel;
break;
}
OffsetCurrentPixel += 4; // None of the pixels had its flag set so we can jump ahead 4 pixels...
if (OffsetCurrentPixel >= OffsetEndRow)
{ // Done scanning this row/column, without finding a separation line.
OutOffsetLineStart = OutOffsetLineEnd = OffsetEndRow - 1;
return 0;
}
}
FindEndOffset:
// Now look for the second extremity of the line.
// We could use a SSE-based optimization as the one above, but it remains to be seen if it would help significantly the performance.
// (discontinuity lines tend to be short)
UINT LineLength = 1;
++OffsetCurrentPixel;
while ((OffsetCurrentPixel <= OffsetEndRow) && ((WorkBuffer[OffsetCurrentPixel] & EdgeFlag) != 0))
{
++LineLength;
++OffsetCurrentPixel;
}
OutOffsetLineEnd = OffsetCurrentPixel - 1;
return LineLength;
}
好像這段代碼我稍微修改了下，源代碼有對OffsetCurrentPixel不是4的整數倍做判斷，其實那個判斷是基於默認使用(__m128i *)這種方式加載變量到SSE寄存器是采用的_mm_load_si128（即16字節對齊有關），如果顯式的使用_mm_loadu_si128則無需那樣寫。
我們看下這裏的技巧主要在於_mm_movemask_ps的使用，在第一步不連續性的標記過程中，我們把邊緣比較分別設置在低31位（水平邊緣）和低30（垂直邊緣）位中，因此，如果是尋找垂直邊緣是，我們把整體向左移動移位，這個時候在他們強制轉換為浮點數，此時_mm_movemask_ps就會根據XMM寄存中每個32位的浮點數的sign返回一個值，如果返回值的都為0，則表示4各數都為正數，說明這4個位置都沒有我們需要尋找的標記，如果結果不為0，這個時候我們應該從低位到高位尋找到第一個不為0的位，他對應的索引就是第一個標記所在的位置，這種位尋找的函數在Intel裏正好有提供，即本例的_BitScanForward函數。
這個方式在很多的優化或計算中都可以借鑒，確實是個不錯的方法。
那麽當找到第一個位置的標記後，我們就需要找到最後一個標記，一般情況下連續的統一標記不會太長，因此後續的尋找不需要借助SSE處理。
還有一個使用了SSE優化的地方就在於最後一步的融合地方，如下所示函數：
inline void MixColors(UINT* ColorBuffer, UINT OffsetDst, float Weight1, UINT Offset1, float Weight2, UINT Offset2)
{
// Load pixels and convert the bytes to dwords
__m128i Col1i = _mm_cvtsi32_si128(ColorBuffer[Offset1]);
__m128i Col2i = _mm_cvtsi32_si128(ColorBuffer[Offset2]);
__m128i Zero = _mm_setzero_si128();
Col1i = _mm_unpacklo_epi16(_mm_unpacklo_epi8(Col1i, Zero), Zero);
Col2i = _mm_unpacklo_epi16(_mm_unpacklo_epi8(Col2i, Zero), Zero);
// Convert to floats so the pixels can be multplied with the weights
__m128 Col1f = _mm_cvtepi32_ps(Col1i);
__m128 Col2f = _mm_cvtepi32_ps(Col2i);
Col1f = _mm_mul_ps(Col1f, _mm_set1_ps(Weight1));
Col2f = _mm_mul_ps(Col2f, _mm_set1_ps(Weight2));
Col1f = _mm_add_ps(Col1f, Col2f);
// Go back to byte from float
__m128i Coli = _mm_cvttps_epi32(Col1f);
Coli = _mm_packs_epi32(Coli, Coli);
Coli = _mm_packus_epi16(Coli, Coli);
// Store the weighted pixel
ColorBuffer[OffsetDst] = (ColorBuffer[OffsetDst] & 0xFF000000) | (_mm_cvtsi128_si32(Coli) & 0x00FFFFFF);
}
這個其實也很簡單，就是4個字節數據乘以4個浮點數，然後累加，最後結果保存為字節數，並且要保留部分字節不變的一個過程，有興趣的朋友可以自己研讀下。
那麽在Intel的代碼中，還利用了TBB進行優化，將計算分成了多任務處理過程，另外，考慮垂直方向的cache命中率問題，將垂直計算過程使用轉置後在調用水平方向的算法進行處理。
那麽後面我在談下這個代碼裏幾個細節的東西。
（1）代碼裏有ComputeUpperBounds和ComputeLowerBounds函數，他們的主要作用是輔助確認模式，代碼本身不難，但是如果你只去讀函數本身，你會發現他的判斷邏輯似乎說不通，類似下面這句這樣的代碼:
if (((WorkBuffer[BelowOffset] & OrthoEdgeFlag) != 0) && ((WorkBuffer[BelowOffset] & EdgeFlag) != 0))
你會感覺到應該不會出現這種情況啊，我在這裏也看了半天，後來才發現，調用他們的代碼在參數上動了手腳，比如下面這個：
ComputeUpperBounds(ui0, ui1, uh0, uh1,
ColorBuffer,
OffsetLineStart - 1, OffsetLineEnd, SeparationLineLength, StepToPreviousRow, BufferPitch, BufferPitch * Height, EdgeFlag);
我們註意OffsetLineStart - 1這裏的減1，原來他並不是從找到的線條的第一個點開始，而是向前一個點，這樣這個模式就容易理解了。
（2）在MLAABlendHTask以及MLAABlendVTask中，我們註意到都沒有處理最後一行或一列像素，這是因為最後一行或最後一列按照前面的邊緣計算模式，都只可能出現一種模式， 比如，最後一行，其水平邊緣肯定是步成立的。因此，我們在這一行找不到Z或U這種模式，也就沒有必要進行處理。
（3）在Intel的代碼中，有個計算ComputeHc函數：
inline float ComputeHc(UINT* WorkBuffer, UINT PrimaryEdgeLength, UINT OffsetC1, UINT OffsetC2, UINT OffsetD2, UINT OffsetD3)
{
int C2 = SumColor(WorkBuffer[OffsetC2]);
int C1 = SumColor(WorkBuffer[OffsetC1]);
int D3 = SumColor(WorkBuffer[OffsetD3]);
int D2 = SumColor(WorkBuffer[OffsetD2]);
int D3D2Diff = D3 - D2;
int C1C2Diff = C1 - C2;
int D2C2Diff = D2 - C2;
return float(PrimaryEdgeLength * D2C2Diff + C1C2Diff + D3D2Diff) / (PrimaryEdgeLength * (C1C2Diff - D3D2Diff) + C1C2Diff + D3D2Diff);
}
我曾經在修改代碼過程中，把UINT PrimaryEdgeLength改為INT PrimaryEdgeLength，結果對一個測試圖像得到的結果就會發生不錯，如下所示：
第一張為原圖，第二張使用UINT的結果圖，第三張為改為int後的結果圖，很明顯在紅色方框部分的圓弧處依舊有較為明顯的毛刺出現，而整體的代碼只是個數據類型發生了改動。
我曾經弄了很久才發現這個問題出現在那個函數裏，但是發現了後卻一直不知道問題出現在那裏，當我把參數改為int類型，輸出了所有的PrimaryEdgeLength值，確認這裏面確實沒有任何的0或者負數，那按照我的想法不應該會出現這個問題，後來仔細搜索了下，原來問題還是處在數據類型上，因為編譯器在同時存在無符號和有符號數計算時，是會將有符號數強制轉換為無符號數的，即所謂的升格處理，上述代碼裏的D3D2Diff 都有可能有負數存在，因此，這中升格處理會使得計算結果完全偏離了我們的預想。
那在自己想一想，我們確認應該在這裏使用int類型才對，但為什麽我們卻獲得了不理想的效果呢，這裏我認為是原文作者的一個觀點有問題，即HC的計算裏，原文有這句話：
To state the requirements for smooth transition between two shapes, we compute the same blended color twice, using parameters of
each shape:
即要保持平滑，而對黑白圖像，HC就是固定值0.5f.
在Intel的代碼裏，存在好幾處類似這樣的代碼：
if (BelowOffset + StepToNextRow < BufferSizeInPixels)
{
OutH1 = ComputeHc(WorkBuffer, LineLength - NSteps, BelowOffset + StepToNextRow + 1, BelowOffset + 1, BelowOffset, CurrentOffset);
}
else
{ // We are too close to the end of the buffer to be able to do the ComputeHc calculation, so use default hc value of 0.5
OutH1 = 0.5;
}
if ((OutH1 > 0) && (OutH1 < 1))
{ // If the computed value of hc is in the acceptable range, then we‘re done else keep walking
OutS1 = CurrentOffset;
break;
}
即如果通過ComputeHc計算的值不在0和1之間，就取值0.5，那麽前面的使用UINT的情況基本計算出的值都離譜的很，很難在0和1之間，因此，程序最後都相當於直接取Hc等於0.5，而使用int時結果不理想，說明原作者的部分思路也不靠譜。 因此，我建議這部分直接使用0.5，去掉這個ComputeHc的過程，又能節省時間又能還有不錯的效果，不曉得Intel的工程師有沒有註意到這一點。
(4) 這個算法其實本身的計算量是不大的，因為一幅圖像中需要修改的像素其實不多，另外，由於計算特性，基本上支持InPlace操作。
註意到Intel共享的代碼確實還有很多局限性，只能處理32位，只能處理4個倍數大小的圖像，還會破壞原始圖像的Alpha通道，好像還有其他的不好的內存問題（符合前面的條件的圖片運行也會掛），根據算法思路我重構了下代碼，使得其能處理任意大小即灰度和24位圖像。
第一，連續性檢測。我們為了不破壞原始圖像，重新定義一個Width*Height大小的字節空間來保存邊緣信息。對於灰度圖，可用如下代碼搞定。
int BlockSize = 16, Block = (Width - 1) / BlockSize; // 垂直和水平方向比較同步進行，考慮到垂直方向比較時最後一列不能參與 __m128i T = _mm_set1_epi8(Threshold); for (int Y = 0; Y < Height - 1; Y++) // 水平比較時最後一行像素不能參與 { unsigned char *LinePS = Src + Y * Stride; unsigned char *LinePN = LinePS + Stride; unsigned char *LinePM = Mask + Y * Width; for (int X = 0; X < Block * BlockSize; X += BlockSize) { __m128i SrcH1 = _mm_loadu_si128((__m128i *)(LinePS + X)); //__m128i SrcH2 = _mm_loadu_si128((__m128i *)(LinePS + X + 1)); __m128i SrcH2 = _mm_insert_epi8(_mm_srli_si128(SrcH1, 1), LinePS[X + 16], 15); // 似乎速度有沒啥區別 __m128i SrcV1 = _mm_loadu_si128((__m128i *)(LinePN + X)); __m128i AbsDiffH = _mm_absdiff_epu8(SrcH1, SrcV1); // 水平方向比較，abs(Bottom - Current)，似乎和常人了解的水平方向不一致 __m128i AbsDiffV = _mm_absdiff_epu8(SrcH1, SrcH2);; // 垂直方向比較，abs(Right - Current) __m128i FlagH = _mm_and_si128(_mm_set1_epi8(EdgeFlagH), _mm_cmpge_epu8(AbsDiffH, T)); // 設置水平方向的標誌位 __m128i FlagV = _mm_and_si128(_mm_set1_epi8(EdgeFlagV), _mm_cmpge_epu8(AbsDiffV, T)); // 設置垂直方向的標誌位 _mm_storeu_si128((__m128i *)(LinePM + X), _mm_or_si128(FlagH, FlagV)); // 設置總的標誌位 }
// ...............
}
這裏我另外一種方式來做水平和垂直方向的優化，我們一次性只處理一行代碼，但是在垂直方向比較時，單獨往後加載一個像素，註意SrcH1和SrcH2 其實只有一個像素不同，因此，我嘗試用_mm_insert_epi8來減少內存讀取量，但同時增加了一個移位操作，實際測試和直接使用_mm_loadu_si128讀取速度差異基本可忽略。由於使用的時字節變量來保存邊緣信息，因此可使用字節比較_mm_cmpge_epu8的結果來進行位操作。
對於24位或者32位圖像，這時我的處理方式是：
int BlockSize = 4, Block = (Width - 1) / BlockSize; // 垂直和水平方向比較同步進行，考慮到垂直方向比較時最後一列不能參與 __m128i T = _mm_set1_epi8(Threshold); for (int Y = 0; Y < Height - 1; Y++) // 水平比較時最後一行像素不能參與 { unsigned char *LinePS = Src + Y * Stride; unsigned char *LinePN = LinePS + Stride; unsigned char *LinePM = Mask + Y * Width; for (int X = 0; X < Block * BlockSize; X += BlockSize) { __m128i SrcH1 = _mm_loadu_si128((__m128i *)(LinePS + X * 4)); __m128i SrcH2 = _mm_loadu_si128((__m128i *)(LinePS + X * 4 + 4)); __m128i SrcV1 = _mm_loadu_si128((__m128i *)(LinePN + X * 4)); __m128i AbsDiffH = _mm_absdiff_epu8(SrcH1, SrcV1); // 水平方向比較，abs(Bottom - Current)，似乎和常人了解的水平方向不一致 __m128i AbsDiffV = _mm_absdiff_epu8(SrcH1, SrcH2); // 垂直方向比較，abs(Right - Current) __m128i FlagH = _mm_andnot_si128(_mm_cmpeq_epi32(_mm_cmpge_epu8(AbsDiffH, T), _mm_setzero_si128()), _mm_set1_epi32(EdgeFlagH)); // 設置水平方向的標誌位,註意_mm_andnot_si128是對參數的順序敏感的,(~a) & b __m128i FlagV = _mm_andnot_si128(_mm_cmpeq_epi32(_mm_cmpge_epu8(AbsDiffV, T), _mm_setzero_si128()), _mm_set1_epi32(EdgeFlagV)); // 設置垂直方向的標誌位 _mm_storesi128_4char(LinePM + X, _mm_or_si128(FlagH, FlagV));// 設置總的標誌位 }
// ..........
}
其實前面沒提， Intel代碼的Threshold默認設置為16，其實這是個特殊的值，如果要支持任意的閾值，則不能像ComparePixelsSSE那樣做，只有像16，32，64等這樣幾個特殊的值才可以用（相反數高位連續的二進制都為1），要支持任意閾值，我們就可以像上面那樣借助_mm_cmpge_epu8和_mm_cmpeq_epi32來共同實現。
在搜索水平線的時候，也同樣需要更換一種技巧，如下所示：
__m128i PixelFlags = _mm_loadu_si128((__m128i*)(Mask + Y)); // 合理加載數據， 不能超出範圍
int HFlags = _mm_movemask_epi8(_mm_cmpeq_epi8(_mm_and_si128(PixelFlags, Flag), Flag)); // _mm_movemask_epi8提取每個字節的最高位來組成一個16位數
if (HFlags != 0) // 說明在這16個元素裏有部分標記是符合條件的
{
unsigned long Index;
_BitScanForward(&Index, HFlags); // 找到第一個位為1值得索引（從低到高位尋找），intel的bsf指令
OutOffsetLineStart = Y + Index;
goto FindEndOffset; // 尋找結束位置
}
基本的原理和之前的類似，只是要改改相關函數。
結合其他的優化技巧（比如轉置使用SSE處理），目前我做到的速度是：720P的32位圖像大約是8ms，1024P的圖像大約是18ms。以上都是單核的處理結果，當然這個的時間其實還和圖像本身的內容和復雜度有關。
從算法本身的角度來說，是支持並行執行的，如果采用雙線程，應該還有個50%左右的提速，用於遊戲後期的渲染應該是夠了。
算法還有一些有意思的特征，由於識別模式的一些問題，對於同一個圖像內容，旋轉一定角度後，其處理的結果並不一定一樣，比如下圖：
第三圖是對原圖旋轉180度後，進行本算法處理，然後在旋轉180度後得到的結果，很明顯，第三圖的結果和第二圖不一樣，而且要好很多。
對於閾值的選擇，這也是個問題，比如下圖：
中間一幅圖的閾值選擇位16，明顯感覺肩膀和衣服下邊緣還有部分鋸齒，而第三個圖的閾值選擇為32，則鋸齒感又少了很多。
那麽由於算法的內在特性和算法研發的背景，對於兩個方向寬度都大於2個像素的鋸齒，算法是無法解決的，比如下面右側處理後的圖基本無效果。
對於本身就已經光滑的邊緣，這個算法一般也不會產生不好的效果，如下圖所示：
左圖的右半部分本身就是很光滑的，而左半部分又鋸齒，處理後的右圖，左半部分鋸齒基本消失，而右半部分基本沒變。
對於黑白圖像，這個的去鋸齒功能也比較顯著。
那麽我找了幾個真正的遊戲裏的畫面，嘗試了下，後續的去鋸齒效果，確實也能獲得很好的結果：
由於這些圖都是網絡上下載的，一般都為JPG格式，本身經過了壓縮，這樣會引入噪音，會對邊緣的識別有一定的影響，而如果是直接處理顯卡輸出的數據，則不會有任何的損失，因此，理論上會取得更好的效果。
測試見附件的SSE_Optimization_Demo的Other菜單。
Demo下載地址：https://files.cnblogs.com/files/Imageshop/SSE_Optimization_Demo.rar
SSE圖像算法優化系列二十四: 基於形態學的圖像後期抗鋸齒算法--MLAA優化研究。