李纯明水平集源码C++复现笔记
一、复现环境
vs工程,windows系统,语言使用C++,算法库使用opencv3.1.0,opencv有一些很好用的矩阵操作函数,防止重复造轮子。
本篇文章不介绍原理推敲,只介绍C++复现中的几个易错点,主要是纽曼边界条件和拉普拉斯算子的重写,并在文章中附上这两个函数的代码。
大佬的源码(Matlab)和原文章很好找,不再做冗余的复制粘贴。
二、复现要点
对于里面的函数,基本的矩阵点乘点除等,直接使用opencv库没有任何问题,需要重写的有:
(1)矩阵正余弦函数;
(2)像素点比较,不要用opencv自带的比较算子,自己写一个很简单,opencv的函数比不出正确结果;
(3)纽曼边界条件需要重写;
(4)梯度和二阶梯度函数需要重写,sobel算子只是个卷积操作,系数不对;
(5)拉普拉斯算子需要重写,Laplacian算子也只是个卷积,系数不对,如果直接使用opencv的算子,最终结果不会收敛,像素数值基本到第十几轮迭代就会达到几十万,再往后就会超过double的精度。
对于纽曼边界条件和拉普拉斯算子复现对程序正确运行影响较大,边界点处理有一定难度,本文主要讲述这两个函数的工程复现。
1、纽曼边界条件(NeumannBoundCond)复现
这个函数保证边界的微分方程有解。
根据Matlab的源码理解,首先求四个顶点,再求四条边。
(1)求取顶点
在图中,用红色的点替代绿色点。
// 修改矩阵四个角的坐标
f.at<double>(0, 0) = g.at<double>(2, 2);
f.at<double>(0, ncol - 1) = g.at<double>(2, ncol - 3);
f.at<double>(nrow - 1, 0) = g.at<double>(nrow - 3, 2);
f.at<double>(nrow - 1, ncol - 1) = g.at<double>(nrow - 3, ncol - 3);
(2)求取上下边界
在图中,用红色的区域代替上边界,下边界同理类推。
// 修改上下边缘的像素值
for (int j = 1; j < ncol - 1; j++) {
f.at<double>(0, j) = g.at<double>(2, j);
f.at<double>(nrow - 1, j) = g.at<double>(nrow - 3, j);
}
(3)左右边界
在图中,用红色的区域代替左边界,右边界同理类推。
// 修改左右边缘的像素值
for (int i = 1; i < nrow - 1; i++) {
f.at<double>(i, 0) = g.at<