steger可以提取固定宽度的血管中心线，需要根据不同宽度调整高斯滤波的sigma，
血管越粗sigma越大；血管越细sigma越小。
但是常见的血管图像中存在多条血管，不同血管（或者同一血管的不同位置）的宽度不一样，这会导致steger算法提取的中心线断裂！
高斯滤波（sigma=2.0）结果：
博主因此提出了多尺度自适应高斯滤波。
多尺度自适应高斯滤波的原理：
多尺度是指，sigma的范围是[3.0,0.1]，对应高斯核尺寸范围是[3，29]。
自适应是指，对于图像中的每一点都会向上下左右四个方向生长，测量上下左右四个方向上的w1、w2、w3和w4，最终取最小的宽度作为该点的宽度w0，然后选取k=2*w0+1和对应的sigma的高斯滤波结果作为该点的高斯滤波结果。
k越大，sigma越小，对应的高斯滤波函数就会越突起陡峭，这就可以使得不同位置的高斯滤波结果是自适应多尺度的高斯滤波结果！
多尺度自适应高斯滤波（sigma∈[3.0,0.1]）结果：
第一代的代码：
自适应条件为：
•alpha: 邻域中“不为0元素比例”
• beta: 邻域中的对称系数
代码效率比较低，后面改进了一下自适应条件和代码结构
#include<opencv2/opencv.hpp>
#include<iostream>
using namespace std;
using namespace cv;
int Get_gray(Mat* im, int col, int row)
{
return static_cast<int>(im->at<uchar>(Point(col, row)));
}
//四周细化算法
void Refine(Mat& image) {
int p[8];
int top = 1, down = 1, right = 1, left = 1;
vector<Point> del;
int grayvalue = 0;
int height = image.rows; //获取图像高度
int width = image.cols; //获取图像宽度
Mat *im = reinterpret_cast<Mat*>((void*)&image); //获取像素点信息
uchar* data;
do {
//cout << del.size() << endl;
vector<Point>().swap(del);
//上下收缩
for (int i = 1; i < height - 1; i++)
{
for (int j = 1; j < width - 1; j++)
{
grayvalue = Get_gray(im, j, i); //获取指定点灰度值
if (grayvalue != 0) //判断中心点是否为前景
{
p[0] = (Get_gray(im, j + 1, i) == 0) ? 0 : 1;
p[1] = (Get_gray(im, j + 1, i - 1) == 0) ? 0 : 1;
p[2] = (Get_gray(im, j, i - 1) == 0) ? 0 : 1;
p[3] = (Get_gray(im, j - 1, i - 1) == 0) ? 0 : 1;
p[4] = (Get_gray(im, j - 1, i) == 0) ? 0 : 1;
p[5] = (Get_gray(im, j - 1, i + 1) == 0) ? 0 : 1;
p[6] = (Get_gray(im, j, i + 1) == 0) ? 0 : 1;
p[7] = (Get_gray(im, j + 1, i + 1) == 0) ? 0 : 1;
if (i < height - 2)
down = (Get_gray(im, j, i + 2) == 0) ? 0 : 1;
else
down = 1;
// 横向直线
if (p[6] && (p[5] || p[7] || p[0] || p[4]) && !(p[1] || p[3]) && p[2] == 0 && down)
{
del.push_back(Point(j, i));
}
if (p[2] && (p[1] || p[3] || p[0] || p[4]) && !(p[5] || p[7]) && p[6] == 0)
{
del.push_back(Point(j, i));
}
}
}
}
for (int i = 1; i < height - 2; i++)
{
grayvalue = Get_gray(im, 0, i);
if (grayvalue != 0)
{
if (Get_gray(im, 0, i - 1) && Get_gray(im, 1, i - 1) && Get_gray(im, 0, i + 1) == 0 && Get_gray(im, 1, i) == 0) //上2，上1，右上1，下1=0，右1=0
{
del.push_back(Point(0, i));
}
if (Get_gray(im, 0, i - 1) == 0 && Get_gray(im, 1, i + 1) && Get_gray(im, 1, i) == 0 && Get_gray(im, 0, i + 2))//上1=0，下1，右下1，右1=0，下2
{
del.push_back(Point(0, i));
}
}
if (grayvalue != 0)
{
if (Get_gray(im, width - 1, i - 1) && Get_gray(im, width - 2, i - 1) && Get_gray(im, width - 1, i + 1) == 0 && Get_gray(im, width - 2, i) == 0) //上2，上1，左上1，下1=0，左1=0
{
del.push_back(Point(width - 1, i));
}
if (Get_gray(im, width - 1, i - 1) == 0 && Get_gray(im, width - 2, i + 1) && Get_gray(im, width - 2, i) == 0 && Get_gray(im, width - 1, i + 2))//上1=0，下1，左下1，左1=0，下2
{
del.push_back(Point(width - 1, i));
}
}
}
//for (int i = 0; i < del.size(); i++)
//{
// uchar* data = image.ptr<uchar>(del[i].y);
// data[del[i].x] = 0;
//}
//左右收缩
for (int i = 1; i < height - 1; i++)
{
for (int j = 1; j < width - 1; j++)
{
grayvalue = Get_gray(im, j, i); //获取指定点灰度值
if (grayvalue != 0) //判断中心点是否为前景
{
p[0] = (Get_gray(im, j + 1, i) == 0) ? 0 : 1;
p[1] = (Get_gray(im, j + 1, i - 1) == 0) ? 0 : 1;
p[2] = (Get_gray(im, j, i - 1) == 0) ? 0 : 1;
p[3] = (Get_gray(im, j - 1, i - 1) == 0) ? 0 : 1;
p[4] = (Get_gray(im, j - 1, i) == 0) ? 0 : 1;
p[5] = (Get_gray(im, j - 1, i + 1) == 0) ? 0 : 1;
p[6] = (Get_gray(im, j, i + 1) == 0) ? 0 : 1;
p[7] = (Get_gray(im, j + 1, i + 1) == 0) ? 0 : 1;
if (j < width - 2)
right = (Get_gray(im, j + 2, i) == 0) ? 0 : 1;
else
right = 1;
//竖直线
if (p[0] && (p[1] || p[7] || p[2] || p[6]) && !(p[3] || p[5]) && p[4] == 0 && right)
{
del.push_back(Point(j, i));
}
if (p[4] && (p[3] || p[5] || p[2] || p[6]) && !(p[1] || p[7]) && p[0] == 0)
{
del.push_back(Point(j, i));
}
}
}
}
for (int j = 1; j < width - 2; j++)
{
grayvalue = Get_gray(im, j, 0);
if (grayvalue != 0)
{
if (Get_gray(im, j - 1, 0) == 0 && Get_gray(im, j + 1, 0) && Get_gray(im, j + 2, 0) && Get_gray(im, j, 1) == 0 && Get_gray(im, j + 1, 1)) //左1=0，右1，右2，下1=0，右下1
{
del.push_back(Point(j, 0));
}
if (Get_gray(im, j - 1, 0) && Get_gray(im, j + 1, 0) == 0 && Get_gray(im, j, 1) == 0 && Get_gray(im, j - 1, 1))//左1，右1=0，下1=0，左下1
{
del.push_back(Point(j, 0));
}
}
}
for (int j = 1; j < width - 2; j++)
{
grayvalue = Get_gray(im, j, height - 1);
if (grayvalue != 0)
{
if (Get_gray(im, j - 1, height - 1) == 0 && Get_gray(im, j + 1, height - 1) && Get_gray(im, j + 2, height - 1) && Get_gray(im, j, height - 2) == 0 && Get_gray(im, j + 1, height - 2)) //左1=0，右1，右2，下1=0，右下1
{
del.push_back(Point(j, height - 1));
}
if (Get_gray(im, j - 1, height - 1) && Get_gray(im, j + 1, height - 1) == 0 && Get_gray(im, j, height - 2) == 0 && Get_gray(im, j - 1, height - 2))//左1，右1=0，下1=0，左下1
{
del.push_back(Point(j, height - 1));
}
}
}
for (int i = 0; i < del.size(); i++)
{
data = image.ptr<uchar>(del[i].y);
data[del[i].x] = 0;
}
}while (!del.empty());
}
int main() {
Mat src0 = imread("test2.bmp", 0);
Mat color_src0 = imread("test2.bmp");
//imshow("src", src0);
//cout << "src: " << endl;
//cout << src0 << endl;
//int num = 13;
//double Ks[13] = { 3,5,7,9,11,13,15,17,19,21,23,25,27 };
//double Sigmas[13] = { 2.0,1.8,1.6,1.4,1.2,1.0,0.9,0.8,0.7,0.6,0.5,0.3,0.1 };
int num = 13;
double Ks[13] = { 3,5,7,9,11,13,15,17,19,21,23,25,27 };
double Sigmas[13] = { 3.0,1.9,1.8,1.7,1.6,1.5,1.4,1.2,1.0,0.8,0.6,0.4,0.1 };
Mat float_src;
src0.convertTo(float_src, CV_32FC1);
vector<Mat> gaussianIms(num);
for (int i = 0; i < num; i++) {
GaussianBlur(float_src, gaussianIms[i], Size(Ks[i], Ks[i]), Sigmas[i], Sigmas[i]);
}
//for (int i = 0; i < num; i++) {
// imshow(to_string(Ks[i]) + " Ks", gaussianIms[i]);
// //cout << to_string(Ks[i]) + " Ks: " << endl;
// //cout << gaussianIms[i] << endl;
//}
//Mat res;
//GaussianBlur(float_src, res, Size(21, 21), 2, 2);
Mat res = Mat::zeros(src0.size(), CV_32FC1);
Mat neibor;
int res_rows = res.rows;
int res_cols = res.cols;
float ks, S, alpha=0, beta = 0;
int lupx, lupy, rdpx, rdpy,rect_w,rect_h;
Mat tmp;
ofstream fout("./text.txt");
for (int i = 0; i < res_rows; i++) {
float* res_ptr = res.ptr<float>(i);
for (int j = 0; j < res_cols; j++) {
for (int l = 0; l < num; l++) {
ks = (Ks[l]-1)/2;
lupx = j - ks;
lupy = i - ks;
rdpx = j + ks;
rdpy = i + ks;
lupx = (lupx < 0) ? 0 : lupx;
lupy = (lupy < 0) ? 0 : lupy;
rdpx = (rdpx >= res_cols) ? res_cols-1 : rdpx;
rdpy = (rdpy >= res_rows) ? res_rows-1 : rdpy;
//cout << i << " " << j << " " << k << " " << Point(lupx, lupy)<<" "<< Rect(lupx, lupy, rdpx - lupx + 1, rdpy - lupy + 1) << endl;
rect_w = rdpx - lupx + 1;
rect_h = rdpy - lupy + 1;
neibor = float_src(Rect(lupx, lupy, rect_w, rect_h));
S = rect_w*rect_h;
if (rect_w == rect_h) {//求解对称系数beta，越不对称越大
tmp = neibor.clone();
flip(tmp, tmp, -1);
tmp = abs(tmp - neibor);
beta = cv::countNonZero(tmp);
beta /= S;
}
alpha = float(cv::countNonZero(neibor)) / S;//自适应调节基数alpha
//cout << i << " " << j << " " << k << " " << Point(lupx, lupy)<<" "<< Point(rdpx, rdpy) << endl;
if ((alpha - beta)<0.7) {//(alpha - beta)<0.7,保证血管中心线位置的高斯激活值远大于其他位置
res_ptr[j] = gaussianIms[l].at<float>(i, j);
//fout << Point(j, i) << " " << Rect(lupx, lupy, rdpx - lupx + 1, rdpy - lupy + 1) << " " << ks * 2 + 1 << " " << alpha << " "<< beta << " " << res_ptr[j] << endl;
//resize(neibor, neibor, Size(100, 100));
//imshow("neibor", neibor);
//waitKey(0);
break;
}
else {
if (l == num - 1) {//确保res图像每个都有值
res_ptr[j] = gaussianIms[l].at<float>(i, j);
//fout << Point(j, i) << " " << Rect(lupx, lupy, rdpx - lupx + 1, rdpy - lupy + 1) << " " << ks * 2 + 1 << " " << alpha << " "<< beta << " " << res_ptr[j] << endl;
}
}
}
}
}
fout.close();
//cout << "res: " << endl;
//cout << res << endl;
//system("pause");
GaussianBlur(res, res, Size(5, 5), 2, 2);
Mat img = res.clone();
//一阶偏导数
Mat m1, m2;
m1 = (Mat_<float>(1, 2) << 1, -1); //x偏导
m2 = (Mat_<float>(2, 1) << 1, -1); //y偏导
Mat dx, dy;
filter2D(img, dx, CV_32FC1, m1);
filter2D(img, dy, CV_32FC1, m2);
//二阶偏导数
Mat m3, m4, m5;
m3 = (Mat_<float>(1, 3) << 1, -2, 1); //二阶x偏导
m4 = (Mat_<float>(3, 1) << 1, -2, 1); //二阶y偏导
m5 = (Mat_<float>(2, 2) << 1, -1, -1, 1); //二阶xy偏导
Mat dxx, dyy, dxy;
filter2D(img, dxx, CV_32FC1, m3);
filter2D(img, dyy, CV_32FC1, m4);
filter2D(img, dxy, CV_32FC1, m5);
//filter2D(dxy, dxy, CV_32FC1, m2);
//hessian矩阵
double maxD = -1;
int imgcol = img.cols;
int imgrow = img.rows;
Mat centerline = Mat::zeros(res.size(), CV_8UC1);
//Mat TMap = Mat(res.size(), CV_32FC1);
//Mat tmp_exImg = Mat(srcIm.rows, srcIm.cols, CV_32FC4);//存放对应中心线上点的直径和法线信息
double nx, ny;
double fmaxD = 0;
Mat eValue;
Mat eVectors;
for (int j = 0; j < imgrow; j++) {
const uchar* inData = src0.ptr<uchar>(j);// 每一行图像的指针
//float* inData2 = TMap.ptr<float>(j);// 每一行图像的指针
for (int i = 0; i < imgcol; i++) {
if (inData[i] > 0) {
Mat hessian(2, 2, CV_32FC1);
hessian.at<float>(0, 0) = dxx.at<float>(j, i);
hessian.at<float>(0, 1) = dxy.at<float>(j, i);
hessian.at<float>(1, 0) = dxy.at<float>(j, i);
hessian.at<float>(1, 1) = dyy.at<float>(j, i);
eigen(hessian, eValue, eVectors);
if (fabs(eValue.at<float>(0, 0)) >= fabs(eValue.at<float>(1, 0))) {
nx = eVectors.at<float>(0, 0);
ny = eVectors.at<float>(0, 1);
fmaxD = eValue.at<float>(0, 0);
}
else {
nx = eVectors.at<float>(1, 0);
ny = eVectors.at<float>(1, 1);
fmaxD = eValue.at<float>(1, 0);
}
double t = -(nx*dx.at<float>(j, i) + ny * dy.at<float>(j, i)) / (nx*nx*dxx.at<float>(j, i) + 2 * nx*ny*dxy.at<float>(j, i) + ny * ny*dyy.at<float>(j, i));
//inData2[i] = t;
if (fabs(t*nx) <= 0.5 && fabs(t*ny) <= 0.5) {
centerline.at<uchar>(j, i) = 255;
}
}
}
}
int dilation_size = 1;
Mat element = getStructuringElement(MORPH_ELLIPSE, Size(2 * dilation_size + 1, 2 * dilation_size + 1),Point(dilation_size, dilation_size));
morphologyEx(centerline, centerline, MORPH_CLOSE,element);
Refine(centerline);
for (int j = 0; j < imgrow; j++) {
const uchar* inData = centerline.ptr<uchar>(j);// 每一行图像的指针
for (int i = 0; i < imgcol; i++) {
if (inData[i] > 0) {
color_src0.at<Vec3b>(j, i)[0] = 0;
color_src0.at<Vec3b>(j, i)[1] = 0;
color_src0.at<Vec3b>(j, i)[2] = 255;
}
}
}
imshow("color_src0", color_src0);
imwrite("color_src0.bmp", color_src0);
imshow("centerline", centerline);
//imshow("res", res);
//fout << res << endl;
//fout.close();
//cout << "roi:" << endl;
//cout << res(Rect(168,120,40,32)) << endl;
//cout << res(Rect(112, 171, 40, 35)) << endl;
//cout << TMap(Rect(157, 157, 51, 73)) << endl;
//cout << TMap(Rect(181, 183, 3, 1)) << endl;
//for (int i = 0; i < num; i++) {
// cout << endl << endl;
// cout << "***gaussianIms_" + to_string(Sigmas[i]) + "_" + to_string(Ks[i]) << endl;
// cout << gaussianIms[i](Rect(157, 157, 51, 73)) << endl;
//}
waitKey(0);
return 0;
}
第二代的代码：
自适应条件为：
对于图像中的每一点都会向上下左右四个方向生长，测量上下左右四个方向上的w1、w2、w3和w4，最终取最小的宽度作为该点的宽度w0，然后选取k=2*w0+1和对应的sigma的高斯滤波结果作为该点的高斯滤波结果
#include<opencv2/opencv.hpp>
#include<iostream>
using namespace std;
using namespace cv;
int Get_gray(Mat* im, int col, int row)
{
return static_cast<int>(im->at<uchar>(Point(col, row)));
}
//四周细化算法
void Refine(Mat& image) {
int p[8];
int top = 1, down = 1, right = 1, left = 1;
vector<Point> del;
int grayvalue = 0;
int height = image.rows; //获取图像高度
int width = image.cols; //获取图像宽度
Mat *im = reinterpret_cast<Mat*>((void*)&image); //获取像素点信息
uchar* data;
do {
//cout << del.size() << endl;
vector<Point>().swap(del);
//上下收缩
for (int i = 1; i < height - 1; i++)
{
for (int j = 1; j < width - 1; j++)
{
grayvalue = Get_gray(im, j, i); //获取指定点灰度值
if (grayvalue != 0) //判断中心点是否为前景
{
p[0] = (Get_gray(im, j + 1, i) == 0) ? 0 : 1;
p[1] = (Get_gray(im, j + 1, i - 1) == 0) ? 0 : 1;
p[2] = (Get_gray(im, j, i - 1) == 0) ? 0 : 1;
p[3] = (Get_gray(im, j - 1, i - 1) == 0) ? 0 : 1;
p[4] = (Get_gray(im, j - 1, i) == 0) ? 0 : 1;
p[5] = (Get_gray(im, j - 1, i + 1) == 0) ? 0 : 1;
p[6] = (Get_gray(im, j, i + 1) == 0) ? 0 : 1;
p[7] = (Get_gray(im, j + 1, i + 1) == 0) ? 0 : 1;
if (i < height - 2)
down = (Get_gray(im, j, i + 2) == 0) ? 0 : 1;
else
down = 1;
// 横向直线
if (p[6] && (p[5] || p[7] || p[0] || p[4]) && !(p[1] || p[3]) && p[2] == 0 && down)
{
del.push_back(Point(j, i));
}
if (p[2] && (p[1] || p[3] || p[0] || p[4]) && !(p[5] || p[7]) && p[6] == 0)
{
del.push_back(Point(j, i));
}
}
}
}
for (int i = 1; i < height - 2; i++)
{
grayvalue = Get_gray(im, 0, i);
if (grayvalue != 0)
{
if (Get_gray(im, 0, i - 1) && Get_gray(im, 1, i - 1) && Get_gray(im, 0, i + 1) == 0 && Get_gray(im, 1, i) == 0) //上2，上1，右上1，下1=0，右1=0
{
del.push_back(Point(0, i));
}
if (Get_gray(im, 0, i - 1) == 0 && Get_gray(im, 1, i + 1) && Get_gray(im, 1, i) == 0 && Get_gray(im, 0, i + 2))//上1=0，下1，右下1，右1=0，下2
{
del.push_back(Point(0, i));
}
}
if (grayvalue != 0)
{
if (Get_gray(im, width - 1, i - 1) && Get_gray(im, width - 2, i - 1) && Get_gray(im, width - 1, i + 1) == 0 && Get_gray(im, width - 2, i) == 0) //上2，上1，左上1，下1=0，左1=0
{
del.push_back(Point(width - 1, i));
}
if (Get_gray(im, width - 1, i - 1) == 0 && Get_gray(im, width - 2, i + 1) && Get_gray(im, width - 2, i) == 0 && Get_gray(im, width - 1, i + 2))//上1=0，下1，左下1，左1=0，下2
{
del.push_back(Point(width - 1, i));
}
}
}
//for (int i = 0; i < del.size(); i++)
//{
// uchar* data = image.ptr<uchar>(del[i].y);
// data[del[i].x] = 0;
//}
//左右收缩
for (int i = 1; i < height - 1; i++)
{
for (int j = 1; j < width - 1; j++)
{
grayvalue = Get_gray(im, j, i); //获取指定点灰度值
if (grayvalue != 0) //判断中心点是否为前景
{
p[0] = (Get_gray(im, j + 1, i) == 0) ? 0 : 1;
p[1] = (Get_gray(im, j + 1, i - 1) == 0) ? 0 : 1;
p[2] = (Get_gray(im, j, i - 1) == 0) ? 0 : 1;
p[3] = (Get_gray(im, j - 1, i - 1) == 0) ? 0 : 1;
p[4] = (Get_gray(im, j - 1, i) == 0) ? 0 : 1;
p[5] = (Get_gray(im, j - 1, i + 1) == 0) ? 0 : 1;
p[6] = (Get_gray(im, j, i + 1) == 0) ? 0 : 1;
p[7] = (Get_gray(im, j + 1, i + 1) == 0) ? 0 : 1;
if (j < width - 2)
right = (Get_gray(im, j + 2, i) == 0) ? 0 : 1;
else
right = 1;
//竖直线
if (p[0] && (p[1] || p[7] || p[2] || p[6]) && !(p[3] || p[5]) && p[4] == 0 && right)
{
del.push_back(Point(j, i));
}
if (p[4] && (p[3] || p[5] || p[2] || p[6]) && !(p[1] || p[7]) && p[0] == 0)
{
del.push_back(Point(j, i));
}
}
}
}
for (int j = 1; j < width - 2; j++)
{
grayvalue = Get_gray(im, j, 0);
if (grayvalue != 0)
{
if (Get_gray(im, j - 1, 0) == 0 && Get_gray(im, j + 1, 0) && Get_gray(im, j + 2, 0) && Get_gray(im, j, 1) == 0 && Get_gray(im, j + 1, 1)) //左1=0，右1，右2，下1=0，右下1
{
del.push_back(Point(j, 0));
}
if (Get_gray(im, j - 1, 0) && Get_gray(im, j + 1, 0) == 0 && Get_gray(im, j, 1) == 0 && Get_gray(im, j - 1, 1))//左1，右1=0，下1=0，左下1
{
del.push_back(Point(j, 0));
}
}
}
for (int j = 1; j < width - 2; j++)
{
grayvalue = Get_gray(im, j, height - 1);
if (grayvalue != 0)
{
if (Get_gray(im, j - 1, height - 1) == 0 && Get_gray(im, j + 1, height - 1) && Get_gray(im, j + 2, height - 1) && Get_gray(im, j, height - 2) == 0 && Get_gray(im, j + 1, height - 2)) //左1=0，右1，右2，下1=0，右下1
{
del.push_back(Point(j, height - 1));
}
if (Get_gray(im, j - 1, height - 1) && Get_gray(im, j + 1, height - 1) == 0 && Get_gray(im, j, height - 2) == 0 && Get_gray(im, j - 1, height - 2))//左1，右1=0，下1=0，左下1
{
del.push_back(Point(j, height - 1));
}
}
}
for (int i = 0; i < del.size(); i++)
{
data = image.ptr<uchar>(del[i].y);
data[del[i].x] = 0;
}
} while (!del.empty());
}
int main() {
Mat src0 = imread("test12.bmp", 0);
Mat color_src0 = imread("test12.bmp");
//imshow("src", src0);
//cout << "src: " << endl;
//cout << src0 << endl;
//int num = 13;
//double Ks[13] = { 3,5,7,9,11,13,15,17,19,21,23,25,27 };
//double Sigmas[13] = { 2.0,1.8,1.6,1.4,1.2,1.0,0.9,0.8,0.7,0.6,0.5,0.3,0.1 };
//int num = 13;
//double Ks[13] = { 3,5,7,9,11,13,15,17,19,21,23,25,27 };
//double Sigmas[13] = { 2.0,1.9,1.8,1.7,1.6,1.5,1.4,1.2,1.0,0.8,0.6,0.4,0.1 };
double dur;
clock_t start, end;
start = clock();
int num = 14;
//double Ks[14] = { 3,5,7,9,11,13,15,17,19,21,23,25,27,29 };
double Ks[14] = { 1,2,3,4,5,6,7,8,9,10,11,12,13,14 };
double Sigmas[14] = { 3.0,2.6,2.2,2.0,1.8,1.6,1.4,1.2,1.0,0.8,0.6,0.3,0.2,0.1 };
Mat float_src;
src0.convertTo(float_src, CV_32FC1);
vector<Mat> gaussianIms(num);
for (int i = 0; i < num; i++) {
GaussianBlur(float_src, gaussianIms[i], Size(Ks[i]*2+1, Ks[i]*2+1), Sigmas[i], Sigmas[i]);
//GaussianBlur(float_src, gaussianIms[i], Size(Ks[i], Ks[i]), Sigmas[i], Sigmas[i]);
}
//for (int i = 0; i < num; i++) {
// imshow(to_string(Ks[i]) + " Ks", gaussianIms[i]);
// //cout << to_string(Ks[i]) + " Ks: " << endl;
// //cout << gaussianIms[i] << endl;
//}
//Mat res;
//GaussianBlur(float_src, res, Size(21, 21), 2, 2);
Mat res = Mat::zeros(src0.size(), CV_32FC1);
int res_rows = res.rows;
int res_cols = res.cols;
//int maxl = 0;
int cnt1, cnt2,l,ks;
//vector<Point> pts = { Point(0,-1),Point(0,1),Point(-1,0),Point(1,0),Point(-1,-1),Point(1,1),Point(-1,1),Point(1,-1) };
vector<Point> pts = { Point(1,0),Point(-1,0),Point(0,1),Point(0,-1)};
Point p,pk;
float* res_ptr = nullptr;
for (int i = 0; i < res_rows; i++) {
res_ptr = res.ptr<float>(i);
for (int j = 0; j < res_cols; j++) {
/*
cnt2 = 1e10;
cnt1 = 0;
k = j+1;
while (k < res_cols&&src0.ptr<uchar>(i)[k] == 255) {
cnt1++;
k++;
}
if (cnt1 < cnt2) {
cnt2 = cnt1;
}
cnt1 = 0;
k = j - 1;
while (k >=0&& src0.ptr<uchar>(i)[k] == 255) {
cnt1++;
k--;
}
if (cnt1 < cnt2) {
cnt2 = cnt1;
}
cnt1 = 0;
k = i + 1;
while (k < res_rows&&src0.ptr<uchar>(k)[j] == 255) {
cnt1++;
k++;
}
if (cnt1 < cnt2) {
cnt2 = cnt1;
}
cnt1 = 0;
k = i - 1;
while (k >= 0 && src0.ptr<uchar>(k)[j] == 255) {
cnt1++;
k--;
}
if (cnt1 < cnt2) {
cnt2 = cnt1;
}
*/
cnt2 = 1e10;
p = Point(j, i);
for (Point& pt : pts) {
cnt1 = 0;
pk = p + (cnt1 + 1)*pt;
while (pk.x >= 0 && pk.x < res_cols&&pk.y >= 0 && pk.y < res_rows&&src0.at<uchar>(pk) >0 && cnt1<cnt2) {
cnt1++;
pk = p + (cnt1 + 1)*pt;
}
if (cnt1 < cnt2) {
cnt2 = cnt1;
}
}
l = 0;
for (; l < num; l++) {
//ks = (Ks[l] - 1) * 2;
if (Ks[l] > cnt2) {
break;
}
}
//if (l > maxl) {
// maxl = l;
//}
res_ptr[j] = gaussianIms[l].at<float>(i, j);
//fout << Point(j, i) << " " << ks * 2 + 1 << " " << res_ptr[j] << endl;
}
}
//fout.close();
//cout << "max k: " << Ks[maxl]*2+1 << endl;
//cout << "res: " << endl;
//cout << res << endl;
//system("pause");
//normalize(res, res, 1, 0, cv::NORM_MINMAX);
GaussianBlur(res, res, Size(3, 3), 1, 1);
//GaussianBlur(res, res, Size(5, 5), 1, 1);
Mat img = res.clone();
//一阶偏导数
Mat m1, m2;
m1 = (Mat_<float>(1, 2) << 1, -1); //x偏导
m2 = (Mat_<float>(2, 1) << 1, -1); //y偏导
Mat dx, dy;
filter2D(img, dx, CV_32FC1, m1);
filter2D(img, dy, CV_32FC1, m2);
//二阶偏导数
Mat m3, m4, m5;
m3 = (Mat_<float>(1, 3) << 1, -2, 1); //二阶x偏导
m4 = (Mat_<float>(3, 1) << 1, -2, 1); //二阶y偏导
m5 = (Mat_<float>(2, 2) << 1, -1, -1, 1); //二阶xy偏导
Mat dxx, dyy, dxy;
filter2D(img, dxx, CV_32FC1, m3);
filter2D(img, dyy, CV_32FC1, m4);
filter2D(img, dxy, CV_32FC1, m5);
//filter2D(dxy, dxy, CV_32FC1, m2);
//hessian矩阵
double maxD = -1;
int imgcol = img.cols;
int imgrow = img.rows;
Mat centerline = Mat::zeros(res.size(), CV_8UC1);
//Mat TMap = Mat(res.size(), CV_32FC1);
//Mat tmp_exImg = Mat(srcIm.rows, srcIm.cols, CV_32FC4);//存放对应中心线上点的直径和法线信息
double nx, ny;
double fmaxD = 0;
Mat eValue;
Mat eVectors;
Mat hessian(2, 2, CV_32FC1);
double t;
for (int j = 0; j < imgrow; j++) {
const uchar* inData = src0.ptr<uchar>(j);// 每一行图像的指针
//float* inData2 = TMap.ptr<float>(j);// 每一行图像的指针
for (int i = 0; i < imgcol; i++) {
if (inData[i] > 0) {
//Mat hessian(2, 2, CV_32FC1);
hessian.at<float>(0, 0) = dxx.at<float>(j, i);
hessian.at<float>(0, 1) = dxy.at<float>(j, i);
hessian.at<float>(1, 0) = dxy.at<float>(j, i);
hessian.at<float>(1, 1) = dyy.at<float>(j, i);
eigen(hessian, eValue, eVectors);
if (fabs(eValue.at<float>(0, 0)) >= fabs(eValue.at<float>(1, 0))) {
nx = eVectors.at<float>(0, 0);
ny = eVectors.at<float>(0, 1);
fmaxD = eValue.at<float>(0, 0);
}
else {
nx = eVectors.at<float>(1, 0);
ny = eVectors.at<float>(1, 1);
fmaxD = eValue.at<float>(1, 0);
}
t = -(nx*dx.at<float>(j, i) + ny * dy.at<float>(j, i)) / (nx*nx*dxx.at<float>(j, i) + 2 * nx*ny*dxy.at<float>(j, i) + ny * ny*dyy.at<float>(j, i));
//inData2[i] = t;
if (fabs(t*nx) <= 0.7 && fabs(t*ny) <= 0.7){
centerline.at<uchar>(j, i) = 255;
}
}
}
}
int dilation_size = 1;
Mat element = getStructuringElement(MORPH_ELLIPSE, Size(2 * dilation_size + 1, 2 * dilation_size + 1), Point(dilation_size, dilation_size));
morphologyEx(centerline, centerline, MORPH_CLOSE, element);
Refine(centerline);
end = clock();
dur = (double)(end - start);
printf("Use Time:%f\n", (dur / CLOCKS_PER_SEC));
for (int j = 0; j < imgrow; j++) {
const uchar* inData = centerline.ptr<uchar>(j);// 每一行图像的指针
for (int i = 0; i < imgcol; i++) {
if (inData[i] > 0) {
color_src0.at<Vec3b>(j, i)[0] = 0;
color_src0.at<Vec3b>(j, i)[1] = 0;
color_src0.at<Vec3b>(j, i)[2] = 255;
}
}
}
imshow("color_src0", color_src0);
imwrite("color_src0.bmp", color_src0);
imshow("centerline", centerline);
//imshow("res", res);
//fout << res << endl;
//fout.close();
//cout << "roi:" << endl;
//cout << res(Rect(169,175,26,17)) << endl;
//cout << res(Rect(176, 130, 26, 22)) << endl;
//cout << res(Rect(151, 191, 26, 20)) << endl;
//cout << res(Rect(112, 171, 40, 35)) << endl;
//cout << TMap(Rect(176, 130, 26, 22)) << endl;
//cout << TMap(Rect(181, 183, 3, 1)) << endl;
//for (int i = 0; i < num; i++) {
// cout << endl << endl;
// cout << "***gaussianIms_" + to_string(Sigmas[i]) + "_" + to_string(Ks[i]) << endl;
// cout << gaussianIms[i](Rect(157, 157, 51, 73)) << endl;
//}
waitKey(0);
return 0;
}
1515
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
