局部极值提取算法
这是http://www.imagepy.org/的作者原创,我只是对其理解之后改进和说明,欢迎大家使用这个小软件!
有需要C++的朋友,可有留邮箱,经测试和改进的版本!
算法原理:(按照程序思路来,当然其中很多可以改进的地方)

第一步:
先进行距离变换(这类问题都是要通过几何操作进行,距离变换就是几何距离)

第二步:
先找全局可能最大值,然后对图像进行标记
全局可能最大值是基础,3X3核在图像上滑动去找全局可能最大值。标记是为了掩膜操作更方便,你看opencv很多函数都有掩膜,PS上面也有那个白色的膜。
图片结果是:背景和最大值标记3,前景标记1,图像边界标记0(为了防止越界)
貌似下面的图片显示有问题,matploylib的显示有时候确实有问题,Debug的时候就和显示不一样,没关系这一步比较简单。

第三步:
对全局可能最大值进行排除,这一步是程序核心!这里以求最大值为例,最小值相反。
-
-
-
-
- 首先对找到的最大值进行排序--从大到小
- 按照从大到小的顺序去寻找每个最大值的领地
- 领地按照“目标范围”和“目标深度”两个原则,具体表现形式为:“不重合”、“不包括”、“领主大于等于领地”
-
-
-
下面以图说明两个大原则:
目标范围:
下图A1这个最大最他所占领的范围由用户输入,假设半径为y,则 x1 = x2 = y ,x1和x2的范围就是A1的范围,P点虽然很低,但是不属于A1的范围,所以不去占领。

目标深度:
下图A1的深度也又用户输入,假设输入为y,则 y1 = y,高度确定之后,那么P点虽然也在山脊上,但是不会被A1所占领。

下面以图为例进一步说明三个具体形式:
注释:三个具体形式是建立在两个原则之上的,所以理解上面是结局算法的基础。
不重合形式:
下图A1和A2有领土纠纷,但是A1峰值高所以先进行领地划分,当A2进行领地划分的时候一旦接触A1的领地就立马被消灭,最后只剩A1的存在。

不包括形式:
下图A1先进行领土划分,过程中直接把A2消灭了,不存在A2的领土划分

领主大于等于领地形式,也可以叫同等峰值先来先得形式:
由于A1先进行领土划分,所以后面的都被消灭了,读到这里你会发现这个程序的BUG,后面再详细谈论这一点

目标深度原则:
要确定和你预定的深度是不是符合,比如你想要山峰之间的深度都是在10cm以上,那么有些不符合的就得去除,去除原则肯定矮的GG

注释:大佬的程序没有几何距离,还有部分Bug,部分经过修正了。

程序有待改进:
这部分可以通过处理---“同样的峰值且山脉相连”的一类极值点,对其平均、覆盖范围最多等操作取一个中间的值!
由于python用的不熟练,而且现在时间太晚了脑子转不动,以后用C++写的时候再改进这一点。




2017.10.17更新:

1.对上面第二步:可能最大值进行补充
请看下面这个图,中间位置肯定不是最大值区域,而且通过上面的第三步过滤也是没办法过滤的,具体看上面原理就知道了,所以需要在第二步进行过滤!
本程序利用如下方法:
A.每一层的最大值都大于外面一层,sum[0]>=sum[1] && sum[1]>sum[2] && sum[2]>sum[3],这里使用到9X9的内核
B.记录每一层最大值出现的位置X0和X1,然后最大值的距离(abs(x0.x - x1.x) <= 2 && abs(x0.y - x1.y) <= 2)小于等于2,这里不严谨,因为每层最大值不止一个。
C.记录最大值出现的数量n_count >= 3
1 float sum[4] = { 0 };//存储最大值进行对比、在这里说不清楚看博客2017.10.17更新
2 sum[0] = img[3][j];
3 Point x0 = Point(0, 0);
4 Point x1 = Point(0, 0);
5 uchar n_count = 0;
6 for (int m = 2; m < 5; m++)
7 {
8 for (int n = -1; n < 2; n++)
9 {
10 if (m == 3 && n == 0) continue;
11 sum[1] = sum[1] < img[m][j + n] ? img[m][j + n] : sum[1];
12 x0 = sum[0] == img[m][j + n] ? Point(m, n) : x0;
13 n_count = sum[0] == img[m][j + n] ? n_count+1 : n_count;
14 //flag = img[3][j + 0] < img[m][j + n] ? true : flag;//如果目标像素小于周围任何一个像素就说明这个一定不是最大值
15 }
16 }
17 for (int m = 1; m < 6; m++)
18 {
19 for (int n = -2; n < 3; n++)
20 {
21 if (2 <= m && m <= 4 && -1 <= n && n <= 1) continue;
22 sum[2] = sum[2] < img[m][j + n] ? img[m][j + n] : sum[2];
23 x1 = sum[0] == img[m][j + n] ? Point(m, n) : x1;
24 n_count = sum[0] == img[m][j + n] ? n_count+1 : n_count;
25 //flag = img[3][j + 0] < img[m][j + n] ? true : flag;//如果目标像素小于周围任何一个像素就说明这个一定不是最大值
26 }
27 }
28 for (int m = 0; m < 7; m++)
29 {
30 for (int n = -3; n < 4; n++)
31 {
32 sum[3] = sum[3] < img[m][j + n] && !(1 <= m && m <= 5 && -2 <=n && n <= 2) ? img[m][j + n] : sum[3];
33 //flag = img[3][j+0] < img[m][j + n] ? true : flag;//如果目标像素小于周围任何一个像素就说明这个一定不是最大值
34 }
35 }
36 x1 = (x1.x == 0 && x1.y == 0) || n_count >= 3 ? x0 : x1;//判断是否存在5X5的最大值(和目标点相同)
37 int tmp = sum[0] >= sum[1] && sum[1] >= sum[2] && sum[2] >= sum[3] && (abs(x0.x - x1.x) <= 2 && abs(x0.y - x1.y) <= 2)
38 ? 0 : FillBlock(src, mask_tmp, Point(j, i));//tmp没意义,就是为了调用后面的函数而已
39 }


2.对第三步进行补充
A.搜索过程遇到边界,那就把这个最大值点去除if (img[i][fill_point[count].x + j] == 2 || msk[i][fill_point[count].x + j] == 0) max_point[k].data = 1;
3.效果图
注释:满足一个条件就判定为最大值!
1 Find_Max(img, mask,0, 5000);

1 Find_Max(img, mask,5000, 0);

1 Find_Max(img, mask,5000, 5000);

1 Find_Max(img, mask,10, 20);

1 Find_Max(img, mask,10, 50);

1 import scipy.ndimage as ndimg
2 import numpy as np
3 from numba import jit
4 import cv2
5
6 def neighbors(shape):
7 dim = len(shape)
8 block = np.ones([3] * dim)
9 block[tuple([1] * dim)] = 0
10 idx = np.where(block > 0)
11 idx = np.array(idx, dtype=np.uint8).T
12 idx = np.array(idx - [1] * dim)
13 acc = np.cumprod((1,) + shape[::-1][:-1])
14 return np.dot(idx, acc[::-1])
15
16
17 @jit # trans index to r, c...
18
19 def idx2rc(idx, acc):
20 rst = np.zeros((len(idx), len(acc)), dtype=np.int16)
21 for i in range(len(idx)):
22 for j in range(len(acc)):
23 rst[i, j] = idx[i] // acc[j]
24 idx[i] -= rst[i, j] * acc[j]
25 return rst
26
27
28 #@jit # fill a node (may be two or more points)
29
30 def fill(img, msk, p, nbs, buf):
31 msk[p] = 3
32 buf[0] = p
33 back = img[p]
34 cur = 0
35 s = 1
36 while cur < s:
37 p = buf[cur]
38 for dp in nbs:
39 cp = p + dp
40 if img[cp] == back and msk[cp] == 1:
41 msk[cp] = 3
42 buf[s] = cp
43 s += 1
44 if s == len(buf):
45 buf[:s - cur] = buf[cur:]
46 s -= cur
47 cur = 0
48 cur += 1
49 #msk[p] = 3
50
51
52 #@jit # my mark
53
54 def mark(img, msk, buf, mode): # mark the array use (0, 1, 2)
55 omark = msk
56 nbs = neighbors(img.shape)
57 idx = np.zeros(1024 * 128, dtype=np.int64)
58 img = img.ravel() # 降维
59 msk = msk.ravel() # 降维
60 s = 0
61 for p in range(len(img)):
62 if msk[p] != 1: continue
63 flag = False
64 for dp in nbs:
65 if mode and img[p + dp] > img[p]:
66 flag = True
67 break
68 elif not mode and img[p + dp] < img[p]:
69 flag = True
70 break
71
72 if flag : continue
73 else : fill(img, msk, p, nbs, buf)
74 idx[s] = p
75 s += 1
76 if s == len(idx): break
77 plt.imshow(omark, cmap='gray')
78 return idx[:s].copy()
79
80
81
82 def filter(img, msk, idx, bur, tor, mode):
83 omark = msk
84 nbs = neighbors(img.shape)
85 acc = np.cumprod((1,) + img.shape[::-1][:-1])[::-1]
86 img = img.ravel()
87 msk = msk.ravel()
88
89 arg = np.argsort(img[idx])[::-1 if mode else 1]
90
91 for i in arg:
92 if msk[idx[i]] != 3:
93 idx[i] = 0
94 continue
95 cur = 0
96 s = 1
97 bur[0] = idx[i]
98 while cur < s:
99 p = bur[cur]
100 if msk[p] == 2:
101 idx[i] = 0
102 break
103
104 for dp in nbs:
105 cp = p + dp
106 if msk[cp] == 0 or cp == idx[i] or msk[cp] == 4: continue
107 if mode and img[cp] < img[idx[i]] - tor: continue
108 if not mode and img[cp] > img[idx[i]] + tor: continue
109 bur[s] = cp
110 s += 1
111 if s == 1024 * 128:
112 cut = cur // 2
113 msk[bur[:cut]] = 2
114 bur[:s - cut] = bur[cut:]
115 cur -= cut
116 s -= cut
117
118 if msk[cp] != 2: msk[cp] = 4
119 cur += 1
120 msk[bur[:s]] = 2
121 #plt.imshow(omark, cmap='gray')
122
123 return idx2rc(idx[idx > 0], acc)
124
125
126 def find_maximum(img, tor, mode=True):
127 msk = np.zeros_like(img, dtype=np.uint8)
128 msk[tuple([slice(1, -1)] * img.ndim)] = 1
129 buf = np.zeros(1024 * 128, dtype=np.int64)
130 omark = msk
131 idx = mark(img, msk, buf, mode)
132 plt.imshow(msk, cmap='gray')
133 idx = filter(img, msk, idx, buf, tor, mode)
134 return idx
135
136
137 if __name__ == '__main__':
138 from scipy.misc import imread
139 from scipy.ndimage import gaussian_filter
140 from time import time
141 import matplotlib.pyplot as plt
142
143 img = cv2.imread('123.png')
144 img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
145 ret2, img = cv2.threshold(img, 0, 255, cv2.THRESH_BINARY | cv2.THRESH_OTSU)
146 img[:] = ndimg.distance_transform_edt(img)
147 plt.imshow(img, cmap='gray')
148 pts = find_maximum(img, 20, True)
149 start = time()
150 pts = find_maximum(img, 10, True)
151 print(time() - start)
152 plt.imshow(img, cmap='gray')
153 plt.plot(pts[:, 1], pts[:, 0], 'y.')
154 plt.show()
C++版本
老版本、不稳定,可以看思路
1 //---fill black value
2 int FillBlock(Mat src, Mat &mask, Point center)
3 {
4 uchar back = src.at<uchar>(center.y, center.x);
5 vector<Point> fill_point;
6 int count = 0, count_mount = 1;
7 fill_point.push_back(center);
8 while (count < count_mount)
9 {
10 vector<uchar*> img;
11 vector<uchar*> msk;
12 for (int i = -1; i < 2; i++)
13 {
14 img.push_back(src.ptr<uchar>(fill_point[count].y + i));
15 msk.push_back(mask.ptr<uchar>(fill_point[count].y + i));
16 }
17 for (size_t i = 0; i < 3; i++)
18 {
19 for (int j = -1; j < 2; j++)
20 {
21 if (img[i][fill_point[count].x + j] == back && !(j == 0 && i == 1) && msk[i][fill_point[count].x + j] == 255)
22 {
23 fill_point.push_back(Point(fill_point[count].x + j, fill_point[count].y + i - 1));
24 msk[i][fill_point[count].x + j] = 1;
25 }
26 }
27 }
28 msk[1][fill_point[count].x] = 1;
29 count_mount = fill_point.size() - 1;
30 fill_point.erase(fill_point.begin());
31 }
32 return 0;
33 }
34 //---cal mask
35 //---@_src
36 //---@mask
37 void MaskImage(InputArray _src, Mat &mask)
38 {
39 Mat src = _src.getMat(),mask_tmp = Mat::zeros(src.size(), CV_8UC1);
40 mask_tmp.setTo(255);
41 Mat rows = Mat::zeros(Size(src.cols, 1), CV_8UC1), cols = Mat::zeros(Size(1, src.rows), CV_8UC1);
42 Mat src_rows_beg = mask_tmp.row(0);
43 Mat src_rows_end = mask_tmp.row(src.rows - 1);
44 Mat src_cols_beg = mask_tmp.col(0);
45 Mat src_cols_end = mask_tmp.col(src.cols - 1);
46 rows.copyTo(src_rows_beg); rows.copyTo(src_rows_end);
47 cols.copyTo(src_cols_beg); cols.copyTo(src_cols_end);
48 for (size_t i = 1; i < src.rows-1; i++)
49 {
50 uchar *img0 = src.ptr<uchar>(i - 1);
51 uchar *img = src.ptr<uchar>(i);
52 uchar *img1 = src.ptr<uchar>(i + 1);
53 uchar *msk = mask_tmp.ptr<uchar>(i);
54 for (size_t j = 1; j < src.cols-1; j++)
55 {
56 bool flag = false;
57 //msk[j] = img[j] == 0 ? 0 : msk[j];
58 if (msk[j] != 255) continue;
59 flag = (img[j] < img[j - 1] || img[j] < img[j + 1]
60 || img[j] < img0[j] || img[j] < img0[j - 1]
61 || img[j] < img0[j + 1] || img[j] < img1[j]
62 || img[j] < img1[j - 1] || img[j] < img1[j + 1])
63 ? true : false;
64 int tmp = flag == true ? FillBlock(src, mask_tmp, Point(j, i)) : 0;
65 }
66 }
67 mask = mask_tmp.clone();
68 }
69 //---filter parts max value
70 //---@
71 //---@
72 //---@gap
73 //---@radius
74
75 vector<Point> Find_Max(InputArray _src, Mat&mask,int gap,int radius)
76 {
77 Mat src = _src.getMat();
78
79 typedef struct MyStruct
80 {
81 Point position;
82 float data;
83 }MyStruct;
84
85 MaskImage(src, mask);
86 vector<MyStruct> max_point;
87 for (size_t i = 0; i < src.rows; i++)
88 {
89 uchar *img = src.ptr<uchar>(i);
90 uchar *msk = mask.ptr<uchar>(i);
91 for (size_t j = 0; j < src.cols; j++)
92 {
93 if (msk[j] != 255) continue;
94 MyStruct my_struct;
95 my_struct.data = img[j];
96 my_struct.position = Point(j, i);
97 max_point.push_back(my_struct);
98 }
99 }
100 for (size_t i = 0; i < max_point.size(); i++)
101 {
102 for (size_t j = i; j < max_point.size(); j++)
103 {
104 MyStruct temp;
105 if (max_point[i].data <= max_point[j].data)
106 {
107 if (max_point[j].data == 0) continue;
108 temp = max_point[j];
109 max_point[j] = max_point[i];
110 max_point[i] = temp;
111 }
112 }
113 }
114 //---find max
115 for (size_t k = 0; k < max_point.size(); k++)//---
116 {
117 uchar back = src.at<uchar>(max_point[k].position.y, max_point[k].position.x);
118 vector<Point> fill_point;
119 int count = 0, count_mount = 1;
120 fill_point.push_back(max_point[k].position);
121
122 while (count < count_mount && max_point[k].data != 1)
123 {
124 vector<uchar*> img;
125 vector<uchar*> msk;
126 for (int i = -1; i < 2; i++)
127 {
128 img.push_back(src.ptr<uchar>(fill_point[count].y + i));
129 msk.push_back(mask.ptr<uchar>(fill_point[count].y + i));
130 }
131 for (int i = 0; i < 3; i++)
132 {
133 for (int j = -1; j < 2; j++)
134 {
135 //---
136 uchar x = pow((max_point[k].position.x - fill_point[count].x + j), 2); //(max_point[k].position.x - img[i][fill_point[count].x + j])*(max_point[k].position.x - img[i][fill_point[count].x + j]);
137 uchar y = pow((max_point[k].position.y - (fill_point[count].y + i - 1)) , 2); // (max_point[k].position.y - img[i][fill_point[count].y + j])*(max_point[k].position.y - img[i][fill_point[count].x + j]);
138 uchar distance = sqrt(x + y);
139 if (img[i][fill_point[count].x + j] <= img[1][fill_point[count].x] - gap
140 || msk[i][fill_point[count].x + j] == 3
141 || msk[i][fill_point[count].x + j] == 0
142 || (j == 0 && i == 1)
143 || distance >= radius) continue;
144 if (img[i][fill_point[count].x + j] == 2) max_point[k].data = 1;
145 msk[i][fill_point[count].x + j] = 3;
146 fill_point.push_back(Point(fill_point[count].x + j, fill_point[count].y + i - 1));
147 count_mount++;
148 }
149 }
150 count++;
151 }
152 if (max_point[k].data == 1)
153 {
154 for (size_t i = 0; i < fill_point.size(); i++)
155 {
156 mask.at<uchar>(fill_point[i]) = 1;
157 }
158 }
159 else
160 {
161 for (size_t i = 0; i < fill_point.size(); i++)
162 {
163 mask.at<uchar>(fill_point[i]) = 2;
164 }
165 max_point[k].data = 255;
166 mask.at<uchar>(max_point[k].position) = 255;
167 }
168 }
169 }
C++版本:
2017.10.17更新
1 int FillBlock(Mat src, Mat &mask, Point center)
2 {
3 uchar back = src.at<uchar>(center.y, center.x);
4 vector<Point> fill_point;
5 int count = 0, count_mount = 1;
6 fill_point.push_back(center);
7 while (count < count_mount)
8 {
9 vector<uchar*> img;
10 vector<uchar*> msk;
11 for (int i = -1; i < 2; i++)
12 {
13 img.push_back(src.ptr<uchar>(fill_point[count].y + i));
14 msk.push_back(mask.ptr<uchar>(fill_point[count].y + i));
15 }
16 for (size_t i = 0; i < 3; i++)
17 {
18 for (int j = -1; j < 2; j++)
19 {
20 if (img[i][fill_point[count].x + j] == back && !(j == 0 && i == 1) && msk[i][fill_point[count].x + j] == 255)
21 {
22 fill_point.push_back(Point(fill_point[count].x + j, fill_point[count].y + i - 1));
23 msk[i][fill_point[count].x + j] = 1;
24 }
25 }
26 }
27 msk[1][fill_point[count].x] = 1;
28 count_mount = fill_point.size() - 1;
29 fill_point.erase(fill_point.begin());
30 }
31 return 0;
32 }
33
34 void MaskImage(InputArray _src, Mat &mask)
35 {
36 Mat src = _src.getMat(),mask_tmp = Mat::zeros(src.size(), CV_8UC1);
37 mask_tmp.setTo(255);
38 Mat rows = Mat::zeros(Size(src.cols, 1), CV_8UC1), cols = Mat::zeros(Size(1, src.rows), CV_8UC1);
39 Mat src_rows_beg = mask_tmp.row(0);
40 Mat src_rows_end = mask_tmp.row(src.rows - 1);
41 Mat src_cols_beg = mask_tmp.col(0);
42 Mat src_cols_end = mask_tmp.col(src.cols - 1);
43 rows.copyTo(src_rows_beg); rows.copyTo(src_rows_end);
44 cols.copyTo(src_cols_beg); cols.copyTo(src_cols_end);
45
46 for (size_t i = 3; i < src.rows-3; i++)
47 {
48 vector<uchar*> img;
49 uchar* msk = mask_tmp.ptr(i);
50 uchar* img1 = src.ptr(i);
51 for (int k = -3; k < 4; k++)
52 {
53 img.push_back(src.ptr<uchar>(k + i));
54 }
55 for (size_t j = 3; j < src.cols-3; j++)
56 {
57 bool flag = false;
58 if (msk[j] != 255) continue;
59 float sum[4] = { 0 };
60 sum[0] = img[3][j];
61 Point x0 = Point(0, 0);
62 Point x1 = Point(0, 0);
63 uchar n_count = 0;
64 for (int m = 2; m < 5; m++)
65 {
66 for (int n = -1; n < 2; n++)
67 {
68 if (m == 3 && n == 0) continue;
69 sum[1] = sum[1] < img[m][j + n] ? img[m][j + n] : sum[1];
70 x0 = sum[0] == img[m][j + n] ? Point(m, n) : x0;
71 n_count = sum[0] == img[m][j + n] ? n_count+1 : n_count;
72 //flag = img[3][j + 0] < img[m][j + n] ? true : flag;//如果目标像素小于周围任何一个像素就说明这个一定不是最大值
73 }
74 }
75 for (int m = 1; m < 6; m++)
76 {
77 for (int n = -2; n < 3; n++)
78 {
79 if (2 <= m && m <= 4 && -1 <= n && n <= 1) continue;
80 sum[2] = sum[2] < img[m][j + n] ? img[m][j + n] : sum[2];
81 x1 = sum[0] == img[m][j + n] ? Point(m, n) : x1;
82 n_count = sum[0] == img[m][j + n] ? n_count+1 : n_count;
83 //flag = img[3][j + 0] < img[m][j + n] ? true : flag;//如果目标像素小于周围任何一个像素就说明这个一定不是最大值
84 }
85 }
86 for (int m = 0; m < 7; m++)
87 {
88 for (int n = -3; n < 4; n++)
89 {
90 sum[3] = sum[3] < img[m][j + n] && !(1 <= m && m <= 5 && -2 <=n && n <= 2) ? img[m][j + n] : sum[3];
91 //flag = img[3][j+0] < img[m][j + n] ? true : flag;
92 }
93 }
94 x1 = (x1.x == 0 && x1.y == 0) || n_count >= 3 ? x0 : x1;
95 int tmp = sum[0] >= sum[1] && sum[1] >= sum[2] && sum[2] >= sum[3] && (abs(x0.x - x1.x) <= 2 && abs(x0.y - x1.y) <= 2)
96 ? 0 : FillBlock(src, mask_tmp, Point(j, i));
97 }
98 }
99 mask = mask_tmp.clone();
100 }
101
102 vector<Point> Find_Max(InputArray _src, Mat&mask,int gap,int radius)
103 {
104 Mat src = _src.getMat();
105
106 typedef struct MyStruct
107 {
108 Point position;
109 float data;
110 }MyStruct;
111
112 MaskImage(src, mask);
113 vector<MyStruct> max_point;
114 for (size_t i = 0; i < src.rows; i++)
115 {
116 uchar *img = src.ptr<uchar>(i);
117 uchar *msk = mask.ptr<uchar>(i);
118 for (size_t j = 0; j < src.cols; j++)
119 {
120 if (msk[j] != 255) continue;
121 MyStruct my_struct;
122 my_struct.data = img[j];
123 my_struct.position = Point(j, i);
124 max_point.push_back(my_struct);
125 }
126 }
127 for (size_t i = 0; i < max_point.size(); i++)
128 {
129 for (size_t j = i; j < max_point.size(); j++)
130 {
131 MyStruct temp;
132 if (max_point[i].data <= max_point[j].data)
133 {
134 if (max_point[j].data == 0) continue;
135 temp = max_point[j];
136 max_point[j] = max_point[i];
137 max_point[i] = temp;
138 }
139 }
140 }
141
142 for (size_t k = 0; k < max_point.size(); k++)//---
143 {
144 uchar back = src.at<uchar>(max_point[k].position.y, max_point[k].position.x);
145 vector<Point> fill_point;
146 int count = 0, count_mount = 1;
147 fill_point.push_back(max_point[k].position);
148
149 while (count < count_mount && max_point[k].data != 1)
150 {
151 vector<uchar*> img;
152 vector<uchar*> msk;
153 for (int i = -1; i < 2; i++)
154 {
155 img.push_back(src.ptr<uchar>(fill_point[count].y + i));
156 msk.push_back(mask.ptr<uchar>(fill_point[count].y + i));
157 }
158 for (int i = 0; i < 3; i++)
159 {
160 for (int j = -1; j < 2; j++)
161 {
162 //---
163 double x = pow((max_point[k].position.x - fill_point[count].x + j), 2); //(max_point[k].position.x - img[i][fill_point[count].x + j])*(max_point[k].position.x - img[i][fill_point[count].x + j]);
164 double y = pow((max_point[k].position.y - (fill_point[count].y + i - 1)), 2); // (max_point[k].position.y - img[i][fill_point[count].y + j])*(max_point[k].position.y - img[i][fill_point[count].x + j]);
165 int distance = sqrt(x + y);
166 if (img[i][fill_point[count].x + j] <= img[0][fill_point[count].x] - gap
167 || msk[i][fill_point[count].x + j] == 3
168 //|| msk[i][fill_point[count].x + j] == 0
169 || (j == 0 && i == 1)
170 || distance >= radius) continue;
171 if (img[i][fill_point[count].x + j] == 2 || msk[i][fill_point[count].x + j] == 0) max_point[k].data = 1;
172 msk[i][fill_point[count].x + j] = 3;
173 fill_point.push_back(Point(fill_point[count].x + j, fill_point[count].y + i - 1));
174 count_mount++;
175 }
176 }
177 count++;
178 }
179 if (max_point[k].data == 1)
180 {
181 for (size_t i = 0; i < fill_point.size(); i++)
182 {
183 mask.at<uchar>(fill_point[i]) = 1;
184 }
185 }
186 else
187 {
188 for (size_t i = 0; i < fill_point.size(); i++)
189 {
190 mask.at<uchar>(fill_point[i]) = 2;
191 }
192 max_point[k].data = 255;
193 mask.at<uchar>(max_point[k].position) = 255;
194 }
195 }
196 vector<Point> print_wjy;
197 for (size_t i = 0; i < mask.rows; i++)
198 {
199 uchar *msk = mask.ptr<uchar>(i);
200 for (size_t j = 0; j < mask.cols; j++)
201 {
202 if (msk[j] == 255)
203 print_wjy.push_back(Point(j, i));
204 }
205
206 }
207 return print_wjy;
208 }