精华内容
下载资源
问答
  • 为求得类似仿真函数的黑箱函数优化问题的全部局部极值点,提出了一种基于适应值曲面分析的新算法。首先,对适应值距离相关系数(fitness distance correlation, FDC)进行了改进,探讨了改进FDC指标与适应值曲面崎岖度...
  • 本次实验利用了基于局部极值的分水岭算法来实现圆斑点的检测。在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

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

    2021-05-26 06:10:43
    通常考察一个算法的性能通常...一般的局部搜索算法一旦陷入局部极值点,算法就在该点处结束,这时得到的可能是一个糟糕的结果。解决的方法就,目标函数差的点,被选中的概率小。考虑归一化问题,使得邻域内所有点被...

    通常考察一个算法的性能通常用局部搜索能力和全局收敛能力这两个指标。局部搜索是指能够无穷接近最优解的能力,而全局收敛能力是指找到全局最优解所在大致位置的能力。局部搜索能力和全局搜索能力,缺一不可。向最优解的导向,对于任何智能算法的性能都是很重要的。 

        局部最优问题(或叫局部峰值局部陷井):现实问题中,f在D上往往有多个局部的极值点。一般的局部搜索算法一旦陷入局部极值点,算法就在该点处结束,这时得到的可能是一个糟糕的结果。解决的方法就,目标函数差的点,被选中的概率小。考虑归一化问题,使得邻域内所有点被选中的概率和为1。总的来说,局部搜索的方法,就是依赖于对解空间进行按邻域搜索。

        起始点问题:一般的局部搜索算法是否能找到全局最优解,与初始点的位置有很大的依赖关系。解决的方法就是随机生成一些初始点,从每个初始点出发进行搜索,找到各自的最优解。再从这些最优解中选择一个最好的结果作为最终的结果。

        学习的重要性:
            1、直接用于无约束的实际问题;
            2、其基本思想和逻辑结构可以推广到约束问题;
            3、约束问题可以转化成无约束问题求解。

        方法分类:
            1、间接法:对简单问题,求解必要条件或充分条件;
            2、迭代算法:
                     零阶法:只需计算函数值 f(x) (直接法)
                     一阶法:需计算 ▽f(x)       (梯度法)
                     二阶法:需计算 ▽2f(x)      (梯度法)

        直接搜索法优点:计算不太复杂;易于实施与快速调试;所需的准备时间较少
    一、单纯形搜索法:

    1962年由Spendley, Hext和Himsworth提出,80年代得到证明。单纯形搜索法是一种无约束最优化的直接方法。单纯形法是求解非线性多元函数、无约束最小化问题的有效方法之一。在许多技术领域内,都取得了有效的成果。

    所谓的单纯形是指n维空间E^n中具有n+1个顶点的凸多面体。比如一维空间中的线段,二维空间中的三角形,三维空间中的四面体等,均为相应空间中的单纯形。单纯形搜索法与其它直接方法相比,基本思想有所不同,在这种方法中,给定维空间E^n中一个单纯形后,求出n+1个顶点上的函数值,确定出有最大函数值的点(称为最高点)和最小函数值的点(称为最低点),然后通过反射、扩展、压缩等方法(几种方法不一定同时使用)求出一个较好点,用它取代最高点,构成新的单纯形,或者通过向最低点收缩形成新的单纯形,用这样的方法逼近极小点。
        单纯形搜索法的基本思想是:先找出一个基本可行解,对它进行鉴别,看是否是最优解;若不是,则按照一定法则转换到另一改进的基本可行解,再鉴别;若仍不是,则再转换,按此重复进行。因基本可行解的个数有限,故经有限次转换必能得出问题的最优解。

    一般解题步骤可归纳如下:

    ①把线性规划问题的约束方程组表达成典范型方程组,找出基本可行解作为初始基本可行解。

    ②若基本可行解不存在,即约束条件有矛盾,则问题无解。

    ③若基本可行解存在,从初始基本可行解作为起点,根据最优性条件和可行性条件,引入非基变量取代某一基变量,找出目标函数值更优的另一基本可行解。

    ④按步骤3进行迭代,直到对应检验数满足最优性条件(这时目标函数值不能再改善),即得到问题的最优解。

    ⑤若迭代过程中发现问题的目标函数值无界,则终止迭代。 

    二、Powell共轭方向法(方向加速法)

    1964年由Powell提出,后经Zangwoll(1967年)和Brent(1973年)改进。该算法有效地利用了迭代过程中的历史信息,建立起能加速收敛的方向。

    理论基础:以二次对称函数为模型进行研究。

    在非线性目标函数中,最简单的是二次函数,故任何对一般函数有效的方法首先应对二次函数有效;在最优点附近,非线性函数可用一个二次函数作近似,故对二次函数使用良好的方法,通常对一般函数也有效;很多实际问题的目标函数是二次函数。

    设A是n×n阶对称正定矩阵,p(0), p(1)为两个n维向量,若 成立p(0)T A p(1) = 0,则称向量p(0)与p(1)为A共轭或A正交,称该两向量的方向为A共轭方向。 

     

    特点:Powell法本质上是以正定二次函数为背景,以共轭方向为基础的一种方法。不同于其他的直接法, Powell法有一套完整的理论体系,故其计算效率高于其他直接法。若每次迭代的前n 个搜索方向都线性无关时,则原始Powell法具有二次终止性 Powell法使用一维搜索,而不是跳跃的探测步。Powell法的搜索方向不一定为下降方向。

    不足:在原始Powell法中,必须保持每次迭代中前n个搜索方向线性无关,否则将永远得不到问题的最优解。

    为了避免出现搜索方向组线性相关的现象,Powell及其他人对原始Powell法进行修正:

     

    三、梯度法

        又名最速下降法,求解无约束多元函数极值的数值方法,早在1847年就已由柯西(Cauchy))提出。它是导出其他更为实用、更为有效的优化方法的理论基础。因此,梯度法是无约束优化方法中最基本的方法之一。该方法选取搜索方向Pκ的出发点是:怎样选取Pk可使ƒ(X)下降得最快?或者说使ƒ(Xκ+λΡκ)-ƒ(Χκ)<0且不等式左式的绝对值尽量大。 

    目标:求出平稳点(满足Ñf(x)=0的x * )。对给定的解,按负梯度方向搜索其邻域解,若邻域中有更优的解(根据给定的评估函数评定),则移动至该邻域解,继续寻找,直至找不到为止。此时,就找到了局部最优解。方法非常简单,但是缺陷也很明显,即陷入到局部最优解之后无法跳出。

     

     

    注意:最速下降法的搜索路径呈直角锯齿形,相邻两个搜索方向是正交的。

     

     

    1、优点:计算简单,需记忆的容量小;对初始点要求低,稳定性高;远离极小点时收敛快,常作为其它方法的第一步。

    2、缺点:收敛速度较慢(线性或不高于线性),原因是最速下降方向只有在该点附近有意义。直线搜索可能会产生一些问题,可能会“之字型”地下降。尤其当目标函数等值面是很扁的椭圆、椭球或类似图形时,收敛更慢。

    四、Newton法

    由最速下降法可知,从全局角度来看,负梯度方向一般不是一个特别好的方向,有没有更好的方向?

    基本思想:利用目标函数f(x)在x(k)处的二阶Taylor展开式去近似目标函数,用二次函数的极小点去逼近目标函数的极小点。

       或  

     

    1、Newton法的最大优点是:Newton法是局部收敛的,当初始点选得合适时收敛很快,具有二阶收敛速度,是目前算法中最快的一种。 对初始点要求高,一般要求初始点离极小点较近,否则不收敛。有时即使是收敛的,但因初始点离极大点或鞍点较近,会收敛于极大点或鞍点。

    2、Newton法的搜索方向-H (x)-1 g(x),称为Newton方向,是一个好方向,对二次函数此方向直指平稳点。如果H(x(k))是正定的,则H(x(k))-1必存在,从而算法是可行的,并且保证求得的平稳点是极小点。但在迭代过程中要求H(x(k))是正定的这一条件不一定能保证,因而它不一定是下降方向。一般在极小点附近的Hesse矩阵容易为正定的。所以基本Newton法在极小点附近才比较有效。

    五、共轭梯度法

    共轭梯度法最早是又Hestenes和Stiefle(1952)提出来的,用于解正定系数矩阵的线性方程组,在这个基础上,Fletcher和Reeves (1964)首先提出了解非线性最优化问题的共轭梯度法。由于共轭梯度法不需要矩阵存储,且有较快的收敛速度和二次终止性等优点,现在共轭梯度法已经广泛地应用与实际问题中。

    共轭梯度法(Conjugate Gradient)是介于最速下降法与牛顿法之间的一个方法,它仅需利用一阶导数信息,但克服了最速下降法收敛慢的缺点,又避免了牛顿法需要存储和计算Hesse矩阵并求逆的缺点,共轭梯度法不仅是解决大型线性方程组最有用的方法之一,也是解大型非线性最优化最有效的算法之一。 

    共轭梯度法是一个典型的共轭方向法,它的每一个搜索方向是互相共轭的,而这些搜索方向d仅仅是负梯度方向与上一次迭代的搜索方向的组合。在各种优化算法中,共轭梯度法是非常重要的一种。其优点是所需存储量小,具有步收敛性,稳定性高,而且不需要任何外来参数。

     

     

    1、共轭梯度法不需要预先估计任何参数就可以计算,每次迭代所需的计算,主要是向量之间的运算,便于并行化。程序简单,对大规模问题很有吸引力。对一般函数为超线性收敛。

    2、共轭梯度法的收敛速度比最速下降法要快得多,同时也避免了要求海塞矩阵的计算。收敛速度依赖于一维搜索的精确性及Hesse矩阵的特征值的分布。

     

     

     

     

     

        局部领域搜索是基于贪婪思想持续地在当前解的领域中进行搜索,虽然算法通用易实现,且容易理解,但其搜索性能完全依赖于领域结构和初解,尤其窥陷入局部极小而无法保证全局优化性。针对局部领域搜索,为了实现全局优化,可尝试的途径有:以可控性概率接受劣解来逃逸局部极小,如模拟退火算法;确定性的局部极小突跳策略,如禁忌策略;扩大领域搜索结构,如TSP的2opt扩展到k-opt;多点并行搜索,如进化计算;变结构领域搜索( Mladenovic et al,1997)等等。

    一、爬山法(Hill-Climbing)

    爬山算法是一种简单的贪心搜索算法,该算法每次从当前解的临近解空间中选择一个最优解作为当前解,直到达到一个局部最优解。爬山算法实现很简单,其主要缺点是会陷入局部最优解,而不一定能搜索到全局最优解。如下图所示:假设C点为当前解,爬山算法搜索到A点这个局部最优解就会停止搜索,因为在A点无论向那个方向小幅度移动都不能得到更优的解。

     

    爬山法是向值增加的方向持续移动到简单循环过程,算法在到达一个“峰顶”时终止,此时相邻状态中没有比该“峰顶”更高的值。爬山法不维护搜索树,当前节点只需要记录当前状态及其目标函数值;爬山法不会前瞻与当前状态不直接相邻的状态的值——“就像健忘的人在大雾中试图登顶珠峰一样”。爬山法从来不会“下山”,只会向值比当前节点好的方向搜索,因而肯定不完备,很容易停留在局部极值上。
    function HillClimbing(problem) return 一个局部最优状态
        输入:problem
        局部变量:current, 一个节点
                  neighbor,一个节点
        current= MakeNode(Initial-State(problem));
        loop do
            neighbor= a highest-valued successor of current ;
            if VALUE[neighbor] <= VALUE[current] then return STATE[current];
            current= neighbor ;
        爬山法又称贪婪局部搜索,只是选择相邻状态中最好的一个。尽管贪婪是七宗罪之一,但是贪婪算法往往能够获得很好的效果。当然,爬山法会遇到以下问题:

    (1)局部极值

    (2)山脊:造成一系列的局部极值

    (3)高原:平坦的局部极值区域——解决办法:继续侧向移动
        随机爬山法:在上山移动中,随机选择下一步,选择的概率随着上山移动到陡峭程度而变化。

    首选爬山法:随机地生成后继节点直到生成一个优于当前节点的后继。

    随机重新开始的爬山法:“如果一开始没有成功,那么尝试,继续尝试”算法通过随机生成的初始状态来进行一系列的爬山法搜索,找到目标时停止搜索。该算法以概率1接近于完备:因为算法最终会生成一个目标状态作为初始状态。如果每次爬山搜索成功的概率为p,则需要重新开始搜索的期望次数为 1/p。

    二、模拟退火算法(Simulated Annealing)

    “模拟退火”算法是源于对热力学中退火过程的模拟,在某一给定初温下,通过缓慢下降温度参数,使算法能够在多项式时间内给出一个近似最优解。

        爬山法是完完全全的贪心法,每次都鼠目寸光的选择一个当前最优解,因此只能搜索到局部的最优值。模拟退火其实也是一种贪心算法,但是它的搜索过程引入了随机因素。模拟退火算法与爬山法类似,但是它没有选择最佳的移动,而是选择随机的移动。如果该移动使情况得到改善,那么接受该移动;否则,算法以某个概率接受该移动。因此有可能会跳出这个局部的最优解,达到全局的最优解。以上图为例,模拟退火算法在搜索到局部最优解A后,会以一定的概率接受到E的移动。也许经过几次这样的不是局部最优的移动后会到达D点,于是就跳出了局部最大值A。

    模拟退火的原理:我们将热力学的理论套用到统计学上,将搜寻空间内每一点想像成空气内的分子;分子的能量,就是它本身的动能;而搜寻空间内的每一点,也像空气分子一样带有“能量”,以表示该点对命题的合适程度。算法先以搜寻空间内一个任意点作起始:每一步先选择一个“邻居”,然后再计算从现有位置到达“邻居”的概率。可以证明,模拟退火算法所得解依概率收敛到全局最优解。

    当陷入局部最优之后,模拟退火算法要求概率根据算法进行过程中,逐步降低(逐渐降低才能趋向稳定)其跳出局部最优的概率,使其越来越趋于稳定。这一特性增加也同样存在缺陷,即对于全局最优解,也存在一定概率跳出,从而使得求解过程不稳定。

    随着邻域的范围的增大,跳出局部极小区域,最终进入全局极小区域的概率越来越大,但是代价是总的迭代次数增加。但是随着邻域范围的增大,会出现所谓的在极值附近来回”振荡”而无法落入极值点的现象。可以推测,随着邻域范围的进一步增大及其随机特性,模拟退火算法将退化到随机寻找极值并进行记录极值的算法。当然这种算法也存在问题,即对于某些情况下,也不易达到全局最优。例如,解空间中仅有两个局部最优,其中一个是全局最优,那么模拟退火似乎并不一定总能进入到全局最优解当中。

    Create initial solution S

    repeat

    for i=1 to iteration-length do

    Generate a random transition from S to Si

    If ( C(S) <= C(Si) ) then

    S=Si

    else if( exp(C(S)-C(Si))/kt > random[0,1) ) then

    S=Si

    Reduce Temperature t

    until ( no change in C(S) )

    C(S): Cost or Loss function of Solution S.

    1、在其他参数合适的情况下,初始解的位置与邻域结构是2个重要因素。

    2、初始解的位置与邻域结构的设计之间对于找到最优解的问题而言存在依赖关系。随着邻域的范围的增大,跳出局部极小区域,最终进入全局极小区域的概率越来越大,同时也可能会产生振荡,无法落入极值区域。

    3、模拟退火本质上也是一种 暴力搜索,只不过是利用随机的性质和局部极值区域连续(也取决于求解问题中邻域的结构)的特性避免了大量计算,进而在较短时间内从某种概率上逼近局部最优值和全局最优值。

    三、禁忌搜索(Tabu Search或Taboo Search,TS)

    禁忌搜索(Tabu Search或Taboo Search,简称TS)的思想最早由Glover(1986)提出,它是对局部领域搜索的一种扩展,是一种全局逐步寻优算法,是对人类智力过程的一种模拟。TS算法通过引入一个灵活的存储结构和相应的禁忌准则来避免迂回搜索,并通过藐视准则来赦免一些被禁忌的优良状态,进而保证多样化的有效探索以最终实现全局优化。相对于模拟退火和遗传算法,TS是又一种搜索特点不同的meta-heuristic算法。

    禁忌搜索是人工智能的一种体现,是局部领域搜索的一种扩展。禁忌搜索最重要的思想是标记对应已搜索的局部最优解的一些对象,并在进一步的迭代搜索中尽量避开这些对象(而不是绝对禁止循环),从而保证对不同的有效搜索途径的探索。禁忌搜索涉及到领域(neighborhood)、禁忌表(tabu list)、禁忌长度(tabu 1ength)、候选解(candidate)、藐视准则(candidate)等概念。

    简单TS算法的基本思想描述如下:

    (1)给定算法参数,随机产生初始解x,置禁忌表为空。

    (2)判断算法终止条件是否满足?若是,则结束算法并输出优化结果;否则,继续以下步骤。

    (3)利用当前解的邻域函数产生其所有(或若干)邻域解,并从中确定若干候选解。

    (4)对候选解判断特赦准则是否满足?若成立,则用满足特赦准则的最佳状态y替代x成为新的当前解,即x=y,并用与y对应的禁忌对象替换最早进入禁忌表的禁忌对象,同时用y替换“best so far”状态,然后转步骤6;否则,继续以下步骤。

    (5)判断候选解对应的各对象的禁忌属性,选择候选解集中非禁忌对象对应的最佳状态为新的当前解,同时用与之对应的禁忌对象替换最早进入禁忌表的禁忌对象元素。

    (6)转步骤(2)。

     

     

    局部搜索,模拟退火,遗传算法,禁忌搜索的形象比喻:为了找出地球上最高的山,一群有志气的兔子们开始想办法。
        1、兔子朝着比现在高的地方跳去。他们找到了不远处的最高山峰。但是这座山不一定是珠穆朗玛峰。这就是局部搜索,它不能保证局部最优值就是全局最优值。
        2、兔子喝醉了。他随机地跳了很长时间。这期间,它可能走向高处,也可能踏入平地。但是,他渐渐清醒了并朝最高方向跳去。这就是模拟退火。
        3、兔子们吃了失忆药片,并被发射到太空,然后随机落到了地球上的某些地方。他们不知道自己的使命是什么。但是,如果你过几年就杀死一部分海拔低的兔子,多产的兔子们自己就会找到珠穆朗玛峰。这就是遗传算法。
        4、兔子们知道一个兔的力量是渺小的。他们互相转告着,哪里的山已经找过,并且找过的每一座山他们都留下一只兔子做记号。他们制定了下一步去哪里寻找的策略。这就是禁忌搜索。 

    展开全文
  • 皇后问题局部极值启发式搜索方法 背景说明:皇后问题是算法领域的著名问题,本篇主要是对皇后问题背景下的局部搜索问题的研究。

    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. 最近新开通了博客 会放一些学习过程中的心得体会(和一些作业或者报告)

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

    展开全文
  • 它是完全矢量化的,而且速度非常快。 设计用于 EMD/HHT 算法
  • 结果表明,Matlab最优化工具箱中的导数算法在多层膜局部优化设计上具有更好的局部极值搜索性能和收敛速度;非导数算法性能较差且收敛时间较长,但具有更多的搜索路径,较适用于设计初期开拓搜索方向。在多层膜反演中...
  • 由于变量的适应度最优与问题的目标函数最优无法达到一致,从而利用极值过程原则的局部搜索算法对TSP问题效果不好,而通过改变变量的适应度,使其与目标函数相关,就能够提高个体解的搜索能力。比较参数取不同值时...
  • sift算法详解及应用

    2015-05-07 09:34:07
    对sift算法的比较详尽的解释:局部极值算法、关键点方向检测,适合初学者。
  • 通用局部搜索算法之模拟退火[转]

    千次阅读 2010-05-26 08:31:00
    通用局部搜索算法之模拟退火[转]设施区位及算法 2009-08-28 20:44:05 阅读24 评论0 字号:大中小 爬山法从来不会“下山”,只会向值比当前节点好的方向搜索,因而肯定不完备,很容易停留在局部极值上;纯粹的随机...
  • 用C++实现了模拟退火算法求多元函数极值,可以避免陷入局部最优解。
  • 遗传算法之求取函数极值 1、前言 在智能控制(刘金琨)这本里面讲了遗传算法求取函数极值的方法,这里给出一些个人理解时的注释,顺带 求解了第10章的课后习题第二题。 遗传算法流程图如下: 2、原书案例 利用遗传...
  • 算法实现过程中,首先用挖掘算法挖掘出已知联配中的不良序列块,然后所有的不良序列块用极值遗传算法重新联配。当初始的序列是用渐进算法联配时,新的求精方法能调整早期的一些错误,充分结合渐进和迭代算法的优点。...
  • 注:由于编码选择二进制编码,所以只能在整数范围进行搜索,也就是说求解到的最优解一定是最优的整数解,如果选择一些映射方法可以将离散问题连续化,但这样就和进化算法本身无关了,所以本文只写了基本的遗传算法 ...
  • 内容:常见的局部优化和全局优化算法。 时间:2019.01.05 梯度下降,牛顿法,共轭梯度法(Python + Matlab),找到空间表面的极值点(3D) 温馨提示: 梯度下降算法速度慢,迭代次数大,最终结果是近似的。 牛顿法...
  • 做数据分析算法,使用MATLAB进行算法研究,使用C#进行工程实现比较合适,目前出现这样的情况,有一个数组,经过某种超分辨算法得到的数据点很稀疏,而且...1、基于局部极值算法从原始数据数组获取局部极值数组(极大...
  • 针对粒子对算法存在过早陷入局部最优导致精度不是很高的问题,建议了一种新的基于粒子对(PPO) 与极值优化(EO)混合算法。该算法利用PPO和EO的优点,借助K-means快速聚类的结果初始化其中一个粒子,并根据一定迭代次数在...
  • 针对非局部平均(NLM)方法对椒盐噪声图像滤波效果较差的问题。通过引入噪声检测结果扩展NLM方法去除图像中椒盐噪声。在噪声检测阶段,利用图像的两个极值Lmin和Lmax, 把图像像素点分为非噪声点和噪声点。在滤波阶段,...
  • 特里斯坦·乌塞尔图像极值查找器... 某些极值可能无法在单个点上进行数学定义,即鞍点下方的示例是像素双峰——这是一个图像离散化的特性,不是算法的问题。 图像边缘不能出现图像极值。 如果噪声是一个问题,您可以
  • 针对粒子对算法存在过早陷入局部最优导致精度不是很高的问题,建议了一种新的基于粒子对(PPO)与极值优化(EO)混合算法。该算法利用PPO和EO的优点,借助K-means快速聚类的结果初始化其中一个粒子,并根据一定迭代次数...
  • 通过老师对算法代码的简单介绍,我对于粒子寻优算法结构有了一个大致的了解 首先根据老师说的原则,c1从小变到大,c2从大变到小,权重w先从大的开始,惯性大的冲得快,后期再往小改。 初始值: 首先更改了c2的...
  • 针对4个benchmark函数对新算法做了测试,并与粒子群优化算法以及已有的几个算法进行了比较,结果表明该算法跳出局部极值点的能力较强、收敛速度更快、寻优精度较高;最后将新算法应用到焊接梁的优化设计问题中,仿真...
  • 为了应用SIFT对图像的描述能力及其鲁棒性,并提高效率,对SIFT的提取算法进行了修改,消除可以引起边缘响应的部分极值点,消除图像细节丰富的局部过邻近点,消除图像背景中的低对比度点,以降低算法复杂度。...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 901
精华内容 360
关键字:

局部极值算法