For https://stackoverflow.com/q/78171957/3962537

#include <chrono>

#include <opencv2/opencv.hpp>

// ===================================================================

namespace original {

double CalculateCriteria(const cv::Mat& im, int th)
{
    cv::Mat thresholdedIm = cv::Mat::zeros(im.size(), CV_64F);
    thresholdedIm.setTo(1, im >= th);

    double nbPixels = static_cast<double>(im.total());
    double nbPixels1 = cv::countNonZero(thresholdedIm);
    double weight1 = nbPixels1 / nbPixels;
    double weight0 = 1 - weight1;

    if (weight1 == 0 || weight0 == 0) {
        return std::numeric_limits<double>::infinity();
    }

    cv::Mat valPixels1, valPixels0;
    im.copyTo(valPixels1, thresholdedIm == 1);
    im.copyTo(valPixels0, thresholdedIm == 0);

    cv::Scalar mean1, stdDev1, mean0, stdDev0;
    cv::meanStdDev(valPixels1, mean1, stdDev1);
    cv::meanStdDev(valPixels0, mean0, stdDev0);

    double var1 = stdDev1[0] * stdDev1[0];
    double var0 = stdDev0[0] * stdDev0[0];

    return weight0 * var0 + weight1 * var1;
}

cv::Mat SelectThreshold(const cv::Mat& img)
{
    int maxVal = 255;
    std::vector<double> criterias;
    for (int th = 0; th <= maxVal; ++th) {
        criterias.push_back(original::CalculateCriteria(img, th));
    }

    auto minElement = std::min_element(criterias.begin(), criterias.end());
    int bestThreshold = static_cast<int>(std::distance(criterias.begin(), minElement));

    cv::Mat binary;
    cv::threshold(img, binary, bestThreshold, 255, cv::THRESH_BINARY);

    return binary;
}

} // namespace original

// ===================================================================

namespace dans1 {

    double CalculateCriteria(const cv::Mat& im, int th)
    {
        cv::Mat thresholdedIm = (im >= th);

        double nbPixels = static_cast<double>(im.total());
        double nbPixels1 = cv::countNonZero(thresholdedIm);
        double weight1 = nbPixels1 / nbPixels;
        double weight0 = 1 - weight1;

        if (weight1 == 0 || weight0 == 0) {
            return std::numeric_limits<double>::infinity();
        }

        cv::Mat valPixels1, valPixels0;
        im.copyTo(valPixels1, thresholdedIm);
        im.copyTo(valPixels0, ~thresholdedIm);

        cv::Scalar mean1, stdDev1, mean0, stdDev0;
        cv::meanStdDev(valPixels1, mean1, stdDev1);
        cv::meanStdDev(valPixels0, mean0, stdDev0);

        double var1 = stdDev1[0] * stdDev1[0];
        double var0 = stdDev0[0] * stdDev0[0];

        return weight0 * var0 + weight1 * var1;
    }

    cv::Mat SelectThreshold(const cv::Mat& img)
    {
        int maxVal = 255;
        std::vector<double> criterias;
        for (int th = 0; th <= maxVal; ++th) {
            criterias.push_back(dans1::CalculateCriteria(img, th));
        }

        auto minElement = std::min_element(criterias.begin(), criterias.end());
        int bestThreshold = static_cast<int>(std::distance(criterias.begin(), minElement));

        cv::Mat binary;
        cv::threshold(img, binary, bestThreshold, 255, cv::THRESH_BINARY);

        return binary;
    }

} // namespace original

// ===================================================================
namespace dans2 {

cv::Mat SelectThreshold(const cv::Mat& img)
{
    int const histSize = 256;
    float range[] = { 0, 256 };
    const float* histRange[] = { range };

    cv::Mat hist;
    cv::calcHist(&img, 1, 0, cv::Mat(), hist, 1, &histSize, histRange, true, false);
    hist.convertTo(hist, CV_64FC1);

    cv::Mat hist_sums = hist.clone();
    for (int i(0); i < histSize; ++i) {
        hist_sums.at<double>(i) *= i;
    }

    cv::Mat hist_rcumm = hist.clone();
    cv::Mat hist_sums_rcumm = hist_sums.clone();

    for (int i(histSize); i > 0; --i) {
        hist_rcumm.at<double>(i - 1) += hist_rcumm.at<double>(i);
        hist_sums_rcumm.at<double>(i - 1) += hist_sums_rcumm.at<double>(i);
    }

    cv::Mat means1 = hist_sums_rcumm / img.total();
    cv::Mat means0 = (hist_sums_rcumm.at<double>(0) - hist_sums_rcumm) / img.total();

    cv::Mat var0 = cv::Mat::zeros(histSize, 1, CV_64FC1);
    cv::Mat var1 = cv::Mat::zeros(histSize, 1, CV_64FC1);
    for (int th(0); th < histSize; ++th) {
        double mean0 = means0.at<double>(th);
        double mean1 = means1.at<double>(th);
        for (int i(0); i < histSize; ++i) {
            if (i >= th) {
                var0.at<double>(th) += hist.at<double>(i) * (0 - mean0) * (0 - mean0);
                var1.at<double>(th) += hist.at<double>(i) * (i - mean1) * (i - mean1);
            } else {
                var0.at<double>(th) += hist.at<double>(i) * (i - mean0) * (i - mean0);
                var1.at<double>(th) += hist.at<double>(i) * (0 - mean1) * (0 - mean1);
            }
        }
    }

    var0 /= img.total();
    var1 /= img.total();

    cv::Mat weights1 = hist_rcumm / img.total();
    cv::Mat weights0 = 1 - weights1;

    cv::Mat criterias = weights0.mul(var0) + weights1.mul(var1);

    cv::Point minLoc;
    cv::minMaxLoc(criterias, nullptr, nullptr, &minLoc);
    double const bestThreshold = minLoc.y;

    cv::Mat binary;
    cv::threshold(img, binary, bestThreshold, 255, cv::THRESH_BINARY);

    return binary;
}

} // namespace dans2

// ===================================================================

int main(int argc, const char* argv[])
{
    cv::Mat src_img = cv::imread("sudoku.png", cv::IMREAD_GRAYSCALE);

    cv::Mat result_img, result_img2, result_img3;

    std::vector<size_t> times, times2, times3;
    for (int i = 0; i < 16; ++i) {
        {
            auto start_t = std::chrono::high_resolution_clock::now();
            result_img = original::SelectThreshold(src_img);
            auto end_t = std::chrono::high_resolution_clock::now();
            auto duration = std::chrono::duration_cast<std::chrono::microseconds>(end_t - start_t);
            times.push_back(duration.count());
        }

        {
            auto start_t = std::chrono::high_resolution_clock::now();
            result_img2 = dans1::SelectThreshold(src_img);
            auto end_t = std::chrono::high_resolution_clock::now();
            auto duration = std::chrono::duration_cast<std::chrono::microseconds>(end_t - start_t);
            times2.push_back(duration.count());
        }

        {
            auto start_t = std::chrono::high_resolution_clock::now();
            result_img3 = dans2::SelectThreshold(src_img);
            auto end_t = std::chrono::high_resolution_clock::now();
            auto duration = std::chrono::duration_cast<std::chrono::microseconds>(end_t - start_t);
            times3.push_back(duration.count());
        }
    }

    std::sort(times.begin(), times.end());
    std::sort(times2.begin(), times2.end());
    std::sort(times3.begin(), times3.end());

    std::cout << times[0] << "   " << times2[0] << "   " << times3[0] << "\n";

    cv::Mat delta2, delta3;
    cv::absdiff(result_img, result_img2, delta2);
    cv::absdiff(result_img, result_img3, delta3);
    std::cout << cv::countNonZero(delta2) << " " << cv::countNonZero(delta3) << "\n";
}