Add Triangle thresholding algorithm

Add Triangle method for automatic threshold computation next to the existing
Otsu's method. Triangle deals better with images whose histogram does not
contain dominant peak.

See paper Zack GW, Rogers WE, Latt SA.: Automatic measurement of sister
chromatid exchange frequency. J Histochem Cytochem. 1977 Jul;25(7):741-53.

modules/cudaarithm/doc/element_operations.rst

@@ -430,7 +430,7 @@ Applies a fixed-level threshold to each array element.
 
     :param maxval: Maximum value to use with  ``THRESH_BINARY`` and  ``THRESH_BINARY_INV`` threshold types.
 
-    :param type: Threshold type. For details, see  :ocv:func:`threshold` . The ``THRESH_OTSU`` threshold type is not supported.
+    :param type: Threshold type. For details, see  :ocv:func:`threshold` . The ``THRESH_OTSU`` and ``THRESH_TRIANGLE`` threshold types are not supported.
 
     :param stream: Stream for the asynchronous version.
 

modules/imgproc/doc/miscellaneous_transformations.rst

@@ -712,11 +712,11 @@ types of thresholding supported by the function. They are determined by ``type``
 
               \texttt{dst} (x,y) =  \fork{0}{if $\texttt{src}(x,y) > \texttt{thresh}$}{\texttt{src}(x,y)}{otherwise}
 
-Also, the special value ``THRESH_OTSU`` may be combined with
-one of the above values. In this case, the function determines the optimal threshold
-value using the Otsu's algorithm and uses it instead of the specified ``thresh`` .
+Also, the special values ``THRESH_OTSU`` or ``THRESH_TRIANGLE`` may be combined with
+one of the above values. In these cases, the function determines the optimal threshold
+value using the Otsu's or Triangle algorithm and uses it instead of the specified ``thresh`` .
 The function returns the computed threshold value.
-Currently, the Otsu's method is implemented only for 8-bit images.
+Currently, the Otsu's and Triangle methods are implemented only for 8-bit images.
 
 
 .. image:: pics/threshold.png

modules/imgproc/include/opencv2/imgproc.hpp

@@ -109,7 +109,8 @@ enum { THRESH_BINARY     = 0, // value = value > threshold ? max_value : 0
        THRESH_TOZERO     = 3, // value = value > threshold ? value : 0
        THRESH_TOZERO_INV = 4, // value = value > threshold ? 0 : value
        THRESH_MASK       = 7,
-       THRESH_OTSU       = 8  // use Otsu algorithm to choose the optimal threshold value
+       THRESH_OTSU       = 8, // use Otsu algorithm to choose the optimal threshold value
+       THRESH_TRIANGLE   = 16  // use Triangle algorithm to choose the optimal threshold value
      };
 
 //! adaptive threshold algorithm

modules/imgproc/include/opencv2/imgproc/types_c.h

@@ -551,8 +551,11 @@ enum
     CV_THRESH_TOZERO      =3,  /* value = value > threshold ? value : 0           */
     CV_THRESH_TOZERO_INV  =4,  /* value = value > threshold ? 0 : value           */
     CV_THRESH_MASK        =7,
-    CV_THRESH_OTSU        =8  /* use Otsu algorithm to choose the optimal threshold value;
+    CV_THRESH_OTSU        =8, /* use Otsu algorithm to choose the optimal threshold value;
                                  combine the flag with one of the above CV_THRESH_* values */
+    CV_THRESH_TRIANGLE    =16  /* use Triangle algorithm to choose the optimal threshold value;
+                                 combine the flag with one of the above CV_THRESH_* values, but not
+                                 with CV_THRESH_OTSU */
 };
 
 /* Adaptive threshold methods */

modules/imgproc/src/thresh.cpp

@@ -986,6 +986,110 @@ getThreshVal_Otsu_8u( const Mat& _src )
     return max_val;
 }
 
+static double
+getThreshVal_Triangle_8u( const Mat& _src )
+{
+    Size size = _src.size();
+    int step = (int) _src.step;
+    if( _src.isContinuous() )
+    {
+        size.width *= size.height;
+        size.height = 1;
+        step = size.width;
+    }
+
+    const int N = 256;
+    int i, j, h[N] = {0};
+    for( i = 0; i < size.height; i++ )
+    {
+        const uchar* src = _src.ptr() + step*i;
+        j = 0;
+        #if CV_ENABLE_UNROLLED
+        for( ; j <= size.width - 4; j += 4 )
+        {
+            int v0 = src[j], v1 = src[j+1];
+            h[v0]++; h[v1]++;
+            v0 = src[j+2]; v1 = src[j+3];
+            h[v0]++; h[v1]++;
+        }
+        #endif
+        for( ; j < size.width; j++ )
+            h[src[j]]++;
+    }
+
+    int left_bound = 0, right_bound = 0, max_ind = 0, max = 0;
+    int temp;
+    bool isflipped = false;
+
+    for( i = 0; i < N; i++ )
+    {
+        if( h[i] > 0 )
+        {
+            left_bound = i;
+            break;
+        }
+    }
+    if( left_bound > 0 )
+        left_bound--;
+
+    for( i = N-1; i > 0; i-- )
+    {
+        if( h[i] > 0 )
+        {
+            right_bound = i;
+            break;
+        }
+    }
+    if( right_bound < N-1 )
+        right_bound++;
+
+    for( i = 0; i < N; i++ )
+    {
+        if( h[i] > max)
+        {
+            max = h[i];
+            max_ind = i;
+        }
+    }
+
+    if( max_ind-left_bound < right_bound-max_ind)
+    {
+        isflipped = true;
+        i = 0, j = N-1;
+        while( i < j )
+        {
+            temp = h[i]; h[i] = h[j]; h[j] = temp;
+            i++; j--;
+        }
+        left_bound = N-1-right_bound;
+        max_ind = N-1-max_ind;
+    }
+
+    double thresh = left_bound;
+    double a, b, dist = 0, tempdist;
+
+    /*
+     * We do not need to compute precise distance here. Distance is maximized, so some constants can
+     * be omitted. This speeds up a computation a bit.
+     */
+    a = max; b = left_bound-max_ind;
+    for( i = left_bound+1; i <= max_ind; i++ )
+    {
+        tempdist = a*i + b*h[i];
+        if( tempdist > dist)
+        {
+            dist = tempdist;
+            thresh = i;
+        }
+    }
+    thresh--;
+
+    if( isflipped )
+        thresh = N-1-thresh;
+
+    return thresh;
+}
+
 class ThresholdRunner : public ParallelLoopBody
 {
 public:
@@ -1086,6 +1190,7 @@ double cv::threshold( InputArray _src, OutputArray _dst, double thresh, double m
 
     Mat src = _src.getMat();
     bool use_otsu = (type & THRESH_OTSU) != 0;
+    bool use_triangle = (type & THRESH_TRIANGLE) != 0;
     type &= THRESH_MASK;
 
     if( use_otsu )
@@ -1094,6 +1199,12 @@ double cv::threshold( InputArray _src, OutputArray _dst, double thresh, double m
         thresh = getThreshVal_Otsu_8u(src);
     }
 
+    if( use_triangle )
+    {
+        CV_Assert( src.type() == CV_8UC1 );
+        thresh = getThreshVal_Triangle_8u(src);
+    }
+
     _dst.create( src.size(), src.type() );
     Mat dst = _dst.getMat();