读源码学算法之poisson matting

Matting中文译为抠图，是图像处理中一个很基本的问题。跟分割不太一样，抠图往往更加精确，甚至连头发丝都要抽取出来，而分割往往很难做到。抠图主要应用在图像合成中，如下所示，抠出的兴趣部分可以合成新的场景。

抠图问题定义为：
$$I=\alpha F+(1-\alpha )B\text{\hspace{1em}(1)}$$
其中， $I$ 表示图像本身，由前景 $F$ 和前景 $B$ 通过透明度参数 $\alpha \in [0,1]$ 合成。
最常见的比如一个玻璃瓶，那么它的颜色基本就是自身的颜色跟背景色的合成。

因此，抠图主要根据输入图像 $I$, 求解右侧的三个未知数$F,B,\alpha$, 很显然这是一个病态问题。
这里主要讲Jian Sun 和 Jiaya Jia发表于SIGGRAPH 2004的经典工作：poisson matting
论文链接：http://kucg.korea.ac.kr/new/seminar/2009/src/PA-09-07-29.pdf
Github代码：https://github.com/xuhongxu96/PoissonMatting

抠图难的部分往往在于前景的边缘，比如人的头发等，除此之外的部分要么是确定的前景，要么是确定的背景，实际上我们只需要对那些不确定的区域求解。因此往往给定一个trimap， 如上图所示，确定的前景用白色表示，确定的背景用黑色表示，剩下的不确定部分使用灰色表示。那么，灰色区域也就是我们需要求解的部分。

首先将(1)看做$F、B、\alpha$ 的函数，那么它的全微分表示为：

$$\nabla I=(F-B)\nabla \alpha +\alpha \nabla F+(1-\alpha )\nabla B\text{\hspace{1em}(2)}$$

假设背景和前景光滑变化，那么：$\nabla F=0,\nabla B=0$， 因此有：$\nabla I=(F-B)\nabla \alpha$, 所以有：

$$\nabla \alpha =\nabla I/(F-B)\text{\hspace{1em}(3)}$$

从而我们可以得到对应的poisson方程：

$$△\alpha =div(\nabla I/(F-B))\text{\hspace{1em}(4)}$$

在图像中，离散的情况下，每个像素仅仅考虑四联通域：

$$△\alpha ={\alpha }_{i,j+1}+{\alpha }_{i,j-1}+{\alpha }_{i+1,j}+{\alpha }_{i-1,j}-4{\alpha }_{i,j}\text{\hspace{1em}(5)}$$

$$di{v}_x(\nabla I/(F-B))=(\nabla I_x(F-B)-\nabla I(F-B{)}_x)/(F-B{)}^2\text{\hspace{1em}(6)}$$

$$di{v}_y(\nabla I/(F-B))=(\nabla I_y(F-B)-\nabla I(F-B{)}_y)/(F-B{)}^2\text{\hspace{1em}(7)}$$

$$div(\nabla I/(F-B))=di{v}_x(\nabla I/(F-B))+di{v}_y(\nabla I/(F-B))\text{\hspace{1em}(8)}$$

于是有：

$$4{\alpha }_{i,j}={\alpha }_{i,j+1}+{\alpha }_{i,j-1}+{\alpha }_{i+1,j}+{\alpha }_{i-1,j}-div(\nabla I/(F-B))\text{\hspace{1em}(9)}$$

接下来就可以迭代求解 ${\alpha }_{i,j}$, 如果其值较大，就将trimap中相应的位置更新为前景，反之更新为背景。

一个问题是为什么不能直接（3）而要转换为（4）所示的poisson方程呢？
我觉得可能的原因是(3)是一个向量函数，(4)是一个标量函数，且（4）的形式方便转换为线性方程组求解。

流程图（参考自github）如下所示：

加注释的源代码如下：

/****************************************************/
// Some useful function
#define I(x, y) (intensity(image(inY(image, y), inX(image, x))))
#define FmB(y, x) (FminusB(inY(FminusB, y), inX(FminusB, x)))

double dist_sqr(cv::Point p1, cv::Point p2) {
    return (p1.x - p2.x) *(p1.x - p2.x) + (p1.y - p2.y)*(p1.y - p2.y);
}
//颜色距离，就是两个颜色的最大值相减
int color_dis(cv::Vec3b p1, cv::Vec3b p2) {
    int t1 = fmax(fmax(p1[0], p1[1]), p1[2]);
    int t2 = fmax(fmax(p2[0], p2[1]), p2[2]);
    return t1 - t2;
}
//防止x,y越界
int inX(cv::Mat &image, int x) {
    if (x < 0) x = 0;
    if (x >= image.cols) x = image.cols - 1;
    return x;
}
int inY(cv::Mat &image, int y) {
    if (y < 0) y = 0;
    if (y >= image.rows) y = image.rows - 1;
    return y;
}
//intensity，三通道的最大值
double intensity(cv::Vec3b v) {  
    return fmax(fmax(v[0], v[1]), v[2]);
}
/***************************************************/

//找出所有边界像素的位置
std::vector<cv::Point> PoissonMatting::findBoundaryPixels(const cv::Mat_<uchar> &trimap, int a, int b){
    std::vector<cv::Point> result;
    for (int x = 1; x < trimap.cols - 1; ++x) {
        for (int y = 1; y < trimap.rows - 1; ++y) {
            if (trimap(y, x) == a) {
                if (trimap(y-1, x)==b || trimap(y+1, x)==b || trimap(y, x-1)==b || trimap(y, x+1)==b) {
                    result.push_back(cv::Point(x, y));
                }
            }
        }
    }
    return result;
}

void PoissonMatting::_matting(cv::Mat _image, cv::Mat _trimap, cv::Mat &_foreground, cv::Mat &_alpha){
    const cv::Mat_<cv::Vec3b> &image = static_cast<const cv::Mat_<cv::Vec3b> &>(_image);
    cv::Mat_<uchar> &trimap = static_cast<cv::Mat_<uchar> &>(_trimap);

    _foreground.create(image.size(), CV_8UC3);
    _alpha.create(image.size(), CV_8UC1);

    cv::Mat_<cv::Vec3b> &foreground = static_cast<cv::Mat_<cv::Vec3b>&>(_foreground);
    cv::Mat_<uchar> &alpha = static_cast<cv::Mat_<uchar>&>(_alpha);

    cv::Mat_<double> FminusB = cv::Mat_<double>::zeros(trimap.rows, trimap.cols);

    for (int times = 0; times < 5; ++times) {
        std::vector<cv::Point> foregroundBoundary = findBoundaryPixels(trimap, 255, 128);
        std::vector<cv::Point> backgroundBoundary = findBoundaryPixels(trimap, 0, 128);

        cv::Mat_<uchar> trimap_blur;
        cv::GaussianBlur(trimap, trimap_blur, cv::Size(9, 9), 0);

        // 构建图像上每个位置的 F-B
        for (int x = 0; x < trimap.cols; ++x) {
            for (int y = 0; y < trimap.rows; ++y) {
                cv::Point current;
                current.x = x;
                current.y = y;
                if (trimap_blur(y, x) == 255) {         //确定的前景部分F-(0,0,0)
                    FminusB(y, x) = color_dis(image(y, x), cv::Vec3b(0, 0, 0));
                } else if (trimap_blur(y, x) == 0) {    //确定的背景部分(0,0,0)-B
                    FminusB(y, x) = color_dis(cv::Vec3b(0, 0, 0), image(y, x));
                } else {
                    // 未知区域的每个位置寻找距离最近的前景像素位置和背景像素位置
                    cv::Point nearestForegroundPoint, nearestBackgroundPoint;
                    double nearestForegroundDistance = 1e9, nearestBackgroundDistance = 1e9;
                    for(cv::Point &p : foregroundBoundary) {
                        double t = dist_sqr(p, current);
                        if (t < nearestForegroundDistance) {
                            nearestForegroundDistance = t;
                            nearestForegroundPoint = p;
                        }
                    }
                    for(cv::Point &p : backgroundBoundary) {
                        double t = dist_sqr(p, current);
                        if (t < nearestBackgroundDistance) {
                            nearestBackgroundDistance = t;
                            nearestBackgroundPoint = p;
                        }
                    }
                    FminusB(y, x) = color_dis(image(nearestForegroundPoint.y, nearestForegroundPoint.x),
                                              image(nearestBackgroundPoint.y, nearestBackgroundPoint.x));
                    if (FminusB(y, x) == 0)
                        FminusB(y, x) = 1e-9;
                }
            }
        }
        // F-B高斯平滑
        cv::GaussianBlur(FminusB, FminusB, cv::Size(9, 9), 0);
        // Solve the Poisson Equation By The Gauss-Seidel Method (Iterative Method)
        for (int times2 = 0; times2 < 300; ++times2) {
            for (int x = 0; x < trimap.cols; ++x) {
                for (int y = 0; y < trimap.rows; ++y) {
                    //白色(F) 黑色(B) 灰色(未知)
                    if (trimap(y, x) == 128) {  
                        // 计算 (▽I/F-B)在所有位置的散度dvgX:(u/v)' = (u'v-uv')/v^2
                        double dvgX = ( (I(x + 1, y) + I(x - 1, y) - 2 * I(x, y)) * FmB(y, x)
                                - (I(x + 1, y) - I(x, y)) * (FmB(y, x + 1) - FmB(y, x)) )
                                / (FmB(y, x) * FmB(y, x));
                        // 计算 (▽I/F-B)在所有位置的散度dvgY:(u/v)' = (u'v-uv')/v^2
                        double dvgY = ( (I(x, y + 1) + I(x, y - 1) - 2 * I(x, y)) * FmB(y, x)
                                - (I(x, y + 1) - I(x, y)) * (FmB(y + 1, x) - FmB(y, x)) )
                                / (FmB(y, x) * FmB(y, x));
                        double dvg = dvgX + dvgY;
                        // a的散度为：a(i+1,j)+a(i-1,j)+a(i,j+1)+a(i,j-1)-4a(i,j) = dvg, 因此得到下面的式子
                        double newAlpha = ((double)alpha(y, x + 1)
                                        + alpha(y, x - 1)
                                        + alpha(y + 1, x)
                                        + alpha(y - 1, x)
                                        - dvg * 255.0) / 4.0;
                        // 根据得到的alpha更新Trimap
                        if (newAlpha > 253)     trimap(y, x) = 255;
                        else if (newAlpha < 3)  trimap(y, x) = 0;
                        if (newAlpha < 0)       newAlpha = 0;
                        if (newAlpha > 255)     newAlpha = 255;
                        // 更新alpha
                        alpha(y, x) = newAlpha;
                    } 
                    else if (trimap(y, x) == 255)  alpha(y, x) = 255;   //前景区域不需计算
                    else if (trimap(y, x) == 0)    alpha(y, x) = 0;     //背景区域不需计算
                }
            }
        }
    }
    // 得到了alpha,因此可以合成红色背景的新图像 I = a Image + (1-a)Red ???
    for (int x = 0; x < alpha.cols; ++x) {
        for (int y = 0; y < alpha.rows; ++y) {
            foreground(y, x) = ((double) alpha(y, x) / 255) * image(y, x) + ((255.0 - alpha(y, x)) / 255 * cv::Vec3b(0, 0, 255));
        }
    }
}

另一篇Matting的经典论文：A Closed Form Solution to Natural Image Matting
的解析参见以下资源：

1. 知乎：Closed Form Solution to Natural Image Matting公式推导 写得很清晰。
2. 知乎：Closed Form Matting 算法抠图以及优化（1）——目标函数
3. 知乎：Closed Form Matting 算法抠图以及优化（2）——拉普拉斯矩阵
4. 知乎：Closed Form Matting 算法抠图以及优化（3）——问题的求解
5. 知乎：Closed Form Matting 算法抠图以及优化 (4)——代码实现
6. CMU老师的纸上推导也很赞！