精华内容
下载资源
问答
  • 这是一个非常快速和简单的局部极值查找器。 它同时计算局部最大值和最小值,并将索引返回到输入向量中的相应值(通过返回第一个点的索引来处理常量区域)。 这不是 Matlab 的 findpeaks 函数的替代品,后者功能更...
  • % 找到真正的局部极值 % ..... x 包含待分析序列的向量, % ..... L 元胞数组 {L(1), L(2)},其中% ..... L(1) 是最大值位置的逻辑向量, % ..... L(2) 是最小值位置的逻辑向量。 在某些情况下,处理时间对用户来说...
  • 局部极值.py

    2019-06-09 11:43:22
    利用python求图像的局部极值,并标注。这是http://www.imagepy.org/的作者原创,我只是对其理解之后改进和说明,欢迎大家使用这个小软件!
  • 返回具有局部极值(最小值、最大值)位置和值的结构。 输出结构要素: minx:最小位置miny:最小值maxx:最大位置maxy:最大值 例子: v=[1 3 5 5 2 -2 -2 -1 8 -10 5]; o1=极值(v) o2=extrems(v,'ysort') o3=...
  • 快速搜索局部极值以获得最大信息系数(MIC)
  • 为求得类似仿真函数的黑箱函数优化问题的全部局部极值点,提出了一种基于适应值曲面分析的新算法。首先,对适应值距离相关系数(fitness distance correlation, FDC)进行了改进,探讨了改进FDC指标与适应值曲面崎岖度...
  • 本文提出了一种新颖的局部极值动态系统,该系统中的每个玩家只需要识别邻居的收益,并在邻居获得最低收益时随机改变其策略即可。 在此更新规则下,针对不同大小的邻域(以其相应的半径r为特征),探讨了合作的演变...
  • 局部极值提取算法 这是http://www.imagepy.org/的作者原创,我只是对其理解之后改进和说明,欢迎大家使用这个小软件! 有需要C++的朋友,可有留邮箱,经测试和改进的版本! 算法原理:(按照程序思路来,当然其中很...

    局部极值提取算法

    这是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 }
    
    View Code

    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 }
    
    展开全文
  • % 约束点可以限制为局部极值。 % [pfunc,pval]=polyfitcon(x,y,n,xc,yc,p0,isLocExt) 返回%pfunc,存储用于拟合的阶数'n'的多项式。 % pval,拟合多项式的权重。 %pfunc = @(x,pval)pval(1)。* x。^(n-1)+...
  • Persistence1D是用于在一维数据中查找局部极值及其持久性的类。 局部最小值和局部最大值根据它们的持久性被提取、配对和排序。 代码在 O(n log n) 时间内运行,其中 n 是输入点的数量。 请检查doc文件夹以获取详细的...
  • 本次实验利用了基于局部极值的分水岭算法来实现圆斑点的检测。在OPENCV中提供了simpleBlobDetector特征检测器来实现这种斑点检测算法,正如它的名称,该算法使用最简单的方式来检测斑点类的特征点,效果较好,设置...

    本次实验利用了基于局部极值的分水岭算法来实现圆斑点的检测。在OPENCV中提供了simpleBlobDetector特征检测器来实现这种斑点检测算法,正如它的名称,该算法使用最简单的方式来检测斑点类的特征点,效果较好,设置较为宽松的参数,就能取得较好的效果。

    一 算法的步骤

    第一步 多次二值化图像

    首先通过一系列连续的阈值把输入的灰度图像转换为一个二值图像的集合,阈值范围为[T1, T2],步长为t,则所有阈值为:
    T1,T1+t,T1+2t,T1+3t,……,T2

    第二步 确定候选圆点

    通过检测每一幅二值图像的边缘的方式提取出每一幅二值图像的连通区域,我们可以认为由边界所围成的不同的连通区域就是该二值图像的斑点。其中,并不是所有的二值图像的连通区域都可以认为是二值图像的斑点,我们通过一些限定条件来得到更准确的斑点。这些限定条件包括颜色,面积和形状,斑点的形状又可以用圆度,偏心率,或凸度来表示,具体参数和计算方法查看参数部分介绍。

    第三步 对图像圆点进行分类,确定目标圆点

    根据所有二值图像斑点的中心坐标对二值图像斑点进行分类,从而形成灰度图像的斑点,属于一类的那些二值图像斑点最终形成灰度图像的斑点,具体来说就是,灰度图像的斑点是由中心坐标间的距离小于阈值的那些二值图像斑点所组成的,即这些二值图像斑点属于该灰度图像斑点;

    二 参数部分

    filterByColor斑点颜色。

    filterByArea斑点面积。连通区域的面积太大和太小都不是斑 。所以我们需要计算连通区域的面积,只有当该面积在我们所设定的最大面积和最小面积之间时,该连通区域才作为斑点被保留下来。

    filterByCircularity 斑点圆度。任意形状的圆度C定义为:
    这里写图片描述
    其中,S和p分别表示该形状的面积和周长,当C为1时,表示该形状是一个完美的圆形,而当C为0时,表示该形状是一个逐渐拉长的多边形。

    filterByInertia斑点惯性率。偏心率是指某一椭圆轨道与理想圆形的偏离程度,长椭圆轨道的偏心率高,而近于圆形的轨道的偏心率低。圆形的偏心率等于0,椭圆的偏心率介于0和1之间,而偏心率等于1表示的是抛物线。直接计算斑点的偏心率较为复杂,但利用图像矩的概念计算图形的惯性率,再由惯性率计算偏心率较为方便。偏心率E和惯性率I之间的关系为:
    这里写图片描述
    因此圆形的惯性率等于1,惯性率越接近1,圆形的程度越高(在实际操作中并没有用到)。
    filterByConvexity 表示斑点凸度。在平面中,凸形图指的是图形的所有部分都在由该图形切线所围成的区域的内部。我们可以用凸度来表示斑点凹凸的程度,凸度V的定义为:
    这里写图片描述
    其中,H表示该斑点的凸壳面积
    minDistBetweenBlobs斑点距离的最小阈值,中心坐标间的距离小于该值就认为为同一个斑点。

    三 源码

    #include "opencv2/opencv.hpp"
    #include "opencv2/core/core.hpp"
    #include "opencv2/highgui/highgui.hpp"
    #include "opencv2/imgproc/imgproc.hpp"
    #include "opencv2/features2d/features2d.hpp" 
    
    
    using namespace cv;
    using namespace std;
    
    int main(int argc, char* argv[])
    {
        Mat srcImage = imread("pos71_run7_img190.jpg");
        Mat srcGrayImage, srcGrayImage1, srcImage1, keyPointImage11, keyPointImage21;
        if (srcImage.channels() == 3)
        {
            cvtColor(srcImage, srcGrayImage, CV_RGB2GRAY);
        }
        else
        {
            srcImage.copyTo(srcGrayImage);
        }
        GaussianBlur(srcGrayImage, srcGrayImage1,Size(11,11),0,0);
        vector<KeyPoint>detectKeyPoint;
        Mat keyPointImage1, keyPointImage2;
    
        SimpleBlobDetector::Params params;
    
        params.filterByArea = true;
        params.minArea = 15;
        params.maxArea = 1500;
    
        params.thresholdStep = 10;    //二值化的阈值步长,即公式1的t
        params.minThreshold = 50 ;
        params.maxThreshold = 200;
    
        params.filterByCircularity = true;  //斑点圆度
        params.minCircularity = 0.8;
        params.maxCircularity = 1;
    
        params.filterByConvexity = false;    //斑点凸度
        params.minConvexity = (float)0.5;
        params.maxConvexity = 1;
    
        params.filterByInertia = false;  //斑点惯性率
        params.minInertiaRatio = (float)0.4;
        params.maxInertiaRatio = 1;
    
        params.filterByColor = true;
        params.blobColor = 255;//白色
    
        params.minRepeatability = 2;//重复的最小次数,
        params.minDistBetweenBlobs = 20; //最小的斑点距离
        Ptr<SimpleBlobDetector> sbd = SimpleBlobDetector::create(params);
        sbd->detect(srcGrayImage1, detectKeyPoint);
        drawKeypoints(srcImage, detectKeyPoint, keyPointImage1, Scalar(0, 255, 0), DrawMatchesFlags::DRAW_RICH_KEYPOINTS);
    
        resize(srcImage, srcImage1, Size(srcImage.cols / 4, srcImage.rows / 4), 0, 0, INTER_LINEAR);
        imshow("src image", srcImage1);
        resize(keyPointImage1, keyPointImage11, Size(keyPointImage1.cols /4, keyPointImage1.rows /4), 0, 0, INTER_LINEAR);
        imshow("keyPoint image1", keyPointImage11);
        imwrite("spos71_run7_img190.jpg", keyPointImage1);
    
        waitKey(0);
        return 0;
    }

    四 实验结果

    这里写图片描述
    ![这里写图片描述]
    这里写图片描述

    五 补充说明

    保存特征点,绘制特征点,在使用中是不透明的,下面做一个简要介绍。

    1.KeyPoint特征点类
    保存特征点各种信息的KeyPoint类在使用的主要属性:

    class KeyPoint
    {
    Point2f  pt;  //特征点坐标
    float  size; //特征点邻域直径
    float  angle; //特征点的方向,值为0~360,负值表示不使用
    float  response; //特征点的响应强度,代表了该点是特征点的程度,可以用于后续处理中特征点排序
    int  octave; //特征点所在的图像金字塔的组
    int  class_id; //用于聚类的id
    }

    主要包含的特征点信息有:位置、邻域直径、特征的方向、响应强度、多尺度信息和分类等。特征点匹配的实现就是通过逐个匹配特征点的这些信息。

    2.drawKeypoints特征点绘制

    opencv提供了一个快速绘制特征点的函数drawKeypoints,函数原型:

    void drawKeypoints( const Mat& image, const vector<KeyPoint>& keypoints, CV_OUT Mat& outImage,  const Scalar&color=Scalar::all(-1), int flags=DrawMatchesFlags::DEFAULT );  

    第一个参数image:原始图像,可以使三通道或单通道图像;
    第二个参数keypoints:特征点向量,向量内每一个元素是一个KeyPoint对象,包含了特征点的各种属性信息;
    第三个参数outImage:特征点绘制的画布图像,可以是原图像;
    第四个参数color:绘制的特征点的颜色信息,默认绘制的是随机彩色;
    第五个参数flags:特征点的绘制模式,其实就是设置特征点的那些信息需要绘制,那些不需要绘制,有以下几种模式可选:
      DEFAULT:只绘制特征点的坐标点,显示在图像上就是一个个小圆点,每个小圆点的圆心坐标都是特征点的坐标。
      DRAW_OVER_OUTIMG:函数不创建输出的图像,而是直接在输出图像变量空间绘制,要求本身输出图像变量就 是一个初始化好了的,size与type都是已经初始化好的变量
      NOT_DRAW_SINGLE_POINTS:单点的特征点不被绘制
      DRAW_RICH_KEYPOINTS:绘制特征点的时候绘制的是一个个带有方向的圆,这种方法同时显示图像的坐,size,和方向,是最能显示特征的一种绘制方式。

    六 参考

    http://blog.csdn.net/zhaocj/article/details/44886475
    http://blog.csdn.net/dcrmg/article/details/52553513

    展开全文
  • 在课堂上,我们介绍了许多边缘感知过滤:双边,WLS,局部极值,扩散图,域变换,局部拉普拉斯算子,L0最小化和引导滤波器。 在此作业中,您有三个选择。 对于第一个选项,您必须实现“本地极值”或“扩散图”(未...
  • 皇后问题局部极值启发式搜索方法 背景说明:皇后问题是算法领域的著名问题,本篇主要是对皇后问题背景下的局部搜索问题的研究。

    8-Queens Problem皇后问题局部极值启发式搜索方法


    Backgrounds

    背景说明:

    皇后问题是算法领域的著名问题,问题的背景是在8*8的国际象棋棋盘上摆满8个皇后使之互相不能攻击(由于皇后在国际象棋规则中横、纵、斜三个方向均可以移动,即摆放要使得皇后两两不在同一横、纵、斜线上)。放松这个问题的变长条件,我们可以将其扩大至N*N的棋盘上摆放N个皇后,使之各自在三个方向上不重复出现。

    大一时刚刚学习编程语言的我们,在初次接到这道题时,主要是运用了递归函数和回溯的思想去解决问题,通过深度优先搜索(depth-first search, DFS)策略,依次增加棋盘上可以摆放的位置。若遇到第k步棋子摆放的所有可能位置都发生冲突,则退回到第k-1步,修改k-1的状态,以达到期望的最终解。


    State space of 8-queens problem

    皇后问题的状态空间:

    状态空间,这里理解成所有摆放的可能性情况。

    皇后问题广义的解的状态空间,在没有过多限制条件的情况下是非常庞大的。

    在64个互不重复的位置上选择无差别的8个位置,总计应该有C(64, 8)=4,426,165,368 (44亿)大小的状态空间. 图示为较为集中的分布情况和可以完全分散开的分布情况。

      

            图 1.1 皇后分布紧凑          图 1.2 皇后分布松散 

    注意这里C(N*N, N)表示在N*N种方案中选择N个的方案数

    不同约束条件下的不同状态空间规模:

    实际上,全局解状态空间的大小也由我们的限制条件决定,下面提出两种解决方案。但总体上策略一致地是,我们用一个长度为8的一维数组表示每一列皇后的位置,这种表示方法首先就压缩掉了很大一部分的解的状态空间——因为相当于即使在横向上有皇后可能会重合,这种每列一个位置数的表示方式放弃了同一列上出现多个皇后的情况。图示为横纵均不重合分散分布和允许行重复的分布情况。  

      

           图 2.1 允许同行冲突 图 2.2 保证横纵均不冲突


    但上述两种表示方法出于求解出尽可能多的局部极值和尽快求解出最优解的考虑,依然是稳妥的。此时状态空间的大小至多已经被压缩至8^8=16,777,216 (1670万)种,相比于44亿已经极大地压缩了解的数量。

    (1).第一种表示方法,我们在允许行重合的情况下,每次操作只在该列上对棋子进行上下移动,并考虑每次移动对于总冲突数的影响。

    此种排列的状态空间为8^8=16,777,216.

    (2).第二种表示方法,我们对1~8的数字进行全排列,得到{a1, a2, …, ai, …, a8}用于表示第i 列上的皇后位置,这样可以保证行列均不重复。

    此种排列的状态空间为8!=40,320. 


    局部极值和全局最值的理解:

    全局最值,是对于每一个搜索问题需要找到的最终结果。对于不同的问题,其问题的离散/连续性可能不同,解的数量、状态分布也不相同。皇后问题本身是一个离散问题(解的状态空间数量是有限个数的,虽然解空间总规模可能非常大),且解的数量也是巨大的。至少对于8皇后,不存在唯一的状态。后续运行程序可知,在不考虑翻转、旋转等价的情况下,有92组互不相同的解。即在横纵均不重合的40,320规模的状态空间(或在更大范围的状态空间中)有92组全局最值。

    局部极值,是我们通过某种途径趋向全局最值时经过的、不能再用当前算法得出相对更进一步的结果(而必须由重启、重新打乱而继续开始一次新的搜索)时,“卡住”的状态。


    图 3 搜索路径与全局最优、局部最优的关系

    (state space为构成搜索路径的n维解空间的子集,objective function为目标函数)

    对于同一个问题,我们考虑f(x1, x2, …, xi, …, xn)=K的目标函数,希望通过启发式搜索的方式寻找全局K的最优值,即在n维解空间上,我们得到由目标函数f(.)决定的决策-目标曲面。在这个n+1维曲面S(X<x1, x2, …, xn>, K)上,不同的约束条件和算法就相当于通过不同的n维路径Ls<x1s, x2s, …, xns>去逼近全局最优K.

    其中Ls是某种启发式搜索算法的搜索路径,对于特定的搜索算法,Ls是n维解状态空间的有向路径。记Ls上的每一步解为Xj, 并设该条路径从某一初始状态X0出发,经过p步可以到达全局或局部最优值(其中任意Xj都是n维状态空间上的点),则该路径Ls可以表示为Ls={X0, X1, …, Xj, …, Xp}. 由于状态空间的不平整性,显然通过不同的算法定义到的局部最优解(局部极值)也是不同的。

    注意这里的定义包括后续目标函数的定义、估值函数的定义、最终的求解目的等等


    Evaluation function

    估值函数:

    对于启发式搜索,我们需要在巨大的解的状态空间中尽可能快地找出尽可能多的全局最优解,这就需要我们在做每一步选择(选择继续扩大解的已知范围,使之尽快铺展至n维或跳转至下一个解的暂留状态)时,能有效评判该解对于“发现最优解”这一目的的贡献意义的函数。或者我们也可以这样理解:搜索算法在在搜索路径的每一个解Xj停留时,并不知道该解是否处于到达全局最优的必经路径上,因此搜算算法会根据一个估值函数,对“该步操作可以到达全局最优或局部最优”这一事件做出概率分析。当搜索路径Ls走到Xj-1这一步时,会对所有可能的后续状态集合{Xj’|Xj’=Xj’1, Xj’2,…}进行估值,最后选出估值最大的一个Xj’作为Xj下一步。后续也据此方法继续走下去。

    当最终走到某一步,估值函数对所有可能的后继都没有最优估值或正向估值时,就说明搜索路径已经到达极值。此时检查目标函数K即可判断当前解是否是全局最优。

    需要说明的是,对于全局最优,我们也可以有不同的定义方法。在本问题中,我们选择用“冲突行、列或斜线上棋子数减一”表示冲突情况,即当某行、列或斜线上皇后棋子数目为0个或1个时,认为没有冲突;皇后棋子数目为2个及以上时,认为冲突数是皇后数目-1. 

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    int getChange(int i, int j) {
        int pi = position[i], pj = position[j];
        // i-pi, j-pj
        // left: i+pi==j+pj, right: i-pi==j-pj
        // left[i+pi-1] left[j+pj-1] right[N-(i-pi)] right[N-(j-pj)]
        if (i+pi == j+pj) {
            // same left line
            return ((left[i+pi-1]>=3)+1)
            + (right[N-(i-pi)]>=2)
            + (right[N-(j-pj)]>=2)
            - ((right[N-(i-pj)]>=1)+1)
            - (left[i+pj-1]>=1)
            - (left[j+pi-1]>=1);
        else if (i-pi == j-pj) {
            // same right line
            return ((right[N-(i-pi)]>=3)+1)
            + (left[i+pi-1]>=2)
            + (left[j+pj-1]>=2)
            - ((left[i+pj-1]>=1)+1)
            - (right[N-(i-pj)]>=1)
            - (right[N-(j-pi)]>=1);
        else {
            // (i, j) no collision
            return (left[i+pi-1]>=2)
            + (left[j+pj-1]>=2)
            + (right[N-(i-pi)]>=2)
            + (right[N-(j-pj)]>=2)
            - (left[i+pj-1]>=1)
            - (left[j+pi-1]>=1)
            - (right[N-(i-pj)]>=1)
            - (right[N-(j-pi)]>=1);
        }
    }
        

    图 4 估值函数

    当发现两列(i, j) 可以被操作时,我们通过计算“交换后”与“交换前”两种状态(两个不同的解的停留点)之间目标函数即总冲突数的变化量,判断该步操作是否应被执行。其中我们考虑了交换前的(i, j)两列处于左斜向冲突(交换后存在右斜向冲突)、右斜向冲突(交换后存在左斜向冲突)、无冲突(交换后彼此依然无冲突)三种情况。 

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    void statistic() {
        clear();
        for (int i=1; i<=N; ++i) {
            int pi=position[i];
            left[i+pi-1]++;
            right[N-(i-pi)]++;
        }
    }
         
    void reCount() {
        statistic();
        Collision = 0;
    leftCollision = 0, rightCollision = 0;
        for (int i=1; i<=2*N-1; ++i) {
            if (left[i]>=2) leftCollision += left[i]-1;
            if (right[i]>=2) rightCollision += right[i]-1;
        }
        Collision = leftCollision + rightCollision;
    }

    图 5 目标函数定义,显然Collision=0时取得全局最优解


    Main Algorithm

    状态空间压缩:

    我们选择交换不重复的列来完成全局最优的搜索结果,这样可以在本身就比较小的状态空间(40,320种)范围内得出结果。并且对于全局最优来说,由于所有的全局最优解均包含在该状态空间中,全局最优的分布率也最高。

    重启策略:

    初始化时生成1到8的升序列并进行时间复杂度O(n)的遍历,每次随机出一个位置,将其与当前扫到的位置做一次交换。之后每次重启仅对现有列进行一次随机交换即可得到一个新的解空间上的X0作为新的一条搜索路径Ls的起点。

    搜索策略: 

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    public:
        void solve(bool ptl) {
            bool moveOn = true;
            int EFT = 0;
            rebootCount = 0;
            reCount(); // initialization process
            startCTime = clock();
            while (!finish()) { // Terminate
                reboot();
                rebootCount++;
                reCount();
                EFT=0; moveOn=true;
                while (moveOn) {
                    moveOn = false;
                    for (int i=1; i<=N; ++i) {
                        for (int j=1; j<i; ++j) {
                            // -------------beginforeach i, j
                            EFT = 0;
                            // for any 'i' col & 'j' col
                            if (collapsed(i, j)) { // (i, j) pair collapsed
                                EFT = getChange(i, j);
                            else EFT = -1;
                            if (EFT>0) { // still trends to minimal
                                moveOn = true;
                                // do position swapping
                                // maintain the array 'left' & 'right'
                                execSwap(i, j);
                            }
                            reCount();
                        // -------------------------- endforeach i, j
                        }
                    }
                }
                if (ptl) {
                    puff();
                    prinf();
                }
            }
            endCTime = clock();
            costCTime = endCTime - startCTime;
        }

    图 6 主要搜索方法

    用两个主要变量控制搜索进程:

    int EFT用来衡量每次可交换列的对(i, j)交换后的估值函数值;

    bool moveOn判断当前循环是否有效,即是否有正向估值产生。

    每当一步搜索不能产生正向估值,即无论如何移动也不能使总的冲突数减小时,一个新的局部极值就可能被找到了。这时我们若判断出总冲突数降至0则可以进一步说明这个极值也是一个全局最优。


    Local search strategy

    贪心与模拟退火:

    爬山问题中,我们判断当前搜索方向是否还会继续趋向于更大值。若当前搜索方向上不再出现上升态势,搜索会停止并驻留在当前解上(即局部最优解),所以爬山策略是一种严格贪心算法,不保证可以找到全局最优。

    模拟退火算法正是通过增加一定几率的反向跳转趋势,使得搜索位置得以“逃脱”局部最优,获得移动到全局最优的可能;通过控制反向跳转几率的收敛,使得解有更大可能性最终停留在全局最优解上。

    局部搜索策略:

    皇后问题中,局部最优解的分布较为稠密,全局最优解的比例也相对较高。这使得模拟退火等一般意义上的概率收敛不一定能取得较好的效果。


    Simple Statistic

    简单统计结果:

    下面给一个运行的例子 

    50000次运行求出最优解,平均每个最优解的求得只需要2.6次重启,即经过2.6个局部最优解到达了全局最优;启动最多的也只有25次;平均耗时极短。  


    程序打印出的一些局部极值示意图微笑




    ps. 最近新开通了博客 会放一些学习过程中的心得体会(和一些作业或者报告)

    如果你觉得有用,欢迎分享,谢谢!

    展开全文
  • 本文提出了一种利用不断调大学习率的方法试图跳出SGD优化过程中的局部极小值或者鞍点的方法,并提出了两种具体的方法:随机漫步起跳法和历史最优起跳法,实验证明相对常规优化方法有一定性能提升。

    /* 版权声明:可以任意转载,转载时请标明文章原始出处和作者信息 .*/

                                                  

                                                                                             黄通文   张俊林( 2016年12月)

                                                                


    注:这篇文章的思路以及本文是我们在2016年底左右做自动发现探索网络结构过程中做的,当时做完发现ICLR 2017有类似论文提出相似的方法,所以没有做更多实验并就把这篇东西搁置了。今天偶然翻到,感觉还是有一定价值,所以现在发出来。


    对于神经网络目标函数优化来说,在参数寻优过程中最常用的算法是梯度下降,目前也出现了很多基于SGD基础上的改进算法,比如带动量的SGD,Adadelta,Adagrad,Adam,RMSProp等梯度下降改进方法,这些算法大多都是针对基础更新公式进行改进:

       


          其中,w是训练过程中的当前参数值,alpha是学习率,grad(w)是此刻的梯度方向。目前的改进算法要么修正学习率alpha要么修正梯度计算方法grad(w),要么两者都修正,比如Adagrad/Adam动态调整学习率,Adam动态调整梯度的更新量或者有些算法加入动量部分。


     一般深度神经网络由于很深的深度以及非线性函数这两个主要因素,导致目标的损失函数是非凸函数,该损失函数曲面中包含许多的驻点,驻点可能是以下三类关键点中的一种:鞍点,局部极小值,全局极小值,其中全局极小值是我们期望训练算法能够找到的,而其它两类是期望能够避免的。一个示意性的损失函数曲面如下图所示:



    不少学习算法在训练过程中,随着误差的减少,迭代次数的增加,步长变化越来越小,训练误差越来越小直到变为零,局部探索到一个驻点,如果这个点是误差曲面的鞍点,即不是最大值也不是最小值的平滑曲面,那么一般结果表现为性能比较差;如果这个驻点是局部极小值,那么表现为性能较好,但不是全局最优值。目前大多数训练算法在碰到无论是鞍点还是局部极值点的时候,因为此刻学习率已经变得非常小,所以会陷入非全局最优值不能自拔。这其实就是目前采用SGD类似的一阶导数方法训练机器学习系统面临的一个很重要的问题。


    目前来看,没有方法能够保证找到非凸函数优化过程中的全局最小值,也就是说,不论你怎么调参将性能调到你能做到的最好,仍然大概率获得了优化过程中的一个鞍点或者局部极小值,那么一个很自然的想法就是:能否在优化过程的末尾,就是已经大概率陷入某个局部极小值或者鞍点的时候,放大学习率,强行让训练算法跳出当前的局部最小值或者鞍点,所谓“世界这么大,我想去看看”,继续寻找可能效果会更好的极值点,这样多次跳出优化过程中的局部最小值,然后比较其中性能最好的极小值,并以那个时刻的参数作为模型参数。当然,这仍然没有办法保证找到全局最小值,只是说多比较一些极值点,矮子里面拔将军的意思。这个思路是去年我们在利用人工生命和遗传算法探索自动生成神经网络结构的不断优化网络性能过程中产生的,当时是看到了SNAPSHOT ENSEMBLES : TRAIN 1,GET M FOR FREE这篇论文受到启发动了这个念头,于是单独测试了一下这个思路,结果表明一般情况下,这种方法是能够有效提升优化效果的,不过提升幅度不算太大(后来看到了这篇论文:SGDR: STOCHASTIC GRADIENT DESCENT RESTARTS 发现基本思路是有点类似的)。


    一.两种跳出驻点的方法


    上面说了基本思路,其实思路很朴实也很简单:常见的训练方法在刚开始的时候学习率设置的比较大,这样可以大幅度快步向极值点迈进,加快收敛速度,随着越来越接近极值点,逐步将学习率调小,这样避免步子太大跳过极值点或者造成优化过程震荡结果不稳定。当优化到性能无法再提升的时候,此时大概率陷入了鞍点或者局部极小值,表现为梯度接近于0而且一般此时学习率已经调整到非常小的程度了。跳出驻点的思路就是说:此时可以重新把已经调整的很小的学习率数值放大,强行逼迫优化算法跳出此刻找到的鞍点或者极值点,在上述示意的损失函数曲面上等于当前在某个山谷的谷底(或者某个平台上,即鞍点),此时迈出一大步跳到比较远的地方,然后继续常规的优化过程,就是不断调小学习率,期望找到附近的另外一个山谷谷底,如此往复若干次,然后找出这些谷底哪个谷底更深,然后采取这个最深的谷底对应的参数作为模型参数。


    所以这里的关键是在优化进行不下去的时候重新调大学习率跳出当前的山谷,根据起跳点的不同,可以有两种不同的方法:


      方法1:随机漫步起跳法

    这个方法比较简单,就是说当前优化走到哪个山谷谷底,就从当前位置起跳,直接调大学习率即可,因为这样类似于走到哪里走不下去就从哪里开始跳到附近的远处,走到哪算哪,所以称之为“随机漫步起跳法”。


      方法2:历史最优起跳法

    上面说过,在这种优化方法过程中,会不断进入某个谷底,那么假设在训练集上优化历史上已经经过K个谷底,在验证集上其对应参数的模型性能是可以知道的,分别为P1,P2……Pk。随机漫步起跳法就是从Pk作为起跳点进行起跳。而历史最优起跳法则从P1,P2……Pk中性能最优的那个山谷谷底开始起跳,假设Pi=Max(P1,P2……Pk),那么从第i个山谷开始起跳,如果跳出后发现第j个山谷的性能Pj要优于Pi,则后面从Pj的位置开始起跳,当从Pi跳过K次都没有发现比Pi性能更好的山谷,此时可以终止优化算法。这个有点像是Pi性能的“一山望着一山高”,然后不断从历史最高峰跳到更高的高山上,可以保证性能应该越来越好。从道理上讲感觉这种历史最优起跳法要好于随机漫步起跳法,但是实验结果表明两种不同的起跳方法性能是差不多的,有时候这个稍好点,有时候那个稍好点,没有哪个方法确定性的比另外一个好。


    二.相关的实验及其结果


    为了检验上述跳出驻点的思路是否有效,我们选择了若干深层CNN模型来测试,这里列出其中某个模型的结果,这个模型结构可能看上去有点奇怪,前面说过这个思路是我们去年在做网络最优结构自动探索任务中产生的,所以这两个比较奇怪的结构是从自动生成的网络结构里面随机选出的,性能比较好。分类数据是cifar10图片分类数据,选择带动量的SGD算法作为基准方法来进行相关实验的验证。


         2.1模型结构及参数配置

    CNN模型结构如下图所示:


    其模型结构为:model structure:('init-29', [('conv', [('nb_filter', 16)]), ('conv',[('nb_filter', 256)]), ('maxp', [('pool_size', 2)]), ('conv', [('nb_filter',256)]), ('conv', [('nb_filter', 256)]), ('conv', [('nb_filter', 256)]),('maxp', [('pool_size', 2)]), ('conv', [('nb_filter', 256)])])")


    后面部分接上fc(512,relu)+dropout(0.2)+fc(10,softmax)的分类层;


    学习算法设置的参数为: 选择使用Nesterov的带动量的随机梯度下降法,初始学习率为0.01,衰减因子为1e-6,momentum因子为0.9。

     

    2.2 Nesterov+动量SGD(固定学习率)实验效果

    在上节所述的模型设置下,从结果看,最好的测试精度是0.9060,模型的误差和精度都变化地比较平滑


    2.3学习率自适应减少进行训练的实验结果

    其它配置类似于2.1.2所述,唯一的不同在于:使用常规策略对学习率进行缩小的改变。改变学习率的策略是如果误差10次不变化,那么学习率缩小为原来的0.1倍。 总迭代次数设置为100次。 训练误差、训练准确率、测试误差、测试准确率的变化曲线为:


     

    从训练结果看,第47次学习率变化,性能上有一次大幅度的提高:


    从训练和测试的结果看,测试精度最高的是在第80轮,达到了0.9144。此时的训练误差0.3921,训练准确率0.9926,测试误差0.3957,测试准确率0.9144。


    由此可见,增加学习率的动态改变,在学习速度和质量上要比原始的学习率恒定的训练方法要好,这里性能上增加近一个百分点。

     

    2.4随机漫步起跳法实验结果

     

    从上面的两组实验结果过程的最后阶段可以看出,误差和准确率变化较小,很可能优化过程陷入一个局部最优值,可能是一个极值,也可能是一个鞍点,不管怎样称呼它,它就是一个驻点。(从结果看,应该是一个局部最优值)。在其它实验条件与2.1.3节所述相同情况下,采取上文所述的随机漫步起跳法强迫优化算法跳出当前局部最优,即在学习率比较低的时候,进行扩大学习率的方式(直接设置为0.1,后面再常规的逐步减少),让其在误差曲面的区域中去寻找更多的驻点。 具体做法是:如果当前发现经过三次学习率变小没有带来准确率的提升,那么在当前训练的点,将学习率恢复到初始的学习率即0.01。 迭代的效果如下训练误差、训练准确率、测试误差、测试准确率的变化曲线如下图为:  


    从训练过程看,从当前迭代的当前位置出发,采用学习率缩小和扩大相结合的方式,目前最高能在第153轮的准确率提高到0.9236,比之前最好的0.9144有一点提升,说明这样的学习率变大策略也是一种可行的思路。


    另外,从图中看,由于我们变大学习率许多次,每次学习率变大都会呈现一个波动形式,故结果表现出锯齿状的波动趋势。


    2.5历史最优起跳法实验结果

       我们做了类似的实验,发现历史最优起跳法与随机漫步起跳法性能相当,有时这个好点有时那个好点,没有发现性能明显差异。

     

    我们随机抽了其它几个最优网络结构自动探索出的深度神经网络结构,表现与上述网络结构类似,而如果基础网络结构本身性能越差,这个差别越大。


    所以可以得出的结论是:优化过程中采用常规动态学习率调整效果要优于固定学习率的方法,而随机漫步起跳法效果略优于常规的动态学习率调整方法,但是历史最优起跳法相对常规动态学习率方法没有发现改善。

     

    相关文献

    1. SGDR: STOCHASTIC GRADIENT DESCENT RESTARTS .ICLR-2017
    2. SNAPSHOT ENSEMBLES : TRAIN 1,GET M FOR FREE

    展开全文
  • 遗传算法是基于达尔文进化论和孟德尔遗传学而发展形成,它是处理非线性模型参数估计的一类通用性较强的寻优方法。遗传算法寻求最优解的基本思想是:从代表问题的可能潜在解集的一个种群出发,此种群由基因编码的一定...
  • MatLab一维数组求解局部极值

    千次阅读 2016-08-04 11:58:02
    MatLab 函数 fingpeaks 对一维数组求解局部极大和极小值徐海蛟老师课堂教学举例说明。clc; clear;% 清屏清空变量figure('Color', 'w');% 背景:白色Data = [1 -2 3 -4 5 -6 7 8 5 4 1 2 -3 -1 -5 9 7 -6 5];plot...
  • 局部极值的小技巧

    2015-11-11 11:42:53
    nbr = np.zeros(len(vec)+1, dtype=bool) nbr[0] = True nbr[1:-1] = np.greater_equal(vec[1:], vec[:-1]) maxmask = (nbr[:-1] & ~nbr[1:])
  • close figure(2) syms x y f=2.*x.^4+y.^4-2.*x.^2-2.*y.^2+3; fx=diff(f,x); fy=diff(f,y); fxx=diff(f,x,2); fyy=diff(fy,y); fxy=diff(fx,y); fx fy fxx fyy fxy [a,b]=solve(fx,fy) ...tmp=2...
  • 从x1x_1x1​出发开始运用梯度下降法,得到局部最小m1m_1m1​ 随机选取x2x_2x2​点作为初始点 从x2x_2x2​出发开始运用梯度下降法,得到局部最小m2m_2m2​ 一直往下继续 最后选到了100个mmm值 选取100个中的最小即可...
  • 它是完全矢量化的,而且速度非常快。 设计用于 EMD/HHT 算法。
  • MatLab 函数 fingpeaks 对一维数组求解局部极大和极小值 徐海蛟老师课堂教学 举例说明。 clc; clear;% 清屏清空变量 figure('Color', 'w');% 背景:白色 Data = [1 -2 3 -4 5 -6 7 8 5 4 1 2 -3 -1 -5 9 7 -6 5]; ...
  • 真的结束于最优点吗?...我们知道,在局部最优点附近,各个维度的导数都接近0,而我们训练模型最常用的梯度下降法又是基于导数与步长的乘积去更新模型参数的,因此一旦陷入了局部最优点,...
  • 多层感知机解决了之前无法模拟异或逻辑的缺陷,同时更多的层数也让网络更能够刻画现实世界中的复杂情形。...但是随着神经网络层数的加深,优化函数越来越容易陷入局部最优解(即过拟合,在训练样本上有很好的拟...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 1,328
精华内容 531
关键字:

局部极值