Demosaic，中文直接發翻譯為去馬賽克， 但是其本質上並不是去除馬賽克，這讓很多第一次接觸這個名詞或者概念的人總會想的更多。因此，一般感測器在採集資訊時一個位置只採集RGB顏色中的一個通道，這樣可以減少採集量，降低成本，由於人的視覺對綠色最為敏感，所以一般情況下，綠色分量會比紅色和藍色分量多一倍的資訊，這樣，根據RGB不同的位置排布，就產生了RGGB、GBRG、GRBG、BGGR四種常用的模式，比如RGGB模式就如下所示：
RGGB BGGR GBRG GRBG
去馬賽克的目的就是從這些確實的資訊中儘量的完美的復原原始資訊，比如上圖第一個點，只有紅色分量是準確的，那就要想辦法透過其他的資訊重構出改點的綠色和藍色分量。
目前，關於這方面的資料也是很多的，這裡我描述下目前我已經最佳化和處理的四個演算法，並分享部分程式碼。
一、雙線性處理
一個最為直接和簡單的想法就是利用領域的相關資訊透過插值來彌補當前缺失的顏色分量，我們以RGGB格式為例，參考上圖（座標都是從0開始，X軸從左向右遞增，Y軸從上向下遞增）。
先不考慮邊緣。
比如（1,1）座標這個點，現在已經有了B這個分量，缺少R和G這個分量，但是他四個對角線上都有R分量，因此，可以使用這個4個畫素的平均值來計算改點的R分量，而在其上下左右四個方向有恰好有4個點的G分量，同理可以用這個四個的平均值來評估當前點的G值。
在考慮（1,2）這個座標點，他的當前有效值是G分量，缺少R和B分量，而他則沒有（1,1）點那麼幸運，他周邊的有效R和B分量都只有2個，因此，只能利用這兩個有效值的平均值來評估該點的R/B分量。
其他的點類似處理。
在考慮邊緣，由於邊緣處的畫素總有某一個方向或某二個方向缺少畫素，因此，可以利用映象的關係把缺少的那一邊取映象的位置資訊來補充。
一個簡單的高效的C++程式碼如下所示：
int IM_DecodeBayerRG8ToBGR_Bilinear_PureC(unsigned char* Src, BitmapData Dest)
{
int Status = IM_STATUS_OK;
int Width = Dest.Width, Height = Dest.Height, Stride = Dest.Stride;
if (((Width & 1) == 1) || ((Height & 1) == 1)) return IM_STATUS_INVALIDPARAMETER;
if ((Width < 8) || (Height < 8)) return IM_STATUS_INVALIDPARAMETER;
unsigned char* RowCopy = (unsigned char*)malloc((Width + 2) * 3 * sizeof(unsigned char));
if (RowCopy == NULL) return IM_STATUS_OUTOFMEMORY;
unsigned char* First = RowCopy;
unsigned char* Second = RowCopy + (Width + 2);
unsigned char* Third = RowCopy + (Width + 2) * 2;
unsigned char* SrcP = (unsigned char*)Src;
Second[0] = SrcP[0];
memcpy(Second + 1, SrcP, Width * sizeof(unsigned char));
Second[Width + 1] = SrcP[Width - 1];
memcpy(First, Second, (Width + 2) * sizeof(unsigned char)); // 第一行和第二行一樣
Third[0] = SrcP[Width];
memcpy(Third + 1, SrcP + Width, Width * sizeof(unsigned char));
Third[Width + 1] = SrcP[Width + Width - 1];
for (int Y = 0; Y < Height; Y++)
{
unsigned char* LinePD = Dest.Scan0 + Y * Stride;
if (Y != 0)
{
unsigned char* Temp = First; First = Second; Second = Third; Third = Temp;
}
if (Y == Height - 1)
{
memcpy(Third, Second, (Width + 2) * sizeof(unsigned char));
}
else
{
Third[0] = SrcP[(Y + 1) * Width];
memcpy(Third + 1, SrcP + (Y + 1) * Width, Width * sizeof(unsigned char));
Third[Width + 1] = SrcP[(Y + 1) * Width + Width - 1];
}
if ((Y & 1) == 0) // 偶數列
{
for (int X = 0; X < Width; X++, LinePD += 3)
{
int P0 = First[X], P1 = First[X + 1], P2 = First[X + 2];
int P3 = Second[X], P4 = Second[X + 1], P5 = Second[X + 2];
int P6 = Third[X], P7 = Third[X + 1], P8 = Third[X + 2];
if ((X & 1) == 0)
{
LinePD[0] = (P0 + P2 + P6 + P8 + 2) >> 2;
LinePD[1] = (P1 + P3 + P5 + P7 + 2) >> 2;
LinePD[2] = P4;
}
else
{
LinePD[0] = (P1 + P7 + 1) >> 1;
LinePD[1] = P4;
LinePD[2] = (P3 + P5 + 1) >> 1;
}
}
}
else
{
for (int X = 0; X < Width; X++, LinePD += 3)
{
int P0 = First[X], P1 = First[X + 1], P2 = First[X + 2];
int P3 = Second[X], P4 = Second[X + 1], P5 = Second[X + 2];
int P6 = Third[X], P7 = Third[X + 1], P8 = Third[X + 2];
if ((X & 1) == 0)
{
LinePD[0] = (P3 + P5 + 1) >> 1;
LinePD[1] = P4;
LinePD[2] = (P1 + P7 + 1) >> 1;
}
else
{
LinePD[0] = P4;
LinePD[1] = (P1 + P3 + P5 + P7 + 2) >> 2;
LinePD[2] = (P0 + P2 + P6 + P8 + 2) >> 2;
}
}
}
}
free(RowCopy);
return IM_STATUS_OK;
}
這個演算法是非常高效的，而且極易使用指令集最佳化，在一臺普通的配置的PC上（12th Gen Intel(R) Core(TM) i7-12700F 2.10 GHz）處理1902*1080的圖，大概只需要2ms，SSE最佳化後可以做到0.6ms，但是這個方法忽略邊緣結構和通道間的相關性，從而容易導致顏色偽彩和影像模糊，比如下面這個經典的測試圖，解碼後的結果有點慘不忍睹。
用彩色顯示RGGB格式圖（即把其他通道的顏色分量設計為0） 雙線性解碼後的結果
在解碼後的水平和垂直方向柵欄中都可以看到明顯的色彩異常。
但是由於這個演算法的極度高效性，在有些要求不高的場合，依舊可以使用該演算法。
二、Hamilton-Adams演算法
這也是個非常經典的演算法，在說明這個演算法之前，必須說明下色差恆定理論。
色差恆定準則與色比恆定準則都是基於顏色通道之間的相關性，目的都是把顏色通道之間的相關性資訊引入顏色插值演算法，提高插值的準確性。色差相比於色比有兩點優勢：
第一，色差的運算簡單，更容易實現。第二， 色比在G通道接近0時誤差較大，色差不存在這類問題。因此，絕大多數顏色插值演算法中使用了色差。
那麼色差恆定理論，其最為核心的思想就是：臨近的兩個彩色畫素，其相同顏色分量之間的差異值應該近似差不多，用公式表示如下：
R(X,Y) - R(X-1,Y) = G(X,Y) - G(X-1,Y) = B(X,Y) - B(X-1,Y)
那這是時候如果我們已經獲取了某個顏色通道的所有位置的值，透過這個公式就很容易推匯出其他位置缺失的值了。
我們還是以上面的（1,1）座標點為例，假定我們已經獲取了所有G通道的資料，也就是說這個（1,1）這個點實際只有R通道資料缺失了（B資料本身就有），這個時候根據顏色恆差理論，應該有
R(1,1) - R(0,0) = G(1,1) - G(0,0) ------> R(1,1) = G(1,1) + R(0,0) - G(0, 0)
實際上滿足這個等式還有（0,2)、(2,0)、（2,2）這三個點（這三個點的紅色分量是準確值），所以為了得到更好的精度，我們可以透過下式最終確定R(1,1)的值。
R(1,1) = (G(1,1) + R(0,0) - G(0, 0) + G(1,1) + R(0,2) - G(0, 2) + G(1,1) + R(2,0) - G(2, 0) + G(1,1) + R(2,2) - G(2, 2)) /4
整理後即為：
R(1,1) = G(1,1) + (R(0,0) + R(0,2) + R(2, 0) + R(2,2) - G(0,0) -G(0,2) -G(2, 0) -G(2,2)) / 4
對於（1，2）這個點，G通道資料本來就有，缺少R和B,那根據顏色恆差理論，應該有
R(1,2) - G(1,2) = R(0,2) - G(0,2)
B(1,2) - G(1,2) = B(1,1) - G(1,1)
同樣的道理，還有
R(1,2) - G(1,2) = R(2,2) - G(2,2)
B(1,2) - G(1,2) = B(1,3) - G(1,3)
類似（1,2）的方法，把他們綜合起來可以得到更為精確的結果：
R(1,2) = G(1,2) + ((R(0,2) + R(2,2) - G(0,2) - G(2,2)) / 2
B(1,2) = G(1,2) + ((B(1,1) + B(1,3) - G(0,2) - G(2,2)) / 2
以上利用顏色恆差理論，就把各個通道之間的資料關聯了起來。那麼前面的計算都有一個前提條件，就是綠色通道的資料都已經知道了。
我們前面說過，綠色通道的資料量本身只缺少了一半，而缺少了那一半中任何一個點的資料都可以用周邊四個點的領域資料填充，因此，如果我們綠色通道就這樣處理，而紅色和藍色則用顏色恆差理論藉助綠色通道的資料後結果如何呢，我們就把這樣做的演算法叫做CDCT把，一個簡單的嚴重這個演算法的程式碼如下所示：
int IM_DecodeBayerRG8ToBGR_CDCT_PureC(unsigned char* Src, BitmapData Dest)
{
int Status = IM_STATUS_OK;
int Width = Dest.Width, Height = Dest.Height, Stride = Dest.Stride;
if (((Width & 1) == 1) || ((Height & 1) == 1)) return IM_STATUS_INVALIDPARAMETER;
if ((Width < 8) || (Height < 8)) return IM_STATUS_INVALIDPARAMETER;
unsigned char* Blue = (unsigned char*)malloc(Width * Height * sizeof(unsigned char));
unsigned char* Green = (unsigned char*)malloc(Width * Height * sizeof(unsigned char));
unsigned char* Red = (unsigned char*)malloc(Width * Height * sizeof(unsigned char));
if ((Blue == NULL) || (Green == NULL) || (Red == NULL))
{
Status = IM_STATUS_OUTOFMEMORY;
goto FreeMemory;
}
// 先直接複製資料，也無需把通道的值單獨提取出來填充，因為後續會把無效的點填充的
memcpy(Blue, Src, Height * Width * sizeof(unsigned char));
memcpy(Green, Src, Height * Width * sizeof(unsigned char));
memcpy(Red, Src, Height * Width * sizeof(unsigned char));
// 因為Green分量佔了1/2畫素，先填充Green畫素
// 因為後續的Green分量涉及到了3*3的領域，對於邊緣的部分，直接使用四個點的平均，單獨提出來,這部分計算量很小，無需加速
for (int X = 0; X < Width; X++)
{
IM_CalcBorderGreen_CDCT(Green, Width, Height, X, 0);
IM_CalcBorderGreen_CDCT(Green, Width, Height, X, Height - 1);
}
for (int Y = 1; Y < Height - 1; Y++)
{
IM_CalcBorderGreen_CDCT(Green, Width, Height, 0, Y);
IM_CalcBorderGreen_CDCT(Green, Width, Height, Width - 1, Y);
}
// 填充Green通道的無效位置的資料
for (int Y = 1; Y < Height - 1; Y++)
{
int Index = Y * Width;
for (int X = 1; X < Width - 1; X++)
{
// 偶數行和偶數列 或者奇數行和奇數列，綠色分量都是無效的
if (((X + Y) & 1) == 0)
{
Green[Index + X] = (Green[Index + X - Width] + Green[Index + X + Width] + Green[Index + X - 1] + Green[Index + X + 1] + 2) / 4;
}
}
}
IM_RGGB_CalcRed_CDCT_PureC(Red, Green, Width, Height);
IM_RGGB_CalcBlue_CDCT_PureC(Blue, Green, Width, Height);
Status = IM_CombineRGB_PureC(Blue, Green, Red, Dest.Scan0, Dest.Width, Dest.Height, Dest.Stride, Width);
if (Status != IM_STATUS_OK) goto FreeMemory;
FreeMemory:
if (Blue != NULL) free(Blue);
if (Green != NULL) free(Green);
if (Red != NULL) free(Red);
return Status;
}
其中 IM_RGGB_CalcRed_CDCT_PureC程式碼如下所示：
void IM_RGGB_CalcRed_CDCT_PureC(unsigned char* Red, unsigned char* Green, int Width, int Height)
{
// R G R G R G
// G B G B G B
// R G R G R G
// G B G B G B
// R G R G R G
// G B G B G B
// 色差恆定原理：色差是色度訊號（R和B分量）與亮度訊號（G分量）的差，在影像很小的範圍內，當前畫素的色差與其周圍點的色差是差不多的，也就是類似下面的說法
// R(I,J)-G(I,J) = R(I,J + 1)-G(I,J + 1),或者寫成R(I,J)-R(I,J+1) = G(I,J)-G(I,J + 1)。這樣利用已經前面已經完成重構的G分量，可以重構出其他未知的R和B分量。
for (int X = 0; X < Width; X++)
{
IM_RGGB_CalcBorderRed_CDCT(Red, Green, Width, Height, X, 0);
IM_RGGB_CalcBorderRed_CDCT(Red, Green, Width, Height, X, Height - 1);
}
for (int Y = 1; Y < Height - 1; Y++)
{
IM_RGGB_CalcBorderRed_CDCT(Red, Green, Width, Height, 0, Y);
IM_RGGB_CalcBorderRed_CDCT(Red, Green, Width, Height, Width - 1, Y);
}
// 填充Red通道的無效位置的資料
for (int Y = 1; Y < Height - 1; Y++)
{
int Index = Y * Width;
for (int X = 1; X < Width - 1; X++)
{
// 偶數行奇數列， 水平方向填充紅色分量
if ((Y & 1) == 0 && (X & 1) == 1)
{
Red[Index + X] = IM_ClampToByte((Red[Index + X - 1] + Red[Index + X + 1] - Green[Index + X - 1] - Green[Index + X + 1] + 1) / 2 + Green[Index + X]);
}
// 奇數行偶數列， 垂直方向填充紅色分量
else if ((Y & 1) == 1 && (X & 1) == 0)
{
Red[Index + X] = IM_ClampToByte((Red[Index + X - Width] + Red[Index + X + Width] - Green[Index + X - Width] - Green[Index + X + Width] + 1) / 2 + Green[Index + X]);
}
// 奇數行奇數列， 水平垂直方向填充紅色分量
else if ((Y & 1) == 1 && (X & 1) == 1)
{
Red[Index + X] = IM_ClampToByte((Red[Index + X - Width - 1] + Red[Index + X - Width + 1] + Red[Index + X + Width - 1] + Red[Index + X + Width + 1] - Green[Index + X - Width - 1] - Green[Index + X - Width + 1] - Green[Index + X + Width - 1] - Green[Index + X + Width + 1] + 2) / 4 + Green[Index + X]);
}
}
}
}
IM_RGGB_CalcBlue_CDCT_PureC是類似的道理。
經過測試，這樣做的結果和直接雙線性相比，基本沒有什麼差異的，所以直接這樣還是不行的，
Hamilton-Adams等人在結合綠色通道的一階導數和周邊紅色和藍色的通道的二階倒數的基礎上，對綠色分量的插值提出瞭如下演算法，這個過程考慮到了畫素之間的邊緣資訊：
當水平方向的梯度大於垂直方向的梯度時，使用垂直方向的有關畫素計算結果，否則使用水平方形的有關值，如果兩者相等，則使用平均值。
實際操作時，都會定義一個閾值，如果水平和垂直的梯度之差的絕對值小於閾值，則使用平均值，如果在閾值之外，再考慮水平和垂直之間的關係。這樣能獲得更為合理的結果。
實際上，仔細看看，上面每個方向的G5的計算也是利用到了顏色恆差理論的，只是只是單獨利用了水平或者垂直方向的畫素而已。
我們分享一下更具上述思路編寫的C++程式碼結果：
int IM_DecodeBayerRG8ToBGR_HamiltonAdams_PureC(unsigned char* Src, BitmapData Dest, int Threshold)
{
int Status = IM_STATUS_OK;
int Width = Dest.Width, Height = Dest.Height, Stride = Dest.Stride;
// 寬度和高度都必須是偶數
if (((Width & 1) == 1) || ((Height & 1) == 1)) return IM_STATUS_INVALIDPARAMETER;
if ((Width < 8) || (Height < 8)) return IM_STATUS_INVALIDPARAMETER;
unsigned char* Blue = (unsigned char*)malloc(Width * Height * sizeof(unsigned char));
unsigned char* Green = (unsigned char*)malloc(Width * Height * sizeof(unsigned char));
unsigned char* Red = (unsigned char*)malloc(Width * Height * sizeof(unsigned char));
if ((Blue == NULL) || (Green == NULL) || (Red == NULL))
{
Status = IM_STATUS_OUTOFMEMORY;
goto FreeMemory;
}
// 先直接複製資料，也無需把通道的值單獨提取出來填充，因為後續會把無效的點填充的
memcpy(Blue, Src, Height * Width * sizeof(unsigned char));
memcpy(Green, Src, Height * Width * sizeof(unsigned char));
memcpy(Red, Src, Height * Width * sizeof(unsigned char));
// 因為Green分量佔了1/2畫素，先填充Green畫素
// 因為後續的Green分量涉及到了5*5的領域，對於邊緣的部分，直接使用四個點的平均，單獨提出來,這部分計算量很小，無需加速
for (int X = 0; X < Width; X++)
{
IM_CalcBorderGreen_CDCT(Green, Width, Height, X, 0);
IM_CalcBorderGreen_CDCT(Green, Width, Height, X, 1);
IM_CalcBorderGreen_CDCT(Green, Width, Height, X, Height - 2);
IM_CalcBorderGreen_CDCT(Green, Width, Height, X, Height - 1);
}
for (int Y = 2; Y < Height - 2; Y++)
{
IM_CalcBorderGreen_CDCT(Green, Width, Height, 0, Y);
IM_CalcBorderGreen_CDCT(Green, Width, Height, 1, Y);
IM_CalcBorderGreen_CDCT(Green, Width, Height, Width - 2, Y);
IM_CalcBorderGreen_CDCT(Green, Width, Height, Width - 1, Y);
}
// 處理剩下的能否安全訪問領域的演算法
for (int Y = 2; Y < Height - 2; Y++)
{
int IndexC = Y * Width;
int IndexN1 = (Y + 1) * Width, IndexN2 = (Y + 2) * Width;
int IndexP1 = (Y - 1) * Width, IndexP2 = (Y - 2) * Width;
unsigned char* Sample = (Y & 1) == 0 ? Red : Blue;
for (int X = 2; X < Width - 2; X++)
{
// 只有當X和Y都是偶數或者都是奇數時才需要處理
if (((X + Y) & 1) == 0)
{
// 周邊藍色或者紅色分量的二階導數
int SecDH = 2 * Sample[IndexC + X] - Sample[IndexC + X + 2] - Sample[IndexC + X - 2];
int SecDV = 2 * Sample[IndexC + X] - Sample[IndexN2 + X] - Sample[IndexP2 + X];
// 加上綠色分量的一階導數得到梯度
int GradH = IM_Abs(Green[IndexC + X - 1] - Green[IndexC + X + 1]) + IM_Abs(SecDH);
int GradV = IM_Abs(Green[IndexP1 + X] - Green[IndexN1 + X]) + IM_Abs(SecDV);
// 如果垂直或者水平的梯度差不多，則計算周邊的平均值
if (IM_Abs(GradV - GradH) < Threshold)
Green[IndexC + X] = IM_ClampToByte((Green[IndexP1 + X] + Green[IndexN1 + X] + Green[IndexC + X - 1] + Green[IndexC + X + 1] + 2) / 4 + (SecDH + SecDV + 4) / 8);
// 如果水平差異小一些，則利用水平方向的平均值
else if (GradH < GradV)
Green[IndexC + X] = IM_ClampToByte((Green[IndexC + X - 1] + Green[IndexC + X + 1] + 1) / 2 + (SecDH + 2) / 4);
// 否則利用垂直方向的平均值
else
Green[IndexC + X] = IM_ClampToByte((Green[IndexP1 + X] + Green[IndexN1 + X] + 1) / 2 + (SecDV + 2) / 4);
}
}
}
IM_RGGB_CalcRed_CDCT_SSE(Red, Green, Width, Height);
IM_RGGB_CalcBlue_CDCT_SSE(Blue, Green, Width, Height);
Status = IM_CombineRGB_PureC(Blue, Green, Red, Dest.Scan0, Dest.Width, Dest.Height, Dest.Stride, Width);
if (Status != IM_STATUS_OK) goto FreeMemory;
FreeMemory:
if (Blue != NULL) free(Blue);
if (Green != NULL) free(Green);
if (Red != NULL) free(Red);
return Status;
}
我們分享下前面說的CDCT演算法以及HamiltonAdams演算法的結果：
簡單的CDCT結果 HamiltonAdams結果
從上面右圖的結果可以看到，相比於雙線性，水平柵欄處似乎已經看不到色彩異常了，垂直柵欄處的異常點也大為減少，但是還是存在些許瑕疵。
在速度方面，由於直接使用雙線性時，可以直接把資料資料有規律的填充到目的圖中，而且計算量很簡單，而使用色差恆定理論後，由於順序的要求以及一些編碼方面的原因，需要使用一些中間記憶體，而且計算量相對來說大了很多，因此，速度也慢了不少，上述C++的程式碼處理1080P的影像，需要大概7ms，經過SSE最佳化後的程式碼可以達到4ms左右的速度，這個速度在某些實時性要求比較高的場景下還是具有實際價值的。
為了進一步提高效果，在HA演算法的基礎上，後續還是有不少人提出了更多的改進演算法，在 IPOL上，可以找到一篇 A Mathematical Analysis and Implementation of Residual Interpolation Demosaicking Algorithms 的綜述性文章，文章裡列舉了數種型別的處理方式，那個網站是個也提供了線上的DEMO和程式碼，可以看到各自的處理後的結果,如下圖所示：
不過這些都是比較傳統的演算法方面的東西了，而且我看了程式碼和有關效果，感覺這些演算法更偏重於理論，實際操作起來，可能效率完全不夠，同時IPOL上還提供了一篇基於CNN深度學習方面的演算法，效果我看了下，確實還是很不錯的，不過就是感覺速度比較堪憂，有興趣可以在A Study of Two CNN Demosaicking Algorithms這裡瀏覽。
那麼我後面關注到的是IPOL的另外一篇文章: Zhang-Wu Directional LMMSE Image Demosaicking,透過測試發現這個文章的效果還是非常不錯和穩定的，而且分析其程式碼，感覺可以做很大的最佳化工作。
這個文章對R和B通道的處理方式是和HA演算法一樣的，都是先要獲取G通道的全部資料，然後採用色差恆定原理計算R和B，在計算G通道時，從水平和垂直方向計算 (G-R) 和 (G-B) 的顏色差異開始。然後，這些計算被視為對實際色差的噪聲估計，並使用線性最小均方誤差框架將它們最佳化組合。
所謂的最佳化組合可以這樣理解，傳統的HA演算法，在計算G通道時，就是根據水平和垂直方向的梯度，決定最後是使用水平還是垂直方向的計算值，而ZhangWu演算法則不完全使用水平或垂直的資訊，而且根據某些過程計算出水平和垂直方向各自的權重然後融合。
這裡簡單的描述了一下色差噪聲估計的公式：
這幾個公式可以看成是一個簡單的去燥過程，其中f表示原始的噪音資訊，s呢就是一個對f的簡單的低通濾波，比如高斯濾波，然後呢公式9計算s的一定領域大小的平均值， 公式10計算s在這個領域內的方差， 公式11則計算f和s差異平方的均值。
利用以上資訊則可以估算出去除噪音後的估算值 u（公式12）以及對一個的估算誤差（公式13）。
分別對水平和垂直反向上的顏色差異進行這樣的濾波，然後利用上面的結果按照下述公式進行融合：
上面的公式是一個一維的去燥過程，我有空將其改為2維的圖形去燥看看效果如何。
這些公式在論文的配套程式碼裡都有有關的實現，有興趣的朋友可以有去看看程式碼，我將其過程進行過程化，這個函式大概有如下程式碼構成：
// Zhang–Wu Directional LMMSE Image Demosaicking int IM_DecodeBayerRG8ToBGR_ZhangWu_SSE(unsigned char* Src, BitmapData Dest) { int Status = IM_STATUS_OK; int Width = Dest.Width, Height = Dest.Height, Stride = Dest.Stride; // 寬度和高度都必須是偶數 if (((Width & 1) == 1) || ((Height & 1) == 1)) return IM_STATUS_INVALIDPARAMETER; // 寬度或者高度小於8有些領域會溢位 if ((Width < 8) || (Height < 8)) return IM_STATUS_INVALIDPARAMETER; unsigned char* Blue = (unsigned char*)malloc(Width * Height * sizeof(unsigned char)); unsigned char* Green = (unsigned char*)malloc(Width * Height * sizeof(unsigned char)); unsigned char* Red = (unsigned char*)malloc(Width * Height * sizeof(unsigned char)); short* DiffH = (short *)malloc(Width * Height * sizeof(short)); short* DiffV = (short *)malloc(Width * Height * sizeof(short)); short* LowPassH = (short*)malloc(Width * Height * sizeof(short)); short* LowPassV = (short *)malloc(Width * Height * sizeof(short)); if ((Blue == NULL) || (Green == NULL) || (Red == NULL) || (DiffH == NULL) || (DiffV == NULL) || (LowPassH == NULL) || (LowPassV == NULL)) { Status = IM_STATUS_OUTOFMEMORY; goto FreeMemory; } // 先直接複製資料，也無需把通道的值單獨提取出來填充，因為後續會把無效的點填充的 memcpy(Blue, Src, Height * Width * sizeof(unsigned char)); //memcpy(Green, Src, Height * Width * sizeof(unsigned char)); // 綠色通道因為其特殊性，後續會在程式碼裡進行填充，不需要單獨複製資料的 memcpy(Red, Src, Height * Width * sizeof(unsigned char)); // 獲取水平方向上準確訊號和插值訊號的差異 IM_GetHoriDiffSignal_SSE(Src, DiffH, Width, Height); // 獲取垂直方向上準確訊號和插值訊號的差異 IM_GetVertDiffSignal_SSE(Src, DiffV, Width, Height); // 對水平差異進行1*9的高斯濾波 IM_HoriLowPass1X9_SSE(DiffH, LowPassH, Width, Height); // 對垂直差異進行9*1的高斯濾波 IM_VertLowPass9X1_SSE(DiffV, LowPassV, Width, Height); // 透過LMMSE演算法計算完整的綠色通道 Status = IM_RGGB_CalcGreen_ZW_SSE(Src, Green, DiffH, DiffV, LowPassH, LowPassV, Width, Height, 4); if (Status != IM_STATUS_OK) goto FreeMemory; // 使用色差恆定原理計算出Red通道的資料 IM_RGGB_CalcRed_CDCT_SSE(Red, Green, Width, Height); // 使用色差恆定原理計算出Blue通道的資料 IM_RGGB_CalcBlue_CDCT_SSE(Blue, Green, Width, Height); // 把RGB單通道資料組合成RGB彩色影像 Status = IM_CombineRGB_SSE(Blue, Green, Red, Dest.Scan0, Dest.Width, Dest.Height, Dest.Stride, Width); if (Status != IM_STATUS_OK) goto FreeMemory; FreeMemory: if (Blue != NULL) free(Blue); if (Green != NULL) free(Green); if (Red != NULL) free(Red); if (DiffH != NULL) free(DiffH); if (DiffV != NULL) free(DiffV); if (LowPassH != NULL) free(LowPassH); if (LowPassV != NULL) free(LowPassV); return Status; }
具體的實現就不做過多的探討，原始作者的C程式碼效率非常的低下，簡單的測試了下1080P的圖大概需要1分鐘左右，這完全沒有實際的意義的，所以需要進行深度的最佳化，比如水平方向的濾波，原作者採用1*9的濾波器，我稍作改進並用SSE指令最佳化如下：
// 不支援In-Place操作 void IM_HoriLowPass1X9_SSE(short* Src, short* Dest, int Width, int Height) { for (int Y = 0; Y < Height; Y++) { int Index = Y * Width; // 邊緣採用映象的關係 4 3 2 1 0 1 2 3 4 5 6 7 Dest[Index + 0] = ((Src[Index + 4] + Src[Index + 4]) * 4 + (Src[Index + 3] + Src[Index + 3]) * 8 + (Src[Index + 2] + Src[Index + 2]) * 16 + (Src[Index + 1] + Src[Index + 1]) * 23 + Src[Index + 0] * 26 + 64) / 128; Dest[Index + 1] = ((Src[Index + 3] + Src[Index + 5]) * 4 + (Src[Index + 2] + Src[Index + 4]) * 8 + (Src[Index + 1] + Src[Index + 3]) * 16 + (Src[Index + 0] + Src[Index + 2]) * 23 + Src[Index + 1] * 26 + 64) / 128; Dest[Index + 2] = ((Src[Index + 2] + Src[Index + 6]) * 4 + (Src[Index + 1] + Src[Index + 5]) * 8 + (Src[Index + 0] + Src[Index + 4]) * 16 + (Src[Index + 1] + Src[Index + 3]) * 23 + Src[Index + 2] * 26 + 64) / 128; Dest[Index + 3] = ((Src[Index + 1] + Src[Index + 7]) * 4 + (Src[Index + 0] + Src[Index + 6]) * 8 + (Src[Index + 1] + Src[Index + 5]) * 16 + (Src[Index + 2] + Src[Index + 4]) * 23 + Src[Index + 3] * 26 + 64) / 128; // W-8 W-7 W-6 W-5 W-4 W-3 W-2 W-1 W-2 W-3 W-4 W-5 Dest[Index + Width - 4] = ((Src[Index + Width - 8] + Src[Index + Width - 2]) * 4 + (Src[Index + Width - 7] + Src[Index + Width - 1]) * 8 + (Src[Index + Width - 6] + Src[Index + Width - 2]) * 16 + (Src[Index + Width - 5] + Src[Index + Width - 3]) * 23 + Src[Index + Width - 4] * 26 + 64) / 128; Dest[Index + Width - 3] = ((Src[Index + Width - 7] + Src[Index + Width - 3]) * 4 + (Src[Index + Width - 6] + Src[Index + Width - 2]) * 8 + (Src[Index + Width - 5] + Src[Index + Width - 2]) * 16 + (Src[Index + Width - 4] + Src[Index + Width - 2]) * 23 + Src[Index + Width - 3] * 26 + 64) / 128; Dest[Index + Width - 2] = ((Src[Index + Width - 6] + Src[Index + Width - 4]) * 4 + (Src[Index + Width - 5] + Src[Index + Width - 3]) * 8 + (Src[Index + Width - 4] + Src[Index + Width - 3]) * 16 + (Src[Index + Width - 3] + Src[Index + Width - 1]) * 23 + Src[Index + Width - 2] * 26 + 64) / 128; Dest[Index + Width - 1] = ((Src[Index + Width - 5] + Src[Index + Width - 5]) * 4 + (Src[Index + Width - 4] + Src[Index + Width - 4]) * 8 + (Src[Index + Width - 3] + Src[Index + Width - 4]) * 16 + (Src[Index + Width - 2] + Src[Index + Width - 2]) * 23 + Src[Index + Width - 1] * 26 + 64) / 128; int BlockSize = 8, Block = (Width - 8) / BlockSize; for (int X = 4; X < 4 + Block * BlockSize; X += BlockSize) { // (V0 * 4 + V1 * 8 + V2 * 16 + V3 * 23 + V4 * 26 + V5 * 23 + V6 * 16 + V7 * 8 + V8 * 4 + 64) / 128 __m128i V0 = _mm_loadu_si128((__m128i*)(Src + Index + X - 4)); __m128i V1 = _mm_loadu_si128((__m128i*)(Src + Index + X - 3)); __m128i V2 = _mm_loadu_si128((__m128i*)(Src + Index + X - 2)); __m128i V3 = _mm_loadu_si128((__m128i*)(Src + Index + X - 1)); __m128i V4 = _mm_loadu_si128((__m128i*)(Src + Index + X + 0)); __m128i V5 = _mm_loadu_si128((__m128i*)(Src + Index + X + 1)); __m128i V6 = _mm_loadu_si128((__m128i*)(Src + Index + X + 2)); __m128i V7 = _mm_loadu_si128((__m128i*)(Src + Index + X + 3)); __m128i V8 = _mm_loadu_si128((__m128i*)(Src + Index + X + 4)); __m128i V08 = _mm_slli_epi16(_mm_add_epi16(V0, V8), 2); // (V0 + V8) * 4 __m128i V17 = _mm_slli_epi16(_mm_add_epi16(V1, V7), 3); // (V1 + V7) * 8 __m128i V26 = _mm_slli_epi16(_mm_add_epi16(V2, V6), 4); // (V2 + V6) * 16 __m128i V35 = _mm_mullo_epi16(_mm_add_epi16(V3, V5), _mm_set1_epi16(23)); // (V3 + V5) * 23 __m128i V44 = _mm_mullo_epi16(V4, _mm_set1_epi16(26)); // V4 * 26 __m128i Sum = _mm_add_epi16(_mm_add_epi16(_mm_add_epi16(V08, V17), _mm_add_epi16(V26, V35)), V44); __m128i Avg = _mm_srai_epi16(_mm_add_epi16(Sum, _mm_set1_epi16(64)), 7); _mm_storeu_si128((__m128i*)(Dest + Index + X), Avg); } for (int X = 4 + Block * BlockSize; X < Width - 4; X++) { Dest[Index + X] = ((Src[Index + X - 4] + Src[Index + X + 4]) * 4 + (Src[Index + X - 3] + Src[Index + X + 3]) * 8 + (Src[Index + X - 2] + Src[Index + X + 2]) * 16 + (Src[Index + X - 1] + Src[Index + X + 1]) * 23 + Src[Index + X] * 26 + 64) / 128; } } }
這極大的提高了執行速度。
在LMMSE的計算過程中，作者的M大小取的是4，即涉及到的領域大小也是9個畫素，作者在程式碼裡使用硬迴圈實現，實際上這個也就是普通一進一出的統計，完全可以做點最佳化的，特別是垂直方向的迴圈，每次都要跳一行畫素，cachemiss很大，所以透過適當的改動結構，能極大的提高速度。經過最佳化，我們測試這個演算法處理處理一副1080P的影像，SSE版本大概耗時12ms，我自己最佳化後的C++程式碼耗時約27ms(使用了編譯器自己的向量化方式編譯，而且非原始作者的C++程式碼），這個速度應該說在實時系統中還是可以接受的，而且這個過程還是可以進行多執行緒並行化的，如果開兩個執行緒，SSE版本有望做到8ms一幀。
同HA演算法相比，這個演算法得到的結果更加完美，而且有瑕疵的地方更少，如下所示：
Zhang Wu演算法結果 原始影像
我們再分享一組測試及影像：
Mosaic Bilinear
HA Zhang Wu
整體來看，效果最好是的Zhang Wu演算法， 其次是HA， 最一般的就是Bilinear了，但是，也不要一板子拍死Bilinear， 在很多不是很複雜的場景的情況下，使用他得到的效果依舊是可以滿足需求的，關鍵是他真的快啊。
當然，如果對演算法執行效率和效果都要做一個均衡的話，應該說HA演算法就比較適中了。
為了比較方便，我編寫了一個簡易的DEMO，供大家做演算法驗證。
https://files.cnblogs.com/files/Imageshop/DeMosaic.rar?t=1725238639&download=true
如果想時刻關注本人的最新文章，也可關注公眾號：
參考資料：
// https://blog.csdn.net/OdellSwan/article/details/136887148 ISP演算法 | Demosaic(一)
// https://blog.csdn.net/OdellSwan/article/details/137246117 ISP演算法 | Demosaic(二)
// https://cloud.tencent.com/developer/article/1934157 ISP影像處理之Demosaic演算法及相關
// https://zhuanlan.zhihu.com/p/594341024 Demosaic（二）Hamilton & Adams插值演算法
// https://zhuanlan.zhihu.com/p/144651850 Understanding ISP Pipeline - Demosaicking
// https://blog.csdn.net/feiyanjia/article/details/124366793