基於像素的傳統圖像修補算法實現

本文轉載自查看原文 2021-04-25 15:04 449

圖像修補算法可以用來修復圖像中的瑕疵，划痕等，或者移除不需要的內容，比如水印或其他物體。傳統的圖像修補算法有基於像素和基於區域的兩種分類，本文介紹基於像素的傳統圖像修補算法的實現。論文可以在這里找到An image inpainting technique based on the fast marching method。基於像素的圖像修補算法主要步驟包括兩個，一是如何在要修不的圖像區域中確定下一個要修補的像素，即要通過一種方法找到下一個待修補的像素，其周圍有最多的已知像素，保證最好的修補效果。二是如何用周圍的已知像素來估算要修補的像素。

找到要修補的下一個像素的基本思想是從需要進行修補區域的邊界開始，由邊界到中心逐漸填充待修補區域中的所有像素。越靠近邊界的像素自然周圍有越多的已知像素，而越向修補區域中心去，周圍已知像素就越少。當然修補區域的形狀是不確定的，所謂的中心是遠離邊界，向修補區域內部。完成這一步驟的就是所謂的FMM(Fast Marching Method)算法。記待修補的圖像區域為Ω，∂Ω為Ω的邊界。為修補Ω中的每一個像素，從∂Ω的初始位置∂Ω_i開始，向Ω內部前進。下一個帶修補像素的選擇，以diatance的大小順序，由小到大選擇。這里distance指的是distance map的distance，所有已知像素的distance為0。如果我們得到一張包含待修補區域的distance map，distance的大小代表了待修補像素距離已知像素的遠近。以distance的升序排列，即可以得到修補像素的順序。

通常要計算distance map的計算量較大，也有一些近似的快速算法。而FMM算法提供了一種即時計算distance的方法。FMM維護一個所謂的narrow band的像素集合，它就是邊界∂Ω的所有像素。對每一個像素，保存它的distance值T，它的像素值I，還有一個flag值f。f可以是以下三種

BAND：屬於narrow band的像素
KNOWN：在∂Ω外部的像素，已知像素。T，I已知
INSIDE：在∂Ω內部的像素，像素待修補。T，I未知

FMM初始時，對於所有位於邊界∂Ω以及其外部區域的像素賦值T為0，而其內部待修補像素的T值為某一個大數（實際使用1e6）。f值按照上述定義給所有圖像像素賦值。所有BAND像素放入一個以T值升序排列的有序集合，就是所謂的narrow band。因為對這個集合主要涉及到插入刪除操作，實際中我使用的鏈表。

運行時按照下面的偽代碼。需要注意的是，不知道什么原因，論文中有一些我認為明顯錯漏的地方，特別是在偽代碼中。所以我這里的偽代碼和原文中的有所不同。

while (NarrowBand not empty) { extract P(i,j) = head(NarrowBand); f(i,j) = KNOWN; for (k,l) in (i,j-1), (i-1,j), (i+1,j), (i,j+1) { if (f(k,l) == INSIDE) { T(k,l) = min(solve(k-1,l,k,l-1), solve(k+1,l,k,l-1), solve(k-1,l,k,l+1), solve(k+1,l,k,l+1)); Inpaint(k,l); f(k,l) = BAND; Insert (k,l) in NarrowBand; } } }

原文中這段偽代碼最大的問題是對f(k,l) != KNOWN的判斷。不是KNOWN的話，就是BAND或者INSIDE，那么就可能會重復把已經是BAND的像素重新加入到NarrowBand了。所以這里直接改成判斷f(k,l)是否是INSIDE，只有INSIDE的像素才需要更新T並修補。FMM首先從NarrowBand中取出一個有最小T值的BAND像素，修改其flag為KNOWN。然后檢查其臨近的上下左右四聯通像素，如果是INSIDE，則依次做

更新其distance值T。
實際修補其像素值。
修改其flag為BAND，加入NarrowBand。

這幾步做完，實際上把邊界∂Ω向邊界Ω內部前進了一步。注意這里我是先做了更新T值，是因為在實際inpaint像素值時是有可能需要用到它本身T值的（后面會講到）。而在原文中，更新distance值T是放在inpaint之后，那么inpaint時就無法用了。這也是我無法理解，認為原文這里是錯誤的原因。

這里更新(k,l)的T值需要其上下左右四聯通像素的T值，解如下方程

其中，，對y也類似。如果兩個max值都不為0的話，設其中一個max值中的已知量，比如T(i-1,j)為T₁，另一個max值中的已知量比如T(i,j+1)為T₂，上述方程實際是一個一元二次方程，很容易求得其兩個解為。而如果要求解同時大於T₁和T₂，則只有一個解，並且abs(T₁-T₂)<1。所以這里也不需要像原文中的偽代碼來。解得四個值，取最小值。比較簡單，就不寫偽代碼了，直接看后面solve()函數的源代碼。

這樣FMM算法的主要流程就完成了，重復以上過程直到NarrowBand為空，就是從邊界一直深入到Ω內部，到所有像素都修補完畢。

剩下的問題就是如何修補單個像素。待修補像素由其鄰域中所有已知像素的加權和得到。

p為要inpaint的像素，為p的大小為ε的領域。領域中單個像素q對p的貢獻為，為q點的梯度。

w(p,q)為權重，是三個分量的乘積。

其中dir(p,q)是方向分量，其保證越接近distance變化的法線方向（即T的梯度方向）的像素對inpaint像素的貢獻越大。這也是前面FMM中為什么要把對T的更新放在實際inpaint之前的原因。雖然計算T的中央差分可以不用它本身的T值，但是如果其他相鄰的像素有一個是INSIDE的話，那么就需要用到它本身的值了。dst(p,q)也很容易理解，是地理距離分量，q距離p越近，對p的影響越大。lev(p,q)是T等值集合分量，讓越接近通過p的輪廓等值線的像素對inpaint貢獻更大，也就是q距離已知像素和p距離已知像素的差值越小越好，即兩者值越接近越好。在實際應用中，d₀和T₀都使用1。這三個分量結合，共同決定領域內每個像素像素對修補值的影響。

下面給出Inpaint的代碼。PriorityHeap是一個以T值排序的像素集合，內部使用雙向鏈表，方便快速插入刪除。初始化narrow band，用同樣的narrow band對三通道數據做Inpaint。可以對三通道各用一個線程加快處理速度。

enum { KNOWN, BAND, INSIDE, }; typedef struct PaintType { int flag; float T; } PaintType; typedef struct HeapElement { float T; CPoint point; HeapElement(float T, CPoint& p) { this->T = T; point = p; } } HeapElement; class PriorityHeap { protected: list<HeapElement> m_band; public: PriorityHeap() {} PriorityHeap(PriorityHeap& ph) { m_band = ph.m_band; } PriorityHeap(PriorityHeap&& ph) noexcept { m_band = move(ph.m_band); } void Push(float T, CPoint& p) { auto rit = m_band.rbegin(); for (; rit != m_band.rend(); ++rit) if (rit->T <= T) break; m_band.emplace(rit.base(), T, p); } void Push(float T, CPoint&& p) { auto rit = m_band.rbegin(); for (; rit != m_band.rend(); ++rit) if (rit->T <= T) break; m_band.emplace(rit.base(), T, p); } BOOL Pop(CPoint* pp) { if (m_band.size() == 0) return FALSE; *pp = m_band.front().point; m_band.pop_front(); return TRUE; } size_t size() { return m_band.size(); } }; BOOL Inpaint(int radius) { CRect* pRect = nullptr; CRect Mask; // sMask包含要修補的區域
    unique_ptr<BYTE[]> sMask; if (pDibMask) { GetMaskCircumRect(Mask); GetMask(sMask); pRect = &Mask; } else
        return FALSE; long i, j; long Width = GetWidth(); long Height = GetHeight(); long Hstart, Hend, Wstart, Wend; long len; pRect->NormalizeRect(); if (pRect->top < 0) pRect->top = 0; if (pRect->left < 0) pRect->left = 0; if (pRect->bottom > Height) pRect->bottom = Height; if (pRect->right > Width) pRect->right = Width; Hstart = pRect->top; Hend = pRect->bottom; Wstart = pRect->left; Wend = pRect->right; ASSERT(Hstart >= 0 && Hstart <= Hend && Hend <= Height); ASSERT(Wstart >= 0 && Wstart <= Wend && Wend <= Width); unique_ptr<BYTE[]> pRed, pGreen, pBlue; if (GetRGB(pRed, pGreen, pBlue, &len) == FALSE) return FALSE; // initialize narrow band
    unique_ptr<PaintType[]> sPT(new PaintType[len]); PriorityHeap narrow_band; BYTE* pm = sMask.get(); PaintType* ppt = sPT.get(); for (i = 0; i < Hstart - 1; ++i)        // to (Hstart-1) because 1 pixel extension to pRect needs to be checked
 { for (j = 0; j < Width; ++j) { // known
            ppt->flag = KNOWN; ppt->T = 0; ++ppt; } pm += Width; } for (; i < min(Hend + 1, Height); ++i) { for (j = 0; j < Wstart - 1; ++j) { // known
            ppt->flag = KNOWN; ppt->T = 0; ++ppt; ++pm; } for (; j < min(Wend + 1, Width); ++j) { if (*pm == 0) { // check 4 connected pixels
                if (i - 1 >= 0 && *(pm - Width)) { // band
                    ppt->flag = BAND; ppt->T = 0; narrow_band.Push(0, CPoint(j, i)); ++ppt; ++pm; continue; } else if (j - 1 >= 0 && *(pm - 1)) { // band
                    ppt->flag = BAND; ppt->T = 0; narrow_band.Push(0, CPoint(j, i)); ++ppt; ++pm; continue; } else if (j + 1 < Width && *(pm + 1)) { // band
                    ppt->flag = BAND; ppt->T = 0; narrow_band.Push(0, CPoint(j, i)); ++ppt; ++pm; continue; } else if (i + 1 < Height && *(pm + Width)) { // band
                    ppt->flag = BAND; ppt->T = 0; narrow_band.Push(0, CPoint(j, i)); ++ppt; ++pm; continue; } else { // known
                    ppt->flag = KNOWN; ppt->T = 0; ++ppt; ++pm; } } else { // inside
                ppt->flag = INSIDE; ppt->T = 1e6f; ++ppt; ++pm; } } for (; j < Width; ++j) { // known
            ppt->flag = KNOWN; ppt->T = 0; ++ppt; ++pm; } } for (; i < Height; ++i) { for (j = 0; j < Width; ++j) { // known
            ppt->flag = KNOWN; ppt->T = 0; ++ppt; } } PriorityHeap narrow_band2(narrow_band), narrow_band3(narrow_band); unique_ptr<PaintType[]> sPT2, sPT3; sPT2.reset(new PaintType[len]); memcpy(sPT2.get(), sPT.get(), sizeof(PaintType) * len); sPT3.reset(new PaintType[len]); memcpy(sPT3.get(), sPT.get(), sizeof(PaintType) * len); #if 0 Inpaint(pRed.get(), sPT.get(), narrow_band, radius); Inpaint(pGreen.get(), sPT2.get(), narrow_band2, radius); Inpaint(pBlue.get(), sPT3.get(), narrow_band3, radius); #else auto InpaintFunc = bind((void(*)(BYTE*, PaintType*, PriorityHeap&, const int)) & Inpaint, placeholders::_1, placeholders::_2, placeholders::_3, placeholders::_4); auto thread2 = thread(InpaintFunc, pGreen.get(), sPT2.get(), narrow_band2, radius); auto thread3 = thread(InpaintFunc, pBlue.get(), sPT3.get(), narrow_band3, radius); Inpaint(pRed.get(), sPT.get(), narrow_band, radius); thread2.join(); thread3.join(); #endif PutRGB(nullptr, pRed, pGreen, pBlue); return TRUE; } static float solve(const PaintType* pPT, const long offset1, const long offset2) { float sol = 1e6; float T1, T2; T1 = pPT[offset1].T; T2 = pPT[offset2].T; if (pPT[offset1].flag != INSIDE) { if (pPT[offset2].flag != INSIDE) { if (abs(T1 - T2) >= 1.0f) sol = 1 + min(T1, T2); else sol = (T1 + T2 + sqrtf(2 - (T1 - T2) * (T1 - T2))) / 2; } else sol = 1 + T1; } else if (pPT[offset2].flag != INSIDE) sol = 1 + T2; return sol; } void Inpaint(BYTE* pData, PaintType* pPT, PriorityHeap& narrow_band, const int radius) { long Width = GetWidth(); long Height = GetHeight(); long len = Width * Height; unique_ptr<float[]> pfData(new float[len]);        // 保存中間數據，提高精度
    for (long k = 0; k < len; ++k) pfData[k] = pData[k]; long offset; CPoint point; while (narrow_band.Pop(&point)) { offset = point.y * Width + point.x; pPT[offset].flag = KNOWN; for (int ii = 0; ii < 4; ++ii) { long off_i; CPoint point_i(point); switch (ii) { case 0: if (--point_i.y < 0) continue; off_i = offset - Width; break; case 1: if (--point_i.x < 0) continue; off_i = offset - 1; break; case 2: if (++point_i.x >= Width) continue; off_i = offset + 1; break; case 3: if (++point_i.y >= Height) continue; off_i = offset + Width; break; } int& flag = pPT[off_i].flag; if (flag == INSIDE) { bool off_ivalid0 = true, off_ivalid1 = true, off_ivalid2 = true, off_ivalid3 = true; long off_i0 = off_i - Width; if (point_i.y - 1 < 0) off_ivalid0 = false; long off_i1 = off_i - 1; if (point_i.x - 1 < 0) off_ivalid1 = false; long off_i2 = off_i + 1; if (point_i.x + 1 >= Width) off_ivalid2 = false; long off_i3 = off_i + Width; if (point_i.y + 1 >= Height) off_ivalid3 = false; // calculate distance map T
 { float minimum = 1e6; if (off_ivalid1 && off_ivalid0) minimum = min(minimum, solve(pPT, off_i1, off_i0)); if (off_ivalid2 && off_ivalid0) minimum = min(minimum, solve(pPT, off_i2, off_i0)); if (off_ivalid1 && off_ivalid3) minimum = min(minimum, solve(pPT, off_i1, off_i3)); if (off_ivalid2 && off_ivalid3) minimum = min(minimum, solve(pPT, off_i2, off_i3)); pPT[off_i].T = minimum; } // inpaint pixel
 { long i, j; long hstart = point_i.y - radius; long hend = point_i.y + radius; long wstart = point_i.x - radius; long wend = point_i.x + radius; if (hstart < 0) hstart = 0; if (hend > Height - 1) hend = Height - 1; if (wstart < 0) wstart = 0; if (wend > Width - 1) wend = Width - 1; // gradT of the pixel being inpainted
                    float gradTx = 0, gradTy = 0; { if (off_ivalid2 && pPT[off_i2].flag != INSIDE) { if (off_ivalid1 && pPT[off_i1].flag != INSIDE) gradTx = (pPT[off_i2].T - pPT[off_i1].T) * 0.5f; else gradTx = pPT[off_i2].T - pPT[off_i].T; } else { if (off_ivalid1 && pPT[off_i1].flag != INSIDE) gradTx = pPT[off_i].T - pPT[off_i1].T; } if (off_ivalid3 && pPT[off_i3].flag != INSIDE) { if (off_ivalid0 && pPT[off_i0].flag != INSIDE) gradTy = (pPT[off_i3].T - pPT[off_i0].T) * 0.5f; else gradTy = pPT[off_i3].T - pPT[off_i].T; } else { if (off_ivalid0 && pPT[off_i0].flag != INSIDE) gradTy = pPT[off_i].T - pPT[off_i0].T; } } float s_weight = 0, sum = 0; for (i = hstart; i <= hend; ++i) { long row_offset = i * Width; float r_y = float(point_i.y - i); for (j = wstart; j <= wend; ++j) { long coff = row_offset + j; if (coff == off_i) continue;        // skip the pixel being inpainted
                            float r_x = float(point_i.x - j); if (pPT[coff].flag != INSIDE) //if (pPT[coff].flag == KNOWN)
 { float dir = abs(gradTx * r_x + gradTy * r_y); if (dir < 1e-6f) dir = 1e-6f; float dst = 1.0f / (r_y * r_y + r_x * r_x); float lev = 1.0f / (1 + abs(pPT[coff].T - pPT[off_i].T)); float w = dir * dst * lev; float gradIx = 0, gradIy = 0; if (j + 1 < Width && pPT[coff + 1].flag != INSIDE) { if (j - 1 >= 0 && pPT[coff - 1].flag != INSIDE) gradIx = (pfData[coff + 1] - pfData[coff - 1]) * 0.5f; else gradIx = float(pfData[coff + 1] - pfData[coff]); } else { if (j - 1 >= 0 && pPT[coff - 1].flag != INSIDE) gradIx = float(pfData[coff] - pfData[coff - 1]); } if (i + 1 < Height && pPT[coff + Width].flag != INSIDE) { if (i - 1 >= 0 && pPT[coff - Width].flag != INSIDE) gradIy = (pfData[coff + Width] - pfData[coff - Width]) * 0.5f; else gradIy = float(pfData[coff + Width] - pfData[coff]); } else { if (i - 1 >= 0 && pPT[coff - Width].flag != INSIDE) gradIy = float(pfData[coff] - pfData[coff - Width]); } float kx, ky; kx = gradIx * r_x; ky = gradIy * r_y; sum += w * (pfData[coff] + 0.35f * (kx + ky)); s_weight += w; } } } pfData[off_i] = sum / s_weight; CLAMP0255(pfData[off_i]); pData[off_i] = BYTE(roundf(pfData[off_i])); } flag = BAND; // push into narrow band
 narrow_band.Push(pPT[off_i].T, point_i); } } } }