影像修補演算法可以用來修復影像中的瑕疵,劃痕等,或者移除不需要的內容,比如水印或其他物體,傳統的影像修補演算法有基于像素和基于區域的兩種分類,本文介紹基于像素的傳統影像修補演算法的實作,論文可以在這里找到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)為T1,另一個max值中的已知量比如T(i,j+1)為T2,上述方程實際是一個一元二次方程,很容易求得其兩個解為
,而如果要求解同時大于T1和T2,則只有一個解
,并且abs(T1-T2)<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距離已知像素的差值越小越好,即兩者值越接近越好,在實際應用中,d0和T0都使用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);
}
}
}
}
Inpaint
代碼中一個比較繁瑣的地方是對邊界的判斷,另外像素加權和那里,我加了一個0.35的系數,不加的話,比較容易資料溢位,
下面來看一下效果,首先用論文里一個簡單的三單色圖案來試一下,以白色圓環為要修補的區域,修補半徑分別為1和3的結果,可以看到半徑越大,模糊就越嚴重一些,對于這種純色圖案,半徑越小越好,

原圖(640x480) ε=1 ε=3
再看一些例子,

原圖(440x198) 待修補圖

ε=3 ε=5
下面這個特意把要修補的某些區域加粗了,可以看到,inpaint區域厚度越大,最后修補區域的結果越模糊,

原圖(423x435) 待修補圖

ε=3 ε=5
也可以用來去除某些不需要的區域,

原圖(698x319) inpaint模板

ε=3 ε=5
總體來說,基于像素的影像修補方法必然地會讓修補結果模糊,因此它適用于比較小區域的影像修補,背景影像也不能太復雜,否則,修復存在模糊及填充紋理不自然的問題,另外待修補區域的選擇也比較重要,區域的形狀對最終結果的影響也很明顯,

原圖(2048x1269) inpaint模板
ε=3 ε=5
轉載請註明出處,本文鏈接:https://www.uj5u.com/qita/280174.html
標籤:其他
上一篇:這些愛好可以提高編程技能

