斑点检测

     参考：http://blog.csdn.net/abcd1992719g/article/details/27071273

             http://www.cnblogs.com/ronny/p/3895883.html#3036526

            这两位博主都实现了这个算法，里面有很多细节，还是很有收获的。

           

#include "opencv2/core/core.hpp"
#include "opencv2/features2d/features2d.hpp"
#include "opencv2/highgui/highgui.hpp"
#include "opencv2/calib3d/calib3d.hpp"
#include "opencv2/nonfree/nonfree.hpp"
#include<opencv2/highgui/highgui.hpp>
#include<opencv2/imgproc/imgproc.hpp>
#include <iostream>
using namespace cv;
using namespace std;

/******************************Blob特征实现**********************************************/

/*****生成log sigma******/
vector<double> create_sigma(double start, double step, double end)
{
    vector<double> sigma;
    while (start <= end + 1e-8)
    {
        double s = exp(start);
        sigma.push_back(s);
        start += step;
    }

    return sigma;
}

/********生成log算子模板************/
Mat create_LOG(int size, double sigma)
{
    Mat H(size, size, CV_64F);
    int cx = (size - 1) / 2;
    int cy = (size - 1) / 2;
    double sum = 0;
    for (int i = 0; i < H.rows; i++)
    {
        for (int j = 0; j < H.cols; j++)
        {
            int nx = i - cx;
            int ny = j - cy;
            double s = -1 / (3.1415926 * sigma*sigma*sigma*sigma);
            s = s* (1 - (nx*nx + ny*ny) / (2 * sigma*sigma));
            s = s*exp(-(nx*nx + ny*ny) / (2 * sigma*sigma));
            sum += s;
            H.at<double>(i, j) = s;
        }
    }
    double mean = sum / (size*size);
    for (int i = 0; i < H.rows; i++)
    {
        for (int j = 0; j < H.cols; j++)
            H.at<double>(i, j) -= mean;
    }

    return H;
}

struct blob
{
    int x;
    int y;
    double sigma;//sigma用于计算半径r=1.5*sigma
    int intensity;//intensity用于当两个 blob重叠度高时取intensity大的（梯度大）
};

/*****生成LOG卷积后的scale_space******/
vector<Mat> creat_scale_space(Mat& im_in, vector<double>& sigma)
{
    vector<Mat> scale_space;
    for (int i = 0; i < sigma.size(); i++)
    {
        int n = ceil(sigma[i] * 3) * 2 + 1;
        Mat LOG = create_LOG(n, sigma[i]);
        Mat dst;
        //为了使得response一致，我们对做卷积的图片要先乘上sigma^2，这一步叫做scale normalize
        filter2D(im_in.mul(sigma[i] * sigma[i]), dst, -1, LOG, Point(-1, -1));
        scale_space.push_back(dst);
    }

    return scale_space;
}

/*******检测所有可能的blob**********/
vector<struct blob> detect_blobs(Mat& im_in, double thresh, vector<double> et, int padding = 10)
{
    vector<struct blob> blobs;
    Mat addborder, norm;
    /*The function copies the source image into the middle of the destination image.
    The areas to the left, to the right, 
    above and below the copied source image will be filled with extrapolated pixels*/

    copyMakeBorder(im_in, addborder, padding, padding, padding, padding, 
        BORDER_CONSTANT, Scalar(255));

    normalize(addborder, norm, 1, 0, NORM_MINMAX, CV_64F);

    vector<Mat> all_ims = creat_scale_space(norm, et);


    for (int i = 0; i < et.size(); i++)
    {
        Mat grayscale, mx;
        normalize(all_ims[i], grayscale, 0, 255, NORM_MINMAX, CV_8U, Mat());
        Mat structedElement(3, 3, CV_8U, Scalar(1));

        dilate(grayscale, mx, structedElement);//默认就是3*3，第三个参数可省略

        grayscale = (grayscale == mx)&(all_ims[i]>thresh);//取大于threshold并且是周围最大  

        for (int j = 0; j < norm.rows; j++)
        {
            for (int k = 0; k < norm.cols; k++)
            {
                if (grayscale.at<double>(j, k)>0)//是局部最大值点
                {
                    struct blob b;
                    b.x = j - padding;
                    b.y = k - padding;
                    b.sigma = et[i];
                    b.intensity = all_ims[i].at<double>(j, k);
                    blobs.push_back(b);
                }
            }
        }

    }

    return blobs;
}

/*****************删除重叠度大的blob********************/
vector<struct blob> prune_blobs(vector<struct blob>blobs_in)
{
    vector<bool> keep(blobs_in.size(), true);

    for (int i = 0; i < blobs_in.size(); i++)
    {
        for (int j = 0; j < blobs_in.size(); j++)
        {
            if ((i == j) || (keep[i] == false) || (keep[j] == false))//??????????
                continue;

            struct blob blobi = blobs_in[i];
            struct blob blobj = blobs_in[j];

            int xi = blobi.x;
            int yi = blobi.y;
            double ri = blobi.sigma;

            int xj = blobj.x;
            int yj = blobj.y;
            double rj = blobj.sigma;
            //算重叠度
            double d = sqrt((xj - xi)*(xj - xi) + (yj - yi)*(yj - yi));
            double rirj = ri + rj;
            double overlap_percent = rirj / d;

            if (overlap_percent > 2.0)
            {
                if (blobi.intensity > blobj.intensity)
                {
                    keep[i] = true;
                    keep[j] = false;
                }
                else {
                    keep[i] = false;
                    keep[j] = true;
                }
            }
        }
    }

    vector<struct blob> blobs_out;
    for (int i = 0; i < blobs_in.size(); i++)
    if (keep[i])
        blobs_out.push_back(blobs_in[i]);

    return blobs_out;
}
int main()
{
    Mat src, gray;
    src = imread("1.jpg");
    cvtColor(src, gray, CV_RGB2GRAY);
    vector<double> sigma = create_sigma(1.0, 0.2, 3.0);
    vector<struct blob> blobs = detect_blobs(gray, 0.2, sigma);
    vector<struct blob> blobs_trimmed = prune_blobs(blobs);

    for (int i = 0; i < blobs_trimmed.size(); i++)
    {
        struct blob b = blobs_trimmed[i];
        circle(src, Point(b.y, b.x), 1.5*b.sigma, Scalar(0), 2, 8, 0);
    }

    namedWindow("result", 1);
    imshow("result", src);

    waitKey(0);
    return 0;
}