Using histogram information from C++ (Otsu threshold)

[kyron]

How would one access the histogram information of an image under C++. I have to implement Otsu' s threshold and would like to IM's infrastructure use as much as possible (unfortunate that it's not already implemented). This algorithm, although not the 'ultimate best' (check out Fred's ImageMagick Scripts ...does an excellent job at comparing thresholds), it is one that is very useful and referenced in research (am currently working on a signature recognition infrastructure).

... And no, a command line version does not suffice here as we need to do some real-time processing , so don't just point me to the post entitled Otsu' Thresholding Method to Binarize an image and Fred's ImageMagick Scripts

[magick]

You can access any MagickCore API method in the MagickCore namespace. For example, MagickCore::GetImageHistogram(). You could also generate your own histogram by iterating over the pixels with the Magick++ pixels class (see http://www.imagemagick.org/Magick++/Pixels.html).

[kyron]

Well, my obvious intent is not to write what is already there (histogram). Now, to use MagickCore methods (C functions I should say), I need to convert ‘Magick::Image*’ to ‘const MagickCore::Image*’ since that is what is expected by all calls to MagickCore functions. Have any suggestions? Write to blob than load from blob..sounds like excessive overhead, is there a simpler way? At this point, I must admit that writing my own GetHistogram function sounds more attractive, even if I am duplicating some code and adding memory management overhead...

Thanks.

[mjwvb]

It's a topic from 7 years ago so you probably don't need it anymore. Anyway, I needed it for a project of mine and couldn't find a working implementation anywhere, so I had to make it myself. If anybody still needs an Otsu thresholding implementation for magick++ / c++, I want to share the code with you. It's based on a paper by Juan Pablo Balarini, Sergio Nesmachnow (2016, http://www.ipol.im/pub/art/2016/158/), so all credits to them.

It runs on the latest imagemagick v6, so I'm not sure if it also works for v7. A few mutations regarding the implementation of Magick::Color (.redQuantum()) and Magick::PixelPacket will do.

Edit:Optimized code to make it ~3 times faster (reducing the total number of loops by replacing the Magick::colorHistogram() by a manual histogram creation loop)

Code: Select all

// Requires following includes
#include <vector>
#include <cmath>
#include <utility>
#include <Magick++.h>

/**
* Calculate the otsu threshold of given image.
*
* This function will automatically convert the input image to to grayscale when
* it is not provided as such.
*
* Returns the calculated threshold between 0 and QuantumRange
*
* Based on paper "A C++ Implementation of Otsu’s Image Segmentation Method"
* by Juan Pablo Balarini, Sergio Nesmachnow (2016)
* http://www.ipol.im/pub/art/2016/158/
*/
int otsuThreshold(const Magick::Image &image)
{
    int n = QuantumRange;
    long totalPixels = image.columns() * image.rows();

    // Prepare histogram image
    Magick::Image histImage(image);
    if (histImage.colorSpace() != Magick::GRAYColorspace) {
        histImage.colorSpace(Magick::GRAYColorspace);
    }

    // Find number of occurrences per color
    // The histogram should stretch all bins
    std::vector<long> histogram(n+1, 0);

    Magick::Pixels view(histImage);
    const Magick::PixelPacket *px = view.getConst(0, 0, image.columns(), image.rows());

    for (long i = 0; i < totalPixels; i++) {
        histogram[std::round(px->red)]++;
        *px++;
    }

    // Compute threshold
    // Init variables
    int threshold = 0;
    double sum = 0;
    double sumB = 0;
    long q1 = 0;
    long q2 = 0;
    double varMax = 0;

    // Auxiliary value for computing m2
    for (int i = 0; i <= n; i++){
        sum += (double) i * histogram[i];
    }

    for (int i = 0; i <= n; i++) {
        // Update q1
        q1 += histogram[i];
        if (q1 == 0) {
            continue;
        }

        // Update q2
        q2 = totalPixels - q1;
        if (q2 == 0) {
            break;
        }

        // Update m1 and m2
        sumB += (double) i * histogram[i];
        double m1 = sumB / q1;
        double m2 = (sum - sumB) / q2;

        // Update the between class variance
        double varBetween = (double) q1 * (double) q2 * (m1 - m2) * (m1 - m2);

        // Update the threshold if necessary
        if (varBetween > varMax) {
            varMax = varBetween;
            threshold = i;
        }
    }

    return threshold;
}

[snibgo]

Thanks for sharing your code. To help search engines find it, I've edited the topic title to include "Otsu threshold".

[fmw42]

For other interested parties, I have bash ImageMagick shell scripts that do otsu thresholding and several other automatic thresholds. See
http://www.fmwconcepts.com/imagemagick/ ... /index.php
http://www.fmwconcepts.com/imagemagick/index.php (see thresholding section of my scripts table)

However, it would be nice to have the otsu threshold and some of my others directly coded into ImageMagick.

[mjwvb]

I agree, it would be nice to have some of the common auto thresholding techniques in the ImageMagick core. Imo there are many situations wherein auto selection is better than the manual (subjective) selection of a threshold level by humans, or you simply want to provide a default threshold level. In fact, in certain cases you want to completely avoid human interaction at all. It's pretty simple code, but also very useful, so +1 from my side!

Liked your threshold comparison study btw, and its amazing conclusion . It helped me in selecting a proper auto thresholding technique. What I didn't understand was that Otsu, K-Means and Isodata are completely identical for the given inputs. I'm curious for which inputs they would differ...

(PS, I changed the code in my original post to make it a lot faster, by reducing the number of histogram creation loops)

[fmw42]

What I didn't understand was that Otsu, K-Means and Isodata are completely identical for the given inputs. I'm curious for which inputs they would differ...

The threshold values are not identical, but often close. Means and isolate are the same technique implemented in two different ways.

I find that for example in some thresholding techniques where the background is white that triangle thresh works best. But for doing notch filtering kapur (entropy) thresholding does best. If you have an image with truly bimodal distribution, then otsu and means are close, but otsu may be slightly better.

ASIDE: Some results and comparisons can be found in OpenCV and in Python Skimage (scipy scikit-image). See
http://docs.opencv.org/3.0-beta/doc/py_ ... lding.html
http://docs.opencv.org/3.0-beta/doc/tut ... shold.html
http://www.bogotobogo.com/python/OpenCV ... ations.php
http://scikit-image.org/docs/dev/api/sk ... lters.html
https://www.slideshare.net/sharmashrush ... lideshow=1
http://www.theobjects.com/dragonfly/dfh ... ilters.htm