Enable Otsu thresholding for CV_16UC1 images

modules/imgproc/src/thresh.cpp

@@ -1125,31 +1125,24 @@ static bool ipp_getThreshVal_Otsu_8u( const unsigned char* _src, int step, Size
 }
 #endif
 
-static double
+template<typename T, size_t BinsOnStack = 0u>
-getThreshVal_Otsu_8u( const Mat& _src )
+static double getThreshVal_Otsu( const Mat& _src, const Size& size)
 {
-    Size size = _src.size();
+    const int N = std::numeric_limits<T>::max() + 1;
-    int step = (int) _src.step;
+    int i, j;
-    if( _src.isContinuous() )
-    {
-        size.width *= size.height;
-        size.height = 1;
-        step = size.width;
-    }
-
-#ifdef HAVE_IPP
-    unsigned char thresh = 0;
-    CV_IPP_RUN_FAST(ipp_getThreshVal_Otsu_8u(_src.ptr(), step, size, thresh), thresh);
-#endif
-
-    const int N = 256;
-    int i, j, h[N] = {0};
     #if CV_ENABLE_UNROLLED
-    int h_unrolled[3][N] = {};
+    AutoBuffer<int, 4 * BinsOnStack> hBuf(4 * N);
+    #else
+    AutoBuffer<int, BinsOnStack> hBuf(N);
+    #endif
+    memset(hBuf.data(), 0, hBuf.size() * sizeof(int));
+    int* h = hBuf.data();
+    #if CV_ENABLE_UNROLLED
+    int* h_unrolled[3] = {h + N, h + 2 * N, h + 3 * N };
     #endif
     for( i = 0; i < size.height; i++ )
     {
-        const uchar* src = _src.ptr() + step*i;
+        const T* src = _src.ptr<T>(i, 0);
         j = 0;
         #if CV_ENABLE_UNROLLED
         for( ; j <= size.width - 4; j += 4 )
@@ -1177,7 +1170,7 @@ getThreshVal_Otsu_8u( const Mat& _src )
     double mu1 = 0, q1 = 0;
     double max_sigma = 0, max_val = 0;
 
-    for( i = 0; i < N; i++ )
+    for(i = 0; i < N; i++ )
     {
         double p_i, q2, mu2, sigma;
 
@@ -1198,10 +1191,44 @@ getThreshVal_Otsu_8u( const Mat& _src )
             max_val = i;
         }
     }
-
     return max_val;
 }
 
+static double
+getThreshVal_Otsu_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;
+    }
+
+#ifdef HAVE_IPP
+    unsigned char thresh = 0;
+    CV_IPP_RUN_FAST(ipp_getThreshVal_Otsu_8u(_src.ptr(), step, size, thresh), thresh);
+#else
+    CV_UNUSED(step);
+#endif
+
+    return getThreshVal_Otsu<uchar, 256u>(_src, size);
+}
+
+static double
+getThreshVal_Otsu_16u( const Mat& _src )
+{
+    Size size = _src.size();
+    if( _src.isContinuous() )
+    {
+        size.width *= size.height;
+        size.height = 1;
+    }
+
+    return getThreshVal_Otsu<ushort>(_src, size);
+}
+
 static double
 getThreshVal_Triangle_8u( const Mat& _src )
 {
@@ -1526,8 +1553,10 @@ double cv::threshold( InputArray _src, OutputArray _dst, double thresh, double m
     CV_Assert( automatic_thresh != (CV_THRESH_OTSU | CV_THRESH_TRIANGLE) );
     if( automatic_thresh == CV_THRESH_OTSU )
     {
-        CV_Assert( src.type() == CV_8UC1 );
+        int src_type = src.type();
-        thresh = getThreshVal_Otsu_8u( src );
+        CV_CheckType(src_type, src_type == CV_8UC1 || src_type == CV_16UC1, "THRESH_OTSU mode");
+        thresh = src.type() == CV_8UC1 ? getThreshVal_Otsu_8u( src )
+                                       : getThreshVal_Otsu_16u( src );
     }
     else if( automatic_thresh == CV_THRESH_TRIANGLE )
     {

modules/imgproc/test/test_thresh.cpp

@@ -46,7 +46,7 @@ namespace opencv_test { namespace {
 class CV_ThreshTest : public cvtest::ArrayTest
 {
 public:
-    CV_ThreshTest();
+    CV_ThreshTest(int test_type = 0);
 
 protected:
     void get_test_array_types_and_sizes( int test_case_idx, vector<vector<Size> >& sizes, vector<vector<int> >& types );
@@ -57,16 +57,22 @@ protected:
     int thresh_type;
     double thresh_val;
     double max_val;
+    int extra_type;
 };
 
 
-CV_ThreshTest::CV_ThreshTest()
+CV_ThreshTest::CV_ThreshTest(int test_type)
 {
+    CV_Assert( (test_type & CV_THRESH_MASK) == 0 );
     test_array[INPUT].push_back(NULL);
     test_array[OUTPUT].push_back(NULL);
     test_array[REF_OUTPUT].push_back(NULL);
     optional_mask = false;
     element_wise_relative_error = true;
+    extra_type = test_type;
+    // Reduce number of test with automated thresholding
+    if (extra_type != 0)
+        test_case_count = 250;
 }
 
 
@@ -78,6 +84,12 @@ void CV_ThreshTest::get_test_array_types_and_sizes( int test_case_idx,
     cvtest::ArrayTest::get_test_array_types_and_sizes( test_case_idx, sizes, types );
     depth = depth == 0 ? CV_8U : depth == 1 ? CV_16S : depth == 2 ? CV_16U : depth == 3 ? CV_32F : CV_64F;
 
+    if ( extra_type == CV_THRESH_OTSU )
+    {
+        depth = cvtest::randInt(rng) % 2 == 0 ? CV_8U : CV_16U;
+        cn = 1;
+    }
+
     types[INPUT][0] = types[OUTPUT][0] = types[REF_OUTPUT][0] = CV_MAKETYPE(depth,cn);
     thresh_type = cvtest::randInt(rng) % 5;
 
@@ -123,18 +135,73 @@ double CV_ThreshTest::get_success_error_level( int /*test_case_idx*/, int /*i*/,
 void CV_ThreshTest::run_func()
 {
     cvThreshold( test_array[INPUT][0], test_array[OUTPUT][0],
-                 thresh_val, max_val, thresh_type );
+                 thresh_val, max_val, thresh_type | extra_type);
 }
 
 
+static double compute_otsu_thresh(const Mat& _src)
+{
+    int depth = _src.depth();
+    int width = _src.cols, height = _src.rows;
+    const int N = 65536;
+    std::vector<int> h(N, 0);
+    int i, j;
+    double mu = 0, scale = 1./(width*height);
+    for(i = 0; i < height; ++i)
+    {
+        for(j = 0; j < width; ++j)
+        {
+          const int val  = depth == CV_16UC1 ? (int)_src.at<ushort>(i, j) : (int)_src.at<uchar>(i,j);
+          h[val]++;
+        }
+    }
+    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;
+}
+
 static void test_threshold( const Mat& _src, Mat& _dst,
-                            double thresh, double maxval, int thresh_type )
+                            double thresh, double maxval, int thresh_type, int extra_type )
 {
     int i, j;
     int depth = _src.depth(), cn = _src.channels();
     int width_n = _src.cols*cn, height = _src.rows;
     int ithresh = cvFloor(thresh);
     int imaxval, ithresh2;
+    if (extra_type == CV_THRESH_OTSU)
+    {
+        thresh = compute_otsu_thresh(_src);
+        ithresh = cvFloor(thresh);
+    }
 
     if( depth == CV_8U )
     {
@@ -157,7 +224,7 @@ static void test_threshold( const Mat& _src, Mat& _dst,
         imaxval = cvRound(maxval);
     }
 
-    assert( depth == CV_8U || depth == CV_16S || depth == CV_16U || depth == CV_32F || depth == CV_64F );
+    CV_Assert( depth == CV_8U || depth == CV_16S || depth == CV_16U || depth == CV_32F || depth == CV_64F );
 
     switch( thresh_type )
     {
@@ -415,10 +482,11 @@ static void test_threshold( const Mat& _src, Mat& _dst,
 void CV_ThreshTest::prepare_to_validation( int /*test_case_idx*/ )
 {
     test_threshold( test_mat[INPUT][0], test_mat[REF_OUTPUT][0],
-                   thresh_val, max_val, thresh_type );
+                   thresh_val, max_val, thresh_type, extra_type );
 }
 
 TEST(Imgproc_Threshold, accuracy) { CV_ThreshTest test; test.safe_run(); }
+TEST(Imgproc_Threshold, accuracyOtsu) { CV_ThreshTest test(CV_THRESH_OTSU); test.safe_run(); }
 
 BIGDATA_TEST(Imgproc_Threshold, huge)
 {