距离变换及其数学推导
来源:互联网 发布:交互原型设计软件 编辑:程序博客网 时间:2024/06/05 08:47
距离变换
假设二值图像I,其中有前景点集F,背景点集B。则由通过距离变换可以得到距离图D。
其中
由这个简单的想法我们可以得到距离变换的实现代码
double ComputeEuDistance(Point2i x1, Point2i x2){ return norm(x1 - x2);}Mat BruteForceDistanceTransform(Mat binary_img){ vector<Point2i> foreth; vector<Point2i> background; for (size_t i = 0; i < binary_img.rows; i++) { for (size_t j = 0; j < binary_img.cols; j++) { if (binary_img.at<uchar>(j, i)) foreth.emplace_back(j, i); else background.emplace_back(j, i); } } Mat result = Mat::zeros(binary_img.rows, binary_img.cols, CV_8U); for (auto& pf : foreth) { double mind = INFINITY; for (auto& bf : background) { auto d = ComputeEuDistance(pf, bf); mind = min(mind, d); } result.at<uchar>(pf) = mind; } double minimal, maximum; minMaxLoc(result, &minimal, &maximum); result = (result - minimal) * (255.0 / (maximum - minimal)); return result;}
L1-norm和L∞-norm下的加速
显然,上述算法时间复杂度过高,消耗的时间往往是不能接受的。
实际实现的时候,我们可以利用图像在空间上的连续性和在存储上的离散性进行加速。
在1-范数和∞-范数下(也就是所谓的街区距离和棋盘距离),我们可以通过相邻像素点的空间关系得到线性时间复杂度的算法。
1-范数下,我们可以通过对图像进行一次从左上到右下的遍历和一次相反方向的遍历来计算距离变换。第一次遍历的时候对当前前景点p,和它左边与上边的点
第二次遍历的时候对当前前景点p,和它右边与下边的点
由此得到如下代码:
Mat L1DistanceTransform(Mat binary_img){ Mat result = binary_img.clone()/255; for (size_t i = 1; i < result.rows; i++) { for (size_t j = 1; j < result.cols; j++) { if (result.at<uchar>(j, i)) { result.at<uchar>(j, i) = min(result.at<uchar>(j - 1, i), result.at<uchar>(j, i - 1)) + 1; } } } for (int i = result.rows - 2; i >=0; i--) { for (int j = result.cols - 2; j >=0; j--) { if (result.at<uchar>(j, i)) { result.at<uchar>(j, i) = min((int)result.at<uchar>(j, i), (min(result.at<uchar>(j + 1, i), result.at<uchar>(j, i + 1))+1) ); } } } double minimal, maximum; minMaxLoc(result, &minimal, &maximum); result = (result - minimal) * (255.0 / (maximum - minimal)); return result;}
对于
L2-norm下的加速
使用棋盘距离和街区距离可以得到线性复杂度的算法,然而问题是这两种距离都没有旋转不变性,也不符合人类直觉。在许多场合欧氏距离进行变换可以得到更好的效果。
好在使用类似的思路,我们也可以对欧氏距离下的距离变换进行加速。假设前景一点p和它周围的一点q,假设已经算出了q的距离变换且离q最近的点为r,d(p,r)>d(q,r),我们有
于是有
其中
同样的,通过两次遍历我们可以得到图像的距离变换,代码如下:
inline int hqp(const int Rx, const int Ry, const int i){ static int G[][2] = { { 1,0 },{ 1,1 },{ 0,1 },{ 1,1 },{ 1,0 },{ 1,1 },{ 0,1 },{ 1,1 } }; return G[i][0] * (2 * Rx + 1) + G[i][1] * (2 * Ry + 1);}//EDTMat EuclidienDistanceTransform(const Mat& binary_img){ /* q2 q3 q4 q1 p q5 q8 q7 q6 */ pair<Mat, Mat> R{ Mat::zeros(binary_img.rows, binary_img.cols, CV_32S) , Mat::zeros(binary_img.rows, binary_img.cols, CV_32S) }; Mat result = Mat::zeros(binary_img.rows, binary_img.cols, CV_32S); static int G[][2] = { { 1,0 },{ 1,1 },{ 0,1 },{ 1,1 },{ 1,0 },{ 1,1 },{ 0,1 },{ 1,1 } }; static int q[8][2] = { { -1,0 },{ -1, -1 },{ 0,-1 },{ 1,-1 },{ 1,0 },{ 1,1 },{ 0,1 },{ -1,1 } }; int qindex[8]; const int h = result.rows; const int w = result.cols; for (size_t i = 0; i < 8; i++) { qindex[i] = q[i][0] + q[i][1] * w; } int* result_ptr = result.ptr<int>(); int* Rx_ptr = R.first.ptr<int>(); int* Ry_ptr = R.second.ptr<int>(); const unsigned char* binary_img_ptr = binary_img.ptr<uchar>(); for (size_t i = 1; i < h - 1; i++) { int index = w*i; for (size_t j = 1; j < w - 1; j++) { index++; if (binary_img_ptr[index]) { int& fp = result_ptr[index]; fp = INT16_MAX; int& Rpx = Rx_ptr[index]; int& Rpy = Ry_ptr[index]; for (size_t k = 0; k < 4; k++) { int absqindex = index + qindex[k]; auto Rx = Rx_ptr[absqindex]; auto Ry = Ry_ptr[absqindex]; auto hp = result_ptr[absqindex] + hqp(Rx, Ry, k); if (hp< fp) { fp = hp; Rpx = Rx + G[k][0]; Rpy = Ry + G[k][1]; } } } } } for (int i = h - 2; i >= 1; i--) { int index = w*i + w - 1; for (int j = w - 2; j >= 1; j--) { index--; if (binary_img_ptr[index]) { int& fp = result_ptr[index]; int& Rpx = Rx_ptr[index]; int& Rpy = Ry_ptr[index]; for (int k = 4; k < 8; k++) { int absqindex = index + qindex[k]; auto Rx = Rx_ptr[absqindex]; auto Ry = Ry_ptr[absqindex]; auto hp = result_ptr[absqindex] + hqp(Rx, Ry, k); if (hp < fp) { fp = hp; Rpx = Rx + G[k][0]; Rpy = Ry + G[k][1]; } } } } } Mat result1; result.convertTo(result1, CV_32F); sqrt(result1, result1); double minimal, maximum; minMaxLoc(result1, &minimal, &maximum); //cout << minimal << " " << maximum << endl; result1 = (result1 - minimal) * (255.0 / (maximum - minimal)); Mat rr; result1.convertTo(rr, CV_8U); return rr;}
阅读全文
0 0
- 距离变换及其数学推导
- 逻辑回归及其数学推导
- 朴素贝叶斯及其数学推导
- OpenGL ES 矩阵变换及其数学原理
- 矩阵变换:沿任意轴旋转及其推导
- 矩阵变换:沿任意轴旋转及其推导
- 矩阵变换:沿任意轴旋转及其推导
- 阵变换:沿任意轴旋转及其推导
- 矩阵变换:沿任意轴旋转及其推导
- 距离变换
- 距离变换
- 距离变换
- 距离变换
- 距离变换
- 距离变换
- 距离变换
- 距离变换
- 距离变换
- http 505 The server does not support the HTTP protocol version used in the request
- 2017年10月26日训练总结
- MyBatis 3(中文版) 第四章 使用注解配置SQL映射器
- Visual Studio 2017 添加引用报错
- 解决Hadoop启动时,没有启动datanode
- 距离变换及其数学推导
- 读懂正则表达式就这么简单
- 第二次深度实习面试(七鑫易维)
- 剑指offer---数字在排序数组中出现的次数
- Linux学习笔记2——VM虚拟机网络环境设置
- 关于微机原理的一些个人见解(2)
- 面试汇总
- python json 爬京东商品评论
- c语言———预处理,结构体