/* Cumulative image of the multiplication of p1.*p2.
p1 and p2 can be the same pointer and it becomes square cumulative image.
The returned cumulative image MUST be freed by the caller! */
template <class T1, class T2, class T3>
BOOL MultiplyCumImage(T1 *p1, T2 *p2, T3 **ppCum)
{
long i, j;
long Width = GetWidth();
long Height = GetHeight();
long integral_len = (Width + 1)*(Height + 1);
// Only allocate cumulative image memory when *ppCum is NULL
if (*ppCum == NULL)
{
try { *ppCum = new T3[integral_len]; }
catch (CException *pe)
{
pe->Delete();
return FALSE;
}
}
memset(*ppCum, 0, sizeof(T3)*integral_len);
// The cumulative values of the leftmost and the topmost pixels are always 0.
for (i = 1; i <= Height; i++)
{
T3 *prow, *puprow;
prow = *ppCum + i*(Width + 1);
puprow = *ppCum + (i - 1)*(Width + 1);
T3 sum = 0;
long up_row_idx = (i - 1)*Width;
for (j = 1; j <= Width; j++)
{
long idx = up_row_idx + j - 1;
sum += p1[idx] * p2[idx];
prow[j] = puprow[j] + sum;
}
}
return TRUE;
}
这样导向滤波实现的主要问题都解决了,算法步骤如下:
计算引导图G的积分图和平方积分图
计算输入图像P的积分图, P和G的乘积积分图
用上两步得出的积分图计算P和G的均值,G的方差,P和G的协方差,窗口半径为r
然后用(5)(6)式计算系数图A和B
计算A和B的积分图
计算A和B的窗口半径r的均值,并用(8)式计算输出图Q
下面的参考代码中,pData存储输入和输出图像,pGuidedData引导图,radius领域半径
long len = Width*Height;
// Cululative image and square cululative for guided image G
UINT64 *pCum = NULL, *pSquareCum = NULL;
CumImage(pGuidedData, &pCum);
MultiplyCumImage(pGuidedData, pGuidedData, &pSquareCum);
// Allocate memory for a and b
float *pa, *pb;
pa = new float[len];
pb = new float[len];
memset(pa, 0, sizeof(float)*len);
memset(pb, 0, sizeof(float)*len);
UINT64 *pInputCum = NULL, *pGPCum = NULL;
CumImage(pData, &pInputCum);
MultiplyCumImage(pGuidedData, pData, &pGPCum);
int field_size;
UINT64 cum, square_cum;
long uprow, downrow, upidx, downidx; // In cumulative image
long leftcol, rightcol;
float g_mean, g_var; // mean and variance of guided image
long row_idx = 0;
UINT64 p_cum, gp_cum;
float p_mean;
// Calculate a and b
// Since we're going to calculate cumulative image of a and b, we have to calculate the whole image of a and b.
for (i = 0; i < Height; i++)
{
// Check the boundary for radius
if (i < radius) uprow = 0;
else uprow = i - radius;
upidx = uprow*(Width + 1);
if (i + radius >= Height) downrow = Height;
else downrow = i + radius + 1;
downidx = downrow*(Width + 1);
for (j = 0; j < Width; j++)
{
// Check the boundary for radius
if (j < radius) leftcol = 0;
else leftcol = j - radius;
if (j + radius >= Width) rightcol = Width;
else rightcol = j + radius + 1;
field_size = (downrow - uprow)*(rightcol - leftcol);
long p1, p2, p3, p4;
p1 = downidx + rightcol;
p2 = downidx + leftcol;
p3 = upidx + rightcol;
p4 = upidx + leftcol;
// Guided image summary in the field
cum = pCum[p1] - pCum[p2] - pCum[p3] + pCum[p4];
// Guided image square summary in the field
square_cum = pSquareCum[p1] - pSquareCum[p2] - pSquareCum[p3] + pSquareCum[p4];
// Field mean
g_mean = (float)(cum) / field_size;
// Field variance
g_var = float(square_cum) / field_size - g_mean * g_mean;
// Summary of input image in the field
p_cum = pInputCum[p1] - pInputCum[p2] - pInputCum[p3] + pInputCum[p4];
// Input image field mean
p_mean = float(p_cum) / field_size;
// Multiply summary in the field
gp_cum = pGPCum[p1] - pGPCum[p2] - pGPCum[p3] + pGPCum[p4];
long idx = row_idx + j;
pa[idx] = (float(gp_cum) / field_size - g_mean*p_mean) / (g_var + epsilon);
pb[idx] = p_mean - g_mean*pa[idx];
}
row_idx += Width;
}
// not needed after this
delete[] pCum;
delete[] pSquareCum;
delete[] pInputCum;
delete[] pGPCum;
// Cumulative image of a and b
float *pCuma = NULL, *pCumb = NULL;
CumImage(pa, &pCuma);
CumImage(pb, &pCumb);
// Finally calculate the output image q=ag+b
float mean_a, mean_b;
row_idx = Hstart*Width;
for (i = Hstart; i < Hend; i++)
{
// Check the boundary for radius
if (i < radius) uprow = 0;
else uprow = i - radius;
upidx = uprow*(Width + 1);
if (i + radius >= Height) downrow = Height;
else downrow = i + radius + 1;
downidx = downrow*(Width + 1);
for (j = Wstart; j < Wend; j++)
{
// Check the boundary for radius
if (j < radius) leftcol = 0;
else leftcol = j - radius;
if (j + radius >= Width) rightcol = Width;
else rightcol = j + radius + 1;
field_size = (downrow - uprow)*(rightcol - leftcol);
long p1, p2, p3, p4;
p1 = downidx + rightcol;
p2 = downidx + leftcol;
p3 = upidx + rightcol;
p4 = upidx + leftcol;
// Field mean
mean_a = (pCuma[p1] - pCuma[p2] - pCuma[p3] + pCuma[p4]) / field_size;
// Field mean
mean_b = (pCumb[p1] - pCumb[p2] - pCumb[p3] + pCumb[p4]) / field_size;
// New pixel value
long idx = row_idx + j;
int value = int(mean_a*pGuidedData[idx] + mean_b);
CLAMP0255(value);
pData[idx] = value;
}
row_idx += Width;
}
delete[] pa;
delete[] pb;
delete[] pCuma;
delete[] pCumb;