Finding the local maxima/peaks and minima/valleys of histograms

Question

Ok, so I have a histogram (represented by an array of ints), and I'm looking for the best way to find local maxima and minima. Each histogram should have 3 peaks, one of them (the first one) probably much higher than the others.

I want to do several things:

1. Find the first "valley" following the first peak (in order to get rid of the first peak altogether in the picture)

2. Find the optimum "valley" value in between the remaining two peaks to separate the picture

   I already know how to do step 2 by implementing a variant of Otsu. But I'm struggling with step 1

3. In case the valley in between the two remaining peaks is not low enough, I'd like to give a warning.

Also, the image is quite clean with little noise to account for

What would be the brute-force algorithms to do steps 1 and 3? I could find a way to implement Otsu, but the brute-force is escaping me, math-wise. As it turns out, there is more documentation on doing methods like otsu, and less on simply finding peaks and valleys. I am not looking for anything more than whatever gets the job done (i.e. it's a temporary solution, just has to be implementable in a reasonable timeframe, until I can spend more time on it)

I am doing all this in c#

Any help on which steps to take would be appreciated! Thank you so much!

EDIT: some more data:

most histogram are likely to be like the first one, with the first peak representing background.

Arman

Could you give some sample data please? – Moluccas

Does the area around the peaks look like it's normal distributed? You could e.g. fit three independent normal distributions to your data. Then you can use standard deviation to decide on cut off points to identify your peaks and the valleys. – Ofris

What about using a k-means Algortihm with k=3 to obtain 3 different clusters? Each centroid should correspond to one of the peaks, if things go well. – Drag

@Arman have you tried considering the histogram as an image and constructing a convex hull and then proceeding to do the analysis on it ? I remember reading a paper that is about 30 years old now, which does something related to this method I'm describing very shortly. – Kindless

You could try to smooth your histograms first (by running average, or Gaussian convolution) and then apply a numerical derivative (i.e. take differences between neighboring values). Extrema can be detected, for example, by a sign change in the first derivative. – Southernly

@David Thanks, that sounds like something I could try. Unfortunately, I changed my design a bit, so histogram segmentation is not quite as important right now as are other things, but I'll try it out as soon as I can get to it and let you know – Arman

Answer 1

Use peakiness-test. It's a method to find all the possible peak between two local minima, and measure the peakiness based on a formula. If the peakiness higher than a threshold, the peak is accepted.

Source: UCF CV CAP5415 lecture 9 slides

Below is my code:

public static List<int> PeakinessTest(int[] histogram, double peakinessThres)
{
    int j=0;
    List<int> valleys = new List<int> ();

    //The start of the valley
    int vA = histogram[j];
    int P = vA;

    //The end of the valley
    int vB = 0;

    //The width of the valley, default width is 1
    int W = 1;

    //The sum of the pixels between vA and vB
    int N = 0;

    //The measure of the peaks peakiness
    double peakiness=0.0;

    int peak=0;
    bool l = false;

    try
    {
        while (j < 254)
        {

            l = false;
            vA = histogram[j];
            P = vA;
            W = 1;
            N = vA;

            int i = j + 1;

            //To find the peak
            while (P < histogram[i])
            {
                P = histogram[i];
                W++;
                N += histogram[i];
                i++;
            }


            //To find the border of the valley other side
            peak = i - 1;
            vB = histogram[i];
            N += histogram[i];
            i++;
            W++;

            l = true;
            while (vB >= histogram[i])
            {
                vB = histogram[i];
                W++;
                N += histogram[i];
                i++;
            }

                //Calculate peakiness
            peakiness = (1 - (double)((vA + vB) / (2.0 * P))) * (1 - ((double)N / (double)(W * P)));

            if (peakiness > peakinessThres & !valleys.Contains(j))
            {
                //peaks.Add(peak);                        
                valleys.Add(j);
                valleys.Add(i - 1);
            }

            j = i - 1;
        }
    }
    catch (Exception)
    {
        if (l)
        {
            vB = histogram[255];

            peakiness = (1 - (double)((vA + vB) / (2.0 * P))) * (1 - ((double)N / (double)(W * P)));

            if (peakiness > peakinessThres)
                valleys.Add(255);

                //peaks.Add(255);
            return valleys;
        }   
    }

        //if(!valleys.Contains(255))
        //    valleys.Add(255);

    return valleys;
}

Pentobarbital