Patchmatch算法简单实现

最新推荐文章于 2024-08-08 08:08:28 发布

咸鱼半条

最新推荐文章于 2024-08-08 08:08:28 发布

阅读量8.3k

点赞数 8

分类专栏：图像处理

图像处理专栏收录该内容

55 篇文章 7 订阅

订阅专栏

本文是一篇使用OpenCV实现patchmatch算法的笔记，实现环境为Visual Studio 2015 x64 + OpenCV 3.0。笔者并没有很多CV方面的经验，所以对算法的解释可能使用不规范的词汇，望各位指正谅解。

patchmatch算法的核心目的是在两张图片之间快速寻找对应的小区域。在应用上，patchmatch算法可以结合图像重组等技术，来实现诸如图像修复、图片融合、去水印等功能。Adobe Photoshop CS5之后提供的修复画刷、内容感知移动功能就是基于patchmatch算法实现的。

详情可戳 http://v.youku.com/v_show/id_XMTIxOTgyNTcy.html
墙外可戳 https://www.youtube.com/watch?v=fMe19oTz6vk

1.Patchmatch的概念

1.1 我们要做什么

patchmatch中的“patch”指的就是图片的一个小区域。具体在此算法中，一个patch是图片的一个3x3或5x5的小正方形区域。patchmatch顾名思义也就是对两张图片的对应patch完成匹配。以下我们称最终要进行处理的图为“目标图”，提供完整信息的图为“参考图”，并假设目标图与参考图的长宽是一致的。

1.2 相似度的定义

如何算是一个成功的匹配呢？很简单，对于图片A的一个patch，我们只需要找到图片B中与之“距离”最近/“相似度”最高的patch就可以了。所以这里我们需要对patch的距离/相似度加以定义。

在这个Demo中笔者使用负的“绝对距离”（各项直接相减、取绝对值、求和）作为相似性量度。

double sim_abs_diff(Mat& a, Mat& b, int _) {    
    // 此处int是为了函数签名一致性加的，实际无用
    Mat tmp;
    absdiff(a, b, tmp); 
    return 1e4 - sum(sum(tmp))[0];
}

1.3 偏移量

由于参考图与目标图的长宽是一致的，所以两个点a，b之间的偏移量就等于b-a（注意，a和b是两个向量）。对于目标图中所有的patch区域，我们要在参考图中找一个对应的patch匹配。我们使用一个三维数组记录这一匹配，数组中存储的是目标图到参考图各个patch的偏移量。为了避免多维数组传参的麻烦事，我们使用标准库中的vector来完成。

typedef vector< vector < vector <int> > > Vector3i;  // 定义三维vector类型
Vector3i f; f.resize(n_rows);
for(int i=0;i<n_rows;i++){
    f[i].resize(n_cols);
    for(int j=0;j<n_cols;j++) f[i][j].resize(2);
}

例如，目标图中左上角为(30,45)的patch在参考图中的对应patch左上角为(10,10)，则f[30][45] = {-20, -35}

2.PatchMatch算法流程

算法共分为初始化、迭代、图片重建三个过程。

2.1 初始化

初始化过程就是把映射表中的每个元素设置成一个合法随机值，即随机的建立一个目标图到参考图的映射关系。

2.2 迭代

迭代过程分为“传播”(Propagation)和“随机搜索”(Random Search)两个步骤。

2.2.1 传播

以下我们把目标图中左上角坐标为(x, y)的patch记作patch(x, y)，并记其在参考图中对应的patch为match(x, y)。即 f : patch(x, y) -> match(x, y)。

传播过程中，我们按扫描序（从左到右，从上到下）遍历每一个patch，将此patch与match(x, y)，match(x-1, y)，match(x, y-1)做相似性比较，取得最相似的match(x’, y’)，并更新f[y][x] = f[y’][x’]。此更新公式最大的优势在于，一旦一个patch被“正确”匹配，则几乎可以确定右下方的所有像素都可以被传播成“正确”的。

另外，传播过程在奇数次迭代时使用扫描序，偶数次迭代使用逆扫描序，注意代码中的循环语句。

for (int i = row_start; i != row_end; i += step) {
    for (int j = col_start; j != col_end; j += step) {
        float sm[3];

        /* 取得当前正在处理的patch(r, c) */
        /* pick_patch(image, row_topleft, col_topleft, row_offset, col_offset, patch_size) */
        Mat patch = pick_patch(img_dst(Range(i, i + patch_size), Range(j, j + patch_size));

        /* 取得当前映射的match(r, c) */
        Mat ipatch = pick_patch(img_ref, i, j, 0, 0, patch_size);
        sm[0] = sim(patch, ipatch, 1);

        /* 取得r方向和c方向的“邻接”match(r-1, c)和match(r, c-1) */
        /* checkvalid函数用于检查新的偏移量是否合法 */
        if (checkvalid(i, j, f[i - step][j][0], f[i - step][j][1])) {
            Mat xpatch = pick_patch(img_ref, i, j, f[i - step][j][0], f[i - step][j][1], patch_size);
            sm[1] = sim(patch, xpatch, 1);
        } else sm[1] = -1e6f;

        if (checkvalid(i, j, f[i][j - step][0], f[i][j - step][1])) {
            Mat ypatch = pick_patch(img_ref, i, j, f[i][j - step][0], f[i][j - step][1], patch_size);
            sm[2] = sim(patch, ypatch, 1);
        } else sm[2] = -1e6f;

        /* 取得相似度最大的match坐标 */
        int k = util::argmax(sm, 3); // 取argmax
        v[i][j] = sm[k];

        /* 更新映射表 */
        switch (k) {
        case 0: break;
        case 1: f[i][j][0] = f[i - step][j][0]; f[i][j][1] = f[i - step][j][1]; break;
        case 2: f[i][j][0] = f[i][j - step][0]; f[i][j][1] = f[i][j - step][1]; break;
        default: break; // error
        }
    }
}

2.2.2 随机搜索

仅仅与邻接match比较很难得到长距离的正确映射，我们还需要进行一个空间稍大的搜索操作以提供更多的匹配可能。为了提高效率，patchmatch算法使用了随机搜索技术：

按扫描序遍历(x, y)：
对于patch(x, y)，我们首先给出一个和图片尺寸一致的搜索窗口（即搜索半径）。设f[x][y] = (x’, y’)，则我们在参考图上以(x’, y’)为中心，以搜索半径为半边长的矩形内随机取一点(x”, y”)，并以其对应的区域作为可能的match’(x, y)。将patch和match’的相似度与patch和原match的相似度相比，如果match’有更高的相似度，则更新f[x][y]为(x”, y”)；否则不做任何事。每次操作结束后，将搜索半径缩小一半，重复上述过程直至搜索半径不大于1。
以下代码是对patch(x, y)的随机搜索过程。

int r_ws = n_rows, c_ws = n_cols;
float alpha = 0.5f, exp = 0.5f;    // alpha是当前的搜索半径因子，exp是收缩速度
                                   // 窗口尺寸 = min(r_ws*alpha, c_ws*alpha)
while (r_ws*alpha > 1 && c_ws*alpha > 1) {
    /* 确定合法的搜索空间，避免下标越界 */
    int rmin = util::max(0, int(i + f[i][j][0] - r_ws*alpha));
    int rmax = util::min(int(i + f[i][j][0] + r_ws*alpha), n_rows - patch_size);
    int cmin = util::max(0, int(j + f[i][j][1] - c_ws*alpha));
    int cmax = util::min(int(j + f[i][j][1] + c_ws*alpha), n_cols - patch_size);

    if (rmin > rmax) rmin = rmax = f[i][j][0] + i;
    if (cmin > cmax) cmin = cmax = f[i][j][1] + j;

    /* random_range(min, max) 生成一个[min, max)区间的随机实数 */
    int r_offset = int(util::random_range(rmin, rmax)) - i;
    int c_offset = int(util::random_range(cmin, cmax)) - j;

    Mat patch = pick_patch(img_dst, i, j, 0, 0, patch_size);
    Mat cand = pick_patch(img_ref, i, j, r_offset, c_offset, patch_size);

    float similarity = sim(patch, cand, 1);
    if (similarity > v[i][j]) {
        v[i][j] = similarity;
        f[i][j][0] = r_offset; f[i][j][1] = c_offset;
    }
    alpha *= exp;
}

3. 图像重构

在Barnes et al.的论文中使用了一些比较高级的图像重构(image reshuffling)技术，涉及到更多其他论文的内容，时间限制笔者没有实现，仅采用了“复制-粘贴”的图像重构技术。逻辑比较简单，直接贴代码。

// 图片类型是CV_8UC3，BGR三通道对应三维向量Vec3b，各项取值范围[0, 255]整数
#define CHANNEL_TYPE Vec3b

void reconstruct(Vector3i& f, Mat& dst, Mat& ref, int patch_size) {
    int n_rows = dst.rows, n_cols = dst.cols;
    /* 复制-粘贴，考虑到一个非边界的patch只有左上角点不会被后续patch覆盖 */
    /* 所以对于非边界的像素点，直接取match中的像素点即可 */
    for (int i = 0; i < n_rows; i++) {
        for (int j = 0; j < n_cols; j++) {
            int r = f[i][j][0], c = f[i][j][1];
            dst.at<CHANNEL_TYPE>(i, j) = ref.at<CHANNEL_TYPE>(i + r, j + c);
        }
    }

    /* 处理边界情况 */
    int r = n_rows - patch_size - 1;
    int c;
    for(c = 0; c < n_cols - patch_size - 1; c++) {
        Mat last_patch = pick_patch(ref, r, c, f[r][c][0], f[r][c][1], patch_size);
        for (int i = 0; i < patch_size; i++){
            for (int j = 0; j < patch_size; j++) {
                dst.at<CHANNEL_TYPE>(r + i, c + j) = last_patch.at<CHANNEL_TYPE>(i, j);
            }
        }
    }

    c = n_cols - patch_size - 1;
    for(r = 0; r < n_rows - patch_size - 1; r++) {
        Mat last_patch = pick_patch(ref, r, c, f[r][c][0], f[r][c][1], patch_size);
        for (int i = 0; i < patch_size; i++){
            for (int j = 0; j < patch_size; j++) {
                dst.at<CHANNEL_TYPE>(r + i, c + j) = last_patch.at<CHANNEL_TYPE>(i, j);
            }
        }
    }
}