OpenCV threshold with mask

Question

I'm trying to use OpenCV's cv::threshold function (more specific THRESH_OTSU), only that I'd like to do it with a mask (any shape), so that the outside (background) is ignored during calculation.

Image is single channel (as it must be), red color bellow is only to mark an example polygon on an image.

I tried using adaptiveThreshold, but there are a couple of problems that make it inappropriate in my case.

Maite

Use the threshold to create a mask, and then use copyTo with the negated mask on a white image. – Medulla

Thanks, but that won't do. Otsu binarizations uses image average and some other pramateres. If I just mask it and copy to, all black pixels will still be used to calculate those parameters. – Maite

Ah ok, I got now what you mean – Medulla

Answer 1

In general, you can simply compute the threshold using cv::threshold, and then copy the src image on dst using the inverted mask.

// Apply cv::threshold on all image
thresh = cv::threshold(src, dst, thresh, maxval, type);

// Copy original image on inverted mask
src.copyTo(dst, ~mask);

With THRESH_OTSU, however, you also need to compute the threshold value only on the masked image. The following code is a modified version of static double getThreshVal_Otsu_8u(const Mat& _src) in thresh.cpp:

double otsu_8u_with_mask(const Mat1b src, const Mat1b& mask)
{
    const int N = 256;
    int M = 0;
    int i, j, h[N] = { 0 };
    for (i = 0; i < src.rows; i++)
    {
        const uchar* psrc = src.ptr(i);
        const uchar* pmask = mask.ptr(i);
        for (j = 0; j < src.cols; j++)
        {
            if (pmask[j])
            {
                h[psrc[j]]++;
                ++M;
            }
        }
    }

    double mu = 0, scale = 1. / (M);
    for (i = 0; i < N; i++)
        mu += i*(double)h[i];

    mu *= scale;
    double mu1 = 0, q1 = 0;
    double max_sigma = 0, max_val = 0;

    for (i = 0; i < N; i++)
    {
        double p_i, q2, mu2, sigma;

        p_i = h[i] * scale;
        mu1 *= q1;
        q1 += p_i;
        q2 = 1. - q1;

        if (std::min(q1, q2) < FLT_EPSILON || std::max(q1, q2) > 1. - FLT_EPSILON)
            continue;

        mu1 = (mu1 + i*p_i) / q1;
        mu2 = (mu - q1*mu1) / q2;
        sigma = q1*q2*(mu1 - mu2)*(mu1 - mu2);
        if (sigma > max_sigma)
        {
            max_sigma = sigma;
            max_val = i;
        }
    }
    return max_val;
}

You then can wrap all in a function, here called threshold_with_mask, that wraps all different cases for you. If there is no mask, or the mask is all-white, then use cv::threshold. Otherwise, use one of the above mentioned approaches. Note that this wrapper works only for CV_8UC1 images (for simplicity sake, you can easily expand it to work with other types, if needed), and accepts all THRESH_XXX combinations as original cv::threshold.

double threshold_with_mask(Mat1b& src, Mat1b& dst, double thresh, double maxval, int type, const Mat1b& mask = Mat1b())
{
    if (mask.empty() || (mask.rows == src.rows && mask.cols == src.cols && countNonZero(mask) == src.rows * src.cols))
    {
        // If empty mask, or all-white mask, use cv::threshold
        thresh = cv::threshold(src, dst, thresh, maxval, type);
    }
    else
    {
        // Use mask
        bool use_otsu = (type & THRESH_OTSU) != 0;
        if (use_otsu)
        {
            // If OTSU, get thresh value on mask only
            thresh = otsu_8u_with_mask(src, mask);
            // Remove THRESH_OTSU from type
            type &= THRESH_MASK;
        }

        // Apply cv::threshold on all image
        thresh = cv::threshold(src, dst, thresh, maxval, type);

        // Copy original image on inverted mask
        src.copyTo(dst, ~mask);
    }
    return thresh;
}

Here is the full code for reference:

#include <opencv2/opencv.hpp>
#include <iostream>
using namespace std;
using namespace cv;

// Modified from thresh.cpp
// static double getThreshVal_Otsu_8u(const Mat& _src)

double otsu_8u_with_mask(const Mat1b src, const Mat1b& mask)
{
    const int N = 256;
    int M = 0;
    int i, j, h[N] = { 0 };
    for (i = 0; i < src.rows; i++)
    {
        const uchar* psrc = src.ptr(i);
        const uchar* pmask = mask.ptr(i);
        for (j = 0; j < src.cols; j++)
        {
            if (pmask[j])
            {
                h[psrc[j]]++;
                ++M;
            }
        }
    }

    double mu = 0, scale = 1. / (M);
    for (i = 0; i < N; i++)
        mu += i*(double)h[i];

    mu *= scale;
    double mu1 = 0, q1 = 0;
    double max_sigma = 0, max_val = 0;

    for (i = 0; i < N; i++)
    {
        double p_i, q2, mu2, sigma;

        p_i = h[i] * scale;
        mu1 *= q1;
        q1 += p_i;
        q2 = 1. - q1;

        if (std::min(q1, q2) < FLT_EPSILON || std::max(q1, q2) > 1. - FLT_EPSILON)
            continue;

        mu1 = (mu1 + i*p_i) / q1;
        mu2 = (mu - q1*mu1) / q2;
        sigma = q1*q2*(mu1 - mu2)*(mu1 - mu2);
        if (sigma > max_sigma)
        {
            max_sigma = sigma;
            max_val = i;
        }
    }

    return max_val;
}

double threshold_with_mask(Mat1b& src, Mat1b& dst, double thresh, double maxval, int type, const Mat1b& mask = Mat1b())
{
    if (mask.empty() || (mask.rows == src.rows && mask.cols == src.cols && countNonZero(mask) == src.rows * src.cols))
    {
        // If empty mask, or all-white mask, use cv::threshold
        thresh = cv::threshold(src, dst, thresh, maxval, type);
    }
    else
    {
        // Use mask
        bool use_otsu = (type & THRESH_OTSU) != 0;
        if (use_otsu)
        {
            // If OTSU, get thresh value on mask only
            thresh = otsu_8u_with_mask(src, mask);
            // Remove THRESH_OTSU from type
            type &= THRESH_MASK;
        }

        // Apply cv::threshold on all image
        thresh = cv::threshold(src, dst, thresh, maxval, type);

        // Copy original image on inverted mask
        src.copyTo(dst, ~mask);
    }
    return thresh;
}


int main()
{
    // Load an image
    Mat1b img = imread("D:\\SO\\img\\nice.jpg", IMREAD_GRAYSCALE);

    // Apply OpenCV version
    Mat1b cvth;
    double cvth_value = threshold(img, cvth, 100, 255, THRESH_OTSU);

    // Create a binary mask
    Mat1b mask(img.rows, img.cols, uchar(0));
    rectangle(mask, Rect(100, 100, 200, 200), Scalar(255), CV_FILLED);

    // Apply threshold with a mask
    Mat1b th;
    double th_value = threshold_with_mask(img, th, 100, 255, THRESH_OTSU, mask);

    // Show results
    imshow("cv::threshod", cvth);
    imshow("threshold_with_balue", th);
    waitKey();

    return 0;
}

Medulla

Thank you for your very comprehensive answer. I thought an hoped there was an easier way than this, now I feel bad for not doing the heavy lifting myself. – Maite

If you want to speed the code up a bit do a parallel for loop like: src.forEach<uchar>([&](uchar& pixel, const int po[]) -> void {uchar maskPix = mask.at<uchar>(po[0], po[1]); if (maskPix > 0) {h[pixel]++; ++M; }}); – Grania

Speed up is also possible using the OpenCv Histogram Algorithm (Cv2.CalcHist) instead of doing it by yourself. This function handles a mask as input. – Propane

@Propane this code is just a minor fix to the original otsu code. If a call to calcHist is useful for speedup, you would probably find it also in the original code. But you can try and evaluate. Please let me know ;) – Medulla

I replaced the first for-loop with cv::CalcHist(..), and it works great. As fast as the standard cv::threshold(..) algorithm. – Propane