局部極值提取算法
這是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 }