Detecting peaks in images

Question

I got a large set of infrared images of seeds, their sizes vary slightly. And I would like to find them (in thefastest possible way). Below i show zoomed in details of the images i process. After a first noise removal and blob filter this is what i have :

The bright white is just direct reflection of the IR lamp, white pixels never combine (stretch out) over multiple seeds.

To make it more clear i placed a letter on some individual seeds.

The problems i have:

A is a single seed (dirt on the seed) generates a slight dark line.
B the nearby X close to it is at its darkest intersection its still brighter as some other seeds (cannt change brightnes or remove if gray value is lower then a certain value.
C those are 3 seeds close to each other.
The smallest visible seeds above should not becomme even smaller.

I'm not making use of mathlab or openCV, as i work directly with locked image and memory data. i can acces pixel data by array or simple getpixel / putpixel commands. I wrote my own graphics library which is fast enough for live camera data, processing speed currently at around 13ms at around 25ms i enter stream processing lag

I wonder how to separate those 'cloudy' blobs better. I'm thinking to find local maxima over a certain pixels range..but that should see A as one seed, while on B find that B and X are not connected. So I'm not sure here, how such a local peek function or another function should like like. Although i code in C# i looked at other c++ functions as well like dilate etc but thats not it. I also wrote a function to check for slope degree (like if it was a mountain height image) but that couldnt devide areas B and C.

Ok well i made different slope detection code, now i dont look for some degree but just the tilting point over a small range, it works nice on X axis.. but essentially i think it should work on both X and Y here's the new result : It can resolve issue A and B !!!

However it wouldnt differentiate between seeds who are aligned in a vertical row, and and it causes small white noise (non connected lines). at places where is there is nearly nothing to detect. I'm not yet sure on how to do the same (combined) over Y axis to get the tops then erase stuff from a certain distance of the top.. (to seperate).

Using this code just showing the pixel operations of it.

 for (int y =  raw.Height;y>5; y--)
        {slopeUP = false;
            int[] peek = new int[raw.Width];
            for (int x = raw.Width; x> 7; x--)
            {
                int a = raw.GetPixelBleu(x, y);
                int b = raw.GetPixelBleu(x - 1, y);
                int c = raw.GetPixelBleu(x - 2, y);
                int d = raw.GetPixelBleu(x - 11, y);
                int f = raw.GetPixelBleu(x - 12, y);

                if ((f + d) > (a + b))slopeUP  = true;

                if ((f + d) < (a + b))
                {
                    if (slopeUP)
                    {
                        slopeUP = false;
                        peek[x - 6] = 10;
                        //raw.SetPixel(x, y, Color.GreenYellow);
                    }
                    else peek[x - 6] = 0;
                }

            }
            for (int x = raw.Width; x > 7; x--) 
            { if (peek[x-1] > 5) raw.SetPixel(x, y, Color.Lavender); }
     }

Sry to heckle. But could we have the image pre-processing (ie before noise removal and blob filter). TBH i think those are the steps where we are going to be able to help a lot more rather than just trying to manage the blobs you have currently. — Sneaky Polar Bear, Sep 01 '17 at 14:37
Wow sry, re-read your post, the bright white is direct reflection saturation. I thought that was some type of mask or result from a blob filter. — Sneaky Polar Bear, Sep 01 '17 at 14:41
Try hitting it with a pretty large and strong Gaussian. I think this may help to combine A and should give you stronger peaks overall. — Sneaky Polar Bear, Sep 01 '17 at 14:50
yes, and its not troubling me the direct reflection, since those cannot overlap (and are always inside seed area). the local peek filter works but only the X axis.. by setting 'peek[x-6] = c;' i can remove some useless low peeks .. but wont seperate vertical alligned seeds.. i'm close but not there yet.. — Peter, Sep 01 '17 at 14:50
lol ya, I have officially decided I am overcomplicating this. A simple peaks with some light clustering afterwards is all this needs. — Sneaky Polar Bear, Sep 01 '17 at 18:14

S. Huber · Answer 1 · 2018-05-23T17:59:00.110

1

In this SO answer to a similar question I applied persistent homology to find peaks in an image. I took your image, scaled it down to 50%, applied an Guassian blur with radius 20 (in Gimp), and applied the methods described in the other article (click to enlarge):

I only show peaks with persistence (see the other SO answer) of at least 20. The persistence diagram is shown here:

Persistence diagram

The 20-th peak would be the little peak on the top-left corner with a persistence of about 9. By applying a stronger Gaussian filter the image gets more diffuse and the peaks will be come more prominent.

Python code can be found here.

edited May 23 '18 at 17:59

answered Nov 08 '17 at 22:50

S. Huber

933
13
25

Thanks i'm aware of the article my code works faster, as complete blur is an expensive operation, for Y axis the code can be slightly optimized as X areas are allready known a lot pixels can be skipped; i needed speed as the code had to work against a 30fps camera. – user613326 Nov 09 '17 at 23:06

Sneaky Polar Bear · Answer 2 · 2017-09-04T21:11:43.990

So, as far as speed, I am only going off of the image you posted here... on which everything runs blazing fast because it is tiny. Note that I padded the image after binarizing and never un-padded, so you will want to either un-pad or shift your results accordingly. You may not even want to pad, but it allows detection of cut off seeds.

Overview of pipeline: removeSaturation>>gaussian blur>>binarize>>padd>>distanceTransform>>peaks>>clustering

That being said here is my code and results:

void drawText(Mat & image);
void onMouse(int event, int x, int y, int, void*);
Mat bruteForceLocalMax(Mat srcImage, int searchRad);
void zoomPixelImage(Mat sourceImage, int multFactor, string name, bool mouseCallback);
Mat mergeLocalPeaks(Mat srcImage, int mergeRadius);
Mat image;
bool debugDisplays = false;


int main()
{
    cout << "Built with OpenCV " << CV_VERSION << endl;
    TimeStamp precisionClock = TimeStamp();

    image = imread("../Raw_Images/Seeds1.png",0);


    if (image.empty()) { cout << "failed to load image"<<endl; }
    else
    {
        zoomPixelImage(image, 5, "raw data", false);
        precisionClock.labeledlapStamp("image read", true);

        //find max value in image that is not oversaturated
        int maxVal = 0;
        for (int x = 0; x < image.rows; x++)
        {
            for (int y = 0; y < image.cols; y++)
            {
                int val = image.at<uchar>(x, y);
                if (val >maxVal && val !=255)
                {
                    maxVal = val;
                }
            }
        }

        //get rid of oversaturation regions (as they throw off processing)
        image.setTo(maxVal, image == 255);

        if (debugDisplays)
        {zoomPixelImage(image, 5, "unsaturated data", false);}
        precisionClock.labeledlapStamp("Unsaturate Data", true);

        Mat gaussianBlurred = Mat();
        GaussianBlur(image, gaussianBlurred, Size(9, 9), 10, 0);
        if (debugDisplays)
        {zoomPixelImage(gaussianBlurred, 5, "blurred data", false);}
        precisionClock.labeledlapStamp("Gaussian", true);

        Mat binarized = Mat();
        threshold(gaussianBlurred, binarized, 50, 255, THRESH_BINARY);
        if (debugDisplays)
        {zoomPixelImage(binarized, 5, "binarized data", false);}
        precisionClock.labeledlapStamp("binarized", true);

        //pad edges (may or may not be neccesary depending on setup)
        Mat paddedImage = Mat();
        copyMakeBorder(binarized, paddedImage, 1, 1, 1, 1, BORDER_CONSTANT, 0);
        if (debugDisplays)
        {zoomPixelImage(paddedImage, 5, "padded data", false);}
        precisionClock.labeledlapStamp("add padding", true);

        Mat distTrans =  Mat();
        distanceTransform(paddedImage, distTrans, CV_DIST_L1,3,CV_8U);
        if (debugDisplays)
        {zoomPixelImage(distTrans, 5, "distanceTransform", true);}
        precisionClock.labeledlapStamp("distTransform", true);

        Mat peaks = Mat();
        peaks = bruteForceLocalMax(distTrans,10);
        if (debugDisplays)
        {zoomPixelImage(peaks, 5, "peaks", false);}
        precisionClock.labeledlapStamp("peaks", true);

        //attempt to cluster any colocated peaks and find the best clustering count
        Mat mergedPeaks = Mat();
        mergedPeaks = mergeLocalPeaks(peaks, 5);
        if (debugDisplays)
        {zoomPixelImage(mergedPeaks, 5, "peaks final", false);}
        precisionClock.labeledlapStamp("final peaks", true);

        precisionClock.fullStamp(false);
        waitKey(0);
    }
}

void drawText(Mat & image)
{
    putText(image, "Hello OpenCV",
            Point(20, 50),
            FONT_HERSHEY_COMPLEX, 1, // font face and scale
            Scalar(255, 255, 255), // white
            1, LINE_AA); // line thickness and type
}

void onMouse(int event, int x, int y, int, void*)
{
    if (event != CV_EVENT_LBUTTONDOWN)
        return;

    Point pt = Point(x, y);
    std::cout << "x=" << pt.x << "\t y=" << pt.y << "\t value=" << int(image.at<uchar>(y,x)) << "\n";

}

void zoomPixelImage(Mat sourceImage, int multFactor, string name, bool normalized)
{
    Mat zoomed;// = Mat::zeros(sourceImage.rows*multFactor, sourceImage.cols*multFactor, CV_8U);
    resize(sourceImage, zoomed, Size(sourceImage.cols*multFactor, sourceImage.rows*multFactor), sourceImage.cols*multFactor, sourceImage.rows*multFactor, INTER_NEAREST);
    if (normalized) { normalize(zoomed, zoomed, 0, 255, NORM_MINMAX); }
    namedWindow(name);
    imshow(name, zoomed);
}

Mat bruteForceLocalMax(Mat srcImage, int searchRad)
{
    Mat outputArray = Mat::zeros(srcImage.rows, srcImage.cols, CV_8U);
    //global search top
    for (int x = 0; x < srcImage.rows - 1; x++)
    {
        for (int y = 0; y < srcImage.cols - 1; y++)
        {
            bool peak = true;
            float centerVal = srcImage.at<uchar>(x, y);
            if (centerVal == 0) { continue; }
            //local search top
            for (int a = -searchRad; a <= searchRad; a++)
            {
                for (int b = -searchRad; b <= searchRad; b++)
                {
                    if (x + a<0 || x + a>srcImage.rows - 1 || y + b < 0 || y + b>srcImage.cols - 1) { continue; }
                    if (srcImage.at<uchar>(x + a, y + b) > centerVal)
                    {
                        peak = false;
                    }
                    if (peak == false) { break; }
                }
                if (peak == false) { break; }
            }
            if (peak)
            {
                outputArray.at<uchar>(x, y) = 255;
            }

        }
    }
    return outputArray;
}

Mat mergeLocalPeaks(Mat srcImage, int mergeRadius)
{
    Mat outputArray = Mat::zeros(srcImage.rows, srcImage.cols, CV_8U);
    //global search top
    for (int x = 0; x < srcImage.rows - 1; x++)
    {
        for (int y = 0; y < srcImage.cols - 1; y++)
        {
            float centerVal = srcImage.at<uchar>(x, y);
            if (centerVal == 0) { continue; }

            int aveX = x;
            int aveY = y;

            int xCenter = -1;
            int yCenter = -1;

            while (aveX != xCenter || aveY != yCenter)
            {
                xCenter = aveX;
                yCenter = aveY;
                aveX = 0;
                aveY = 0;
                int peakCount = 0;
                //local search top
                for (int a = -mergeRadius; a <= mergeRadius; a++)
                {
                    for (int b = -mergeRadius; b <= mergeRadius; b++)
                    {
                        if (xCenter + a<0 || xCenter + a>srcImage.rows - 1 || yCenter + b < 0 || yCenter + b>srcImage.cols - 1) { continue; }
                        if (srcImage.at<uchar>(xCenter + a, yCenter + b) > 0)
                        {
                            aveX += (xCenter + a);
                            aveY += (yCenter + b);
                            peakCount += 1;
                        }
                    }
                }
                double dCentX = ((double)aveX / (double)peakCount);
                double dCentY = ((double)aveY / (double)peakCount);
                aveX = floor(dCentX);
                aveY = floor(dCentY);
            }
            outputArray.at<uchar>(xCenter, yCenter) = 255;
        }
    }
    return outputArray;
}

speed:

debug images:

results:

Hope this helps! Cheers!

I tried erode but isnt working for me, which has mainly todo that sometimes they are close with a rather bright dark line in between (the B problem), and removing the darker pixels ... but i updated local maxima code and will add it to the question. — Peter, Sep 01 '17 at 14:20
I think that "B" can only be solved by local maxima detection, while "A" demands that such a detection should be over some minimal length. The given image is just one of a few test cases that seamed hard to solve. Also I got versions where some are just barley visible (very dark), still local maxima did detect them without altering much of the graphics. Also that math is fast (simple adding and comparing).. but it needs to be combined with a similair vertical version.. then i think on the cross of horizontal and vertical detection i will get them all. — user613326, Sep 01 '17 at 22:11
for local maxima combined over X and Y axis my thoughts are on how to find the cross marks, (rewriting to y is easy part).. i think Peek should be a 2D array then so a Y pass can check if the X pass had allready written something below it. (its in mind but no access to the code and data till monday) — user613326, Sep 01 '17 at 22:17

Detecting peaks in images

2 Answers2