标签：org 几何 index pts time() else filter 排除 ges
局部极值提取算法
这是http://www.imagepy.org/的作者原创，我只是对其理解之后改进和说明，欢迎大家使用这个小软件！
算法原理：（按照程序思路来，当然其中很多可以改进的地方）
第一步：
先进行距离变换（这类问题都是要通过几何操作进行，距离变换就是几何距离）
第二步：
先找全局可能最大值，然后对图像进行标记
全局可能最大值是基础，3X3核在图像上滑动去找全局可能最大值。标记是为了掩膜操作更方便，你看opencv很多函数都有掩膜，PS上面也有那个白色的膜。
图片结果是：背景和最大值标记3，前景标记1，图像边界标记0（为了防止越界）
貌似下面的图片显示有问题，matploylib的显示有时候确实有问题，Debug的时候就和显示不一样，没关系这一步比较简单。
第三步：
对全局可能最大值进行排除，这一步是程序核心！这里以求最大值为例，最小值相反。
下面以图为例进一步说明三个原则：
不重合原则：
下图A1和A2有领土纠纷，但是A1峰值高所以先进行领地划分，当A2进行领地划分的时候一旦接触A1的领地就立马被消灭，最后只剩A1的存在。
不包括原则：
下图A1先进行领土划分，过程中直接把A2消灭了，不存在A2的领土划分
领主大于等于领地原则，也可以叫同等峰值先来先得原则：
由于A1先进行领土划分，所以后面的都被消灭了，读到这里你会发现这个程序的BUG，后面再详细谈论这一点
程序有待改进：
这部分可以通过处理---“同样的峰值且山脉相连”的一类极值点，对其平均、覆盖范围最多等操作取一个中间的值！
由于python用的不熟练，而且现在时间太晚了脑子转不动，以后用C++写的时候再改进这一点。
1 import scipy.ndimage as ndimg 2 import numpy as np 3 from numba import jit 4 import cv2 5 6 7 # 获得邻居像素坐标，相当于C++的核（支持多维度） 8 def neighbors(shape): 9 dim = len(shape) 10 block = np.ones([3] * dim) 11 block[tuple([1] * dim)] = 0 12 idx = np.where(block > 0) 13 idx = np.array(idx, dtype=np.uint8).T 14 idx = np.array(idx - [1] * dim) 15 acc = np.cumprod((1,) + shape[::-1][:-1]) 16 return np.dot(idx, acc[::-1]) 17 18 19 @jit # trans index to r, c... 20 #坐标转换，一维-->>二维 21 def idx2rc(idx, acc): 22 rst = np.zeros((len(idx), len(acc)), dtype=np.int16) 23 for i in range(len(idx)): 24 for j in range(len(acc)): 25 rst[i, j] = idx[i] // acc[j] 26 idx[i] -= rst[i, j] * acc[j] 27 return rst 28 29 30 #@jit # fill a node (may be two or more points) 31 #这里和寻找山脊线一个思路，不过这里是填充相同水平位置 32 #不懂的可以看看我注释的山脊线算法 33 def fill(img, msk, p, nbs, buf): 34 msk[p] = 3 35 buf[0] = p 36 back = img[p] 37 cur = 0 38 s = 1 39 while cur < s: 40 p = buf[cur] 41 for dp in nbs: 42 cp = p + dp 43 if img[cp] == back and msk[cp] == 1:#对周围相等的像素进行扩充赋值为3 44 msk[cp] = 3 45 buf[s] = cp 46 s += 1 47 #防止buf太大越界 48 if s == len(buf): 49 buf[:s - cur] = buf[cur:] 50 s -= cur 51 cur = 0 52 cur += 1 53 #msk[p] = 3 54 55 56 #@jit # my mark 57 #此时的图片边界四个边界为0，中间是1 58 #图片结果是：背景3，前景1，最大值3 59 def mark(img, msk, buf, mode): # mark the array use (0, 1, 2) 60 omark = msk #查看演示的步骤，因为是浅拷贝所以不需要再去升维度 61 nbs = neighbors(img.shape) 62 idx = np.zeros(1024 * 128, dtype=np.int64) 63 img = img.ravel() # 降维 64 msk = msk.ravel() # 降维 65 s = 0 66 for p in range(len(img)): 67 if msk[p] != 1: continue # 防止越界 68 flag = False # 这里改进，找平原和最高峰（背景和最大值） 69 for dp in nbs: 70 if mode and img[p + dp] > img[p]: #这里改进 71 flag = True 72 break 73 elif not mode and img[p + dp] < img[p]: 74 flag = True 75 break 76 if flag : continue 77 else : fill(img, msk, p, nbs, buf) 78 idx[s] = p 79 s += 1 80 if s == len(idx): break 81 plt.imshow(omark, cmap=‘gray‘) 82 return idx[:s].copy() 83 84 85 #@jit # 3 最大值标记 2 最终殖民地标记 4 临时殖民地标记 86 def filter(img, msk, idx, bur, tor, mode): 87 omark = msk # 查看演示的步骤，因为是浅拷贝所以不需要再去升维度 88 nbs = neighbors(img.shape) 89 acc = np.cumprod((1,) + img.shape[::-1][:-1])[::-1] 90 img = img.ravel() 91 msk = msk.ravel() 92 93 arg = np.argsort(img[idx])[::-1 if mode else 1] #最大值（大到小）和最小值（小到大） 94 95 for i in arg: 96 if msk[idx[i]] != 3: #如果被其他占领那就GG 97 idx[i] = 0 98 continue 99 cur = 0 100 s = 1 101 bur[0] = idx[i] #存储第一个位置，后面会自动存储，其作用是存储占领的地盘 102 while cur < s: 103 p = bur[cur] 104 if msk[p] == 2: #如果当前殖民地和别人有冲突，那就GG，因为先占领的峰值更高 105 idx[i] = 0 106 break 107 108 for dp in nbs: 109 cp = p + dp 110 if msk[cp] == 0 or cp == idx[i] or msk[cp] == 4: continue 111 if mode and img[cp] < img[idx[i]] - tor: continue 112 if not mode and img[cp] > img[idx[i]] + tor: continue 113 bur[s] = cp 114 s += 1 115 if s == 1024 * 128: 116 cut = cur // 2 117 msk[bur[:cut]] = 2 118 bur[:s - cut] = bur[cut:] 119 cur -= cut 120 s -= cut 121 122 if msk[cp] != 2: msk[cp] = 4 #临时标记占领的地盘为4 123 cur += 1 124 msk[bur[:s]] = 2 #占领的地盘标记为2 125 #plt.imshow(omark, cmap=‘gray‘) 126 127 return idx2rc(idx[idx > 0], acc) 128 129 130 def find_maximum(img, tor, mode=True): 131 msk = np.zeros_like(img, dtype=np.uint8) # 返回img大小的零矩阵 132 msk[tuple([slice(1, -1)] * img.ndim)] = 1 # slice是切片数据，slice*img.ndim就是在二维或者三维上面进行切片，这里是让边缘为0，内部填充为1 133 buf = np.zeros(1024 * 128, dtype=np.int64) 134 omark = msk 135 idx = mark(img, msk, buf, mode) 136 plt.imshow(msk, cmap=‘gray‘) 137 idx = filter(img, msk, idx, buf, tor, mode) 138 return idx 139 140 141 if __name__ == ‘__main__‘: 142 from scipy.misc import imread 143 from scipy.ndimage import gaussian_filter 144 from time import time 145 import matplotlib.pyplot as plt 146 147 img = cv2.imread(‘123.jpg‘) 148 img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY) 149 ret2, img = cv2.threshold(img, 0, 255, cv2.THRESH_BINARY | cv2.THRESH_OTSU) 150 img[:] = ndimg.distance_transform_edt(img) 151 pts = find_maximum(img, 20, True) 152 start = time() 153 pts = find_maximum(img, 10, True) 154 print(time() - start) 155 plt.imshow(img, cmap=‘gray‘) 156 plt.plot(pts[:, 1], pts[:, 0], ‘y.‘) 157 plt.show()
标签：org 几何 index pts time() else filter 排除 ges
原文地址：http://www.cnblogs.com/wjy-lulu/p/7416135.html
