vcs.maemo.org Git - opencv/blob - src/cv/cvthresh.cpp

   1 /*M///////////////////////////////////////////////////////////////////////////////////////
   2 //
   3 //  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
   4 //
   5 //  By downloading, copying, installing or using the software you agree to this license.
   6 //  If you do not agree to this license, do not download, install,
   7 //  copy or use the software.
   8 //
   9 //
  10 //                           License Agreement
  11 //                For Open Source Computer Vision Library
  12 //
  13 // Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
  14 // Copyright (C) 2009, Willow Garage Inc., all rights reserved.
  15 // Third party copyrights are property of their respective owners.
  16 //
  17 // Redistribution and use in source and binary forms, with or without modification,
  18 // are permitted provided that the following conditions are met:
  19 //
  20 //   * Redistribution's of source code must retain the above copyright notice,
  21 //     this list of conditions and the following disclaimer.
  22 //
  23 //   * Redistribution's in binary form must reproduce the above copyright notice,
  24 //     this list of conditions and the following disclaimer in the documentation
  25 //     and/or other materials provided with the distribution.
  26 //
  27 //   * The name of the copyright holders may not be used to endorse or promote products
  28 //     derived from this software without specific prior written permission.
  29 //
  30 // This software is provided by the copyright holders and contributors "as is" and
  31 // any express or implied warranties, including, but not limited to, the implied
  32 // warranties of merchantability and fitness for a particular purpose are disclaimed.
  33 // In no event shall the Intel Corporation or contributors be liable for any direct,
  34 // indirect, incidental, special, exemplary, or consequential damages
  35 // (including, but not limited to, procurement of substitute goods or services;
  36 // loss of use, data, or profits; or business interruption) however caused
  37 // and on any theory of liability, whether in contract, strict liability,
  38 // or tort (including negligence or otherwise) arising in any way out of
  39 // the use of this software, even if advised of the possibility of such damage.
  40 //
  41 //M*/
  42
  43 #include "_cv.h"
  44
  45 namespace cv
  46 {
  47
  48 static void
  49 thresh_8u( const Mat& _src, Mat& _dst, uchar thresh, uchar maxval, int type )
  50 {
  51     int i, j;
  52     uchar tab[256];
  53     Size roi = _src.size();
  54     roi.width *= _src.channels();
  55 #if CV_SSE2
  56     __m128i _x80 = _mm_set1_epi8('\x80');
  57     __m128i thresh_u = _mm_set1_epi8(thresh);
  58     __m128i thresh_s = _mm_set1_epi8(thresh ^ 0x80);
  59     __m128i maxval_ = _mm_set1_epi8(maxval);
  60 #endif
  61
  62     if( _src.isContinuous() && _dst.isContinuous() )
  63     {
  64         roi.width *= roi.height;
  65         roi.height = 1;
  66     }
  67
  68     switch( type )
  69     {
  70     case THRESH_BINARY:
  71         for( i = 0; i <= thresh; i++ )
  72             tab[i] = 0;
  73         for( ; i < 256; i++ )
  74             tab[i] = maxval;
  75         break;
  76     case THRESH_BINARY_INV:
  77         for( i = 0; i <= thresh; i++ )
  78             tab[i] = maxval;
  79         for( ; i < 256; i++ )
  80             tab[i] = 0;
  81         break;
  82     case THRESH_TRUNC:
  83         for( i = 0; i <= thresh; i++ )
  84             tab[i] = (uchar)i;
  85         for( ; i < 256; i++ )
  86             tab[i] = thresh;
  87         break;
  88     case THRESH_TOZERO:
  89         for( i = 0; i <= thresh; i++ )
  90             tab[i] = 0;
  91         for( ; i < 256; i++ )
  92             tab[i] = (uchar)i;
  93         break;
  94     case THRESH_TOZERO_INV:
  95         for( i = 0; i <= thresh; i++ )
  96             tab[i] = (uchar)i;
  97         for( ; i < 256; i++ )
  98             tab[i] = 0;
  99         break;
 100     default:
 101         CV_Error( CV_StsBadArg, "Unknown threshold type" );
 102     }
 103
 104     for( i = 0; i < roi.height; i++ )
 105     {
 106         const uchar* src = (const uchar*)(_src.data + _src.step*i);
 107         uchar* dst = (uchar*)(_dst.data + _dst.step*i);
 108         j = 0;
 109
 110     #if CV_SSE2
 111         switch( type )
 112         {
 113         case THRESH_BINARY:
 114             for( ; j <= roi.width - 32; j += 32 )
 115             {
 116                 __m128i v0, v1;
 117                 v0 = _mm_loadu_si128( (const __m128i*)(src + j) );
 118                 v1 = _mm_loadu_si128( (const __m128i*)(src + j + 16) );
 119                 v0 = _mm_cmpgt_epi8( _mm_xor_si128(v0, _x80), thresh_s );
 120                 v1 = _mm_cmpgt_epi8( _mm_xor_si128(v1, _x80), thresh_s );
 121                 v0 = _mm_and_si128( v0, maxval_ );
 122                 v1 = _mm_and_si128( v1, maxval_ );
 123                 _mm_storeu_si128( (__m128i*)(dst + j), v0 );
 124                 _mm_storeu_si128( (__m128i*)(dst + j + 16), v1 );
 125             }
 126
 127             for( ; j <= roi.width - 8; j += 8 )
 128             {
 129                 __m128i v0 = _mm_loadl_epi64( (const __m128i*)(src + j) );
 130                 v0 = _mm_cmpgt_epi8( _mm_xor_si128(v0, _x80), thresh_s );
 131                 v0 = _mm_and_si128( v0, maxval_ );
 132                 _mm_storel_epi64( (__m128i*)(dst + j), v0 );
 133             }
 134             break;
 135
 136         case THRESH_BINARY_INV:
 137             for( ; j <= roi.width - 32; j += 32 )
 138             {
 139                 __m128i v0, v1;
 140                 v0 = _mm_loadu_si128( (const __m128i*)(src + j) );
 141                 v1 = _mm_loadu_si128( (const __m128i*)(src + j + 16) );
 142                 v0 = _mm_cmpgt_epi8( _mm_xor_si128(v0, _x80), thresh_s );
 143                 v1 = _mm_cmpgt_epi8( _mm_xor_si128(v1, _x80), thresh_s );
 144                 v0 = _mm_andnot_si128( v0, maxval_ );
 145                 v1 = _mm_andnot_si128( v1, maxval_ );
 146                 _mm_storeu_si128( (__m128i*)(dst + j), v0 );
 147                 _mm_storeu_si128( (__m128i*)(dst + j + 16), v1 );
 148             }
 149
 150             for( ; j <= roi.width - 8; j += 8 )
 151             {
 152                 __m128i v0 = _mm_loadl_epi64( (const __m128i*)(src + j) );
 153                 v0 = _mm_cmpgt_epi8( _mm_xor_si128(v0, _x80), thresh_s );
 154                 v0 = _mm_andnot_si128( v0, maxval_ );
 155                 _mm_storel_epi64( (__m128i*)(dst + j), v0 );
 156             }
 157             break;
 158
 159         case THRESH_TRUNC:
 160             for( ; j <= roi.width - 32; j += 32 )
 161             {
 162                 __m128i v0, v1;
 163                 v0 = _mm_loadu_si128( (const __m128i*)(src + j) );
 164                 v1 = _mm_loadu_si128( (const __m128i*)(src + j + 16) );
 165                 v0 = _mm_subs_epu8( v0, _mm_subs_epu8( v0, thresh_u ));
 166                 v1 = _mm_subs_epu8( v1, _mm_subs_epu8( v1, thresh_u ));
 167                 _mm_storeu_si128( (__m128i*)(dst + j), v0 );
 168                 _mm_storeu_si128( (__m128i*)(dst + j + 16), v1 );
 169             }
 170
 171             for( ; j <= roi.width - 8; j += 8 )
 172             {
 173                 __m128i v0 = _mm_loadl_epi64( (const __m128i*)(src + j) );
 174                 v0 = _mm_subs_epu8( v0, _mm_subs_epu8( v0, thresh_u ));
 175                 _mm_storel_epi64( (__m128i*)(dst + j), v0 );
 176             }
 177             break;
 178
 179         case THRESH_TOZERO:
 180             for( ; j <= roi.width - 32; j += 32 )
 181             {
 182                 __m128i v0, v1;
 183                 v0 = _mm_loadu_si128( (const __m128i*)(src + j) );
 184                 v1 = _mm_loadu_si128( (const __m128i*)(src + j + 16) );
 185                 v0 = _mm_and_si128( v0, _mm_cmpgt_epi8(_mm_xor_si128(v0, _x80), thresh_s ));
 186                 v1 = _mm_and_si128( v1, _mm_cmpgt_epi8(_mm_xor_si128(v1, _x80), thresh_s ));
 187                 _mm_storeu_si128( (__m128i*)(dst + j), v0 );
 188                 _mm_storeu_si128( (__m128i*)(dst + j + 16), v1 );
 189             }
 190
 191             for( ; j <= roi.width - 8; j += 8 )
 192             {
 193                 __m128i v0 = _mm_loadl_epi64( (const __m128i*)(src + j) );
 194                 v0 = _mm_and_si128( v0, _mm_cmpgt_epi8(_mm_xor_si128(v0, _x80), thresh_s ));
 195                 _mm_storel_epi64( (__m128i*)(dst + j), v0 );
 196             }
 197             break;
 198
 199         case THRESH_TOZERO_INV:
 200             for( ; j <= roi.width - 32; j += 32 )
 201             {
 202                 __m128i v0, v1;
 203                 v0 = _mm_loadu_si128( (const __m128i*)(src + j) );
 204                 v1 = _mm_loadu_si128( (const __m128i*)(src + j + 16) );
 205                 v0 = _mm_andnot_si128( _mm_cmpgt_epi8(_mm_xor_si128(v0, _x80), thresh_s ), v0 );
 206                 v1 = _mm_andnot_si128( _mm_cmpgt_epi8(_mm_xor_si128(v1, _x80), thresh_s ), v1 );
 207                 _mm_storeu_si128( (__m128i*)(dst + j), v0 );
 208                 _mm_storeu_si128( (__m128i*)(dst + j + 16), v1 );
 209             }
 210
 211             for( ; j <= roi.width - 8; j += 8 )
 212             {
 213                 __m128i v0 = _mm_loadl_epi64( (const __m128i*)(src + j) );
 214                 v0 = _mm_andnot_si128( _mm_cmpgt_epi8(_mm_xor_si128(v0, _x80), thresh_s ), v0 );
 215                 _mm_storel_epi64( (__m128i*)(dst + j), v0 );
 216             }
 217             break;
 218         }
 219     #endif
 220
 221         for( ; j <= roi.width - 4; j += 4 )
 222         {
 223             uchar t0 = tab[src[j]];
 224             uchar t1 = tab[src[j+1]];
 225
 226             dst[j] = t0;
 227             dst[j+1] = t1;
 228
 229             t0 = tab[src[j+2]];
 230             t1 = tab[src[j+3]];
 231
 232             dst[j+2] = t0;
 233             dst[j+3] = t1;
 234         }
 235
 236         for( ; j < roi.width; j++ )
 237             dst[j] = tab[src[j]];
 238     }
 239 }
 240
 241
 242 static void
 243 thresh_32f( const Mat& _src, Mat& _dst, float thresh, float maxval, int type )
 244 {
 245     int i, j;
 246     Size roi = _src.size();
 247     roi.width *= _src.channels();
 248     const float* src = (const float*)_src.data;
 249     float* dst = (float*)_dst.data;
 250     size_t src_step = _src.step/sizeof(src[0]);
 251     size_t dst_step = _dst.step/sizeof(dst[0]);
 252 #if CV_SSE2
 253     __m128 thresh4 = _mm_set1_ps(thresh), maxval4 = _mm_set1_ps(maxval);
 254 #endif
 255
 256     if( _src.isContinuous() && _dst.isContinuous() )
 257     {
 258         roi.width *= roi.height;
 259         roi.height = 1;
 260     }
 261
 262     switch( type )
 263     {
 264     case THRESH_BINARY:
 265         for( i = 0; i < roi.height; i++, src += src_step, dst += dst_step )
 266         {
 267             j = 0;
 268         #if CV_SSE2
 269             for( ; j <= roi.width - 8; j += 8 )
 270             {
 271                 __m128 v0, v1;
 272                 v0 = _mm_loadu_ps( src + j );
 273                 v1 = _mm_loadu_ps( src + j + 4 );
 274                 v0 = _mm_cmpgt_ps( v0, thresh4 );
 275                 v1 = _mm_cmpgt_ps( v1, thresh4 );
 276                 v0 = _mm_and_ps( v0, maxval4 );
 277                 v1 = _mm_and_ps( v1, maxval4 );
 278                 _mm_storeu_ps( dst + j, v0 );
 279                 _mm_storeu_ps( dst + j + 4, v1 );
 280             }
 281         #endif
 282
 283             for( ; j < roi.width; j++ )
 284                 dst[j] = src[j] > thresh ? maxval : 0;
 285         }
 286         break;
 287
 288     case THRESH_BINARY_INV:
 289         for( i = 0; i < roi.height; i++, src += src_step, dst += dst_step )
 290         {
 291             j = 0;
 292         #if CV_SSE2
 293             for( ; j <= roi.width - 8; j += 8 )
 294             {
 295                 __m128 v0, v1;
 296                 v0 = _mm_loadu_ps( src + j );
 297                 v1 = _mm_loadu_ps( src + j + 4 );
 298                 v0 = _mm_cmple_ps( v0, thresh4 );
 299                 v1 = _mm_cmple_ps( v1, thresh4 );
 300                 v0 = _mm_and_ps( v0, maxval4 );
 301                 v1 = _mm_and_ps( v1, maxval4 );
 302                 _mm_storeu_ps( dst + j, v0 );
 303                 _mm_storeu_ps( dst + j + 4, v1 );
 304             }
 305         #endif
 306
 307             for( ; j < roi.width; j++ )
 308                 dst[j] = src[j] <= thresh ? maxval : 0;
 309         }
 310         break;
 311
 312     case THRESH_TRUNC:
 313         for( i = 0; i < roi.height; i++, src += src_step, dst += dst_step )
 314         {
 315             j = 0;
 316         #if CV_SSE2
 317             for( ; j <= roi.width - 8; j += 8 )
 318             {
 319                 __m128 v0, v1;
 320                 v0 = _mm_loadu_ps( src + j );
 321                 v1 = _mm_loadu_ps( src + j + 4 );
 322                 v0 = _mm_min_ps( v0, thresh4 );
 323                 v1 = _mm_min_ps( v1, thresh4 );
 324                 _mm_storeu_ps( dst + j, v0 );
 325                 _mm_storeu_ps( dst + j + 4, v1 );
 326             }
 327         #endif
 328
 329             for( ; j < roi.width; j++ )
 330                 dst[j] = std::min(src[j], thresh);
 331         }
 332         break;
 333
 334     case THRESH_TOZERO:
 335         for( i = 0; i < roi.height; i++, src += src_step, dst += dst_step )
 336         {
 337             j = 0;
 338         #if CV_SSE2
 339             for( ; j <= roi.width - 8; j += 8 )
 340             {
 341                 __m128 v0, v1;
 342                 v0 = _mm_loadu_ps( src + j );
 343                 v1 = _mm_loadu_ps( src + j + 4 );
 344                 v0 = _mm_and_ps(v0, _mm_cmpgt_ps(v0, thresh4));
 345                 v1 = _mm_and_ps(v1, _mm_cmpgt_ps(v1, thresh4));
 346                 _mm_storeu_ps( dst + j, v0 );
 347                 _mm_storeu_ps( dst + j + 4, v1 );
 348             }
 349         #endif
 350
 351             for( ; j < roi.width; j++ )
 352             {
 353                 float v = src[j];
 354                 dst[j] = v > thresh ? v : 0;
 355             }
 356         }
 357         break;
 358
 359     case THRESH_TOZERO_INV:
 360         for( i = 0; i < roi.height; i++, src += src_step, dst += dst_step )
 361         {
 362             j = 0;
 363         #if CV_SSE2
 364             for( ; j <= roi.width - 8; j += 8 )
 365             {
 366                 __m128 v0, v1;
 367                 v0 = _mm_loadu_ps( src + j );
 368                 v1 = _mm_loadu_ps( src + j + 4 );
 369                 v0 = _mm_and_ps(v0, _mm_cmple_ps(v0, thresh4));
 370                 v1 = _mm_and_ps(v1, _mm_cmple_ps(v1, thresh4));
 371                 _mm_storeu_ps( dst + j, v0 );
 372                 _mm_storeu_ps( dst + j + 4, v1 );
 373             }
 374         #endif
 375             for( ; j < roi.width; j++ )
 376             {
 377                 float v = src[j];
 378                 dst[j] = v <= thresh ? v : 0;
 379             }
 380         }
 381         break;
 382     default:
 383         return CV_Error( CV_StsBadArg, "" );
 384     }
 385 }
 386
 387
 388 static double
 389 getThreshVal_Otsu_8u( const Mat& _src )
 390 {
 391     Size size = _src.size();
 392     if( _src.isContinuous() )
 393     {
 394         size.width *= size.height;
 395         size.height = 1;
 396     }
 397     const int N = 256;
 398     int i, j, h[N] = {0};
 399     for( i = 0; i < size.height; i++ )
 400     {
 401         const uchar* src = _src.data + _src.step*i;
 402         for( j = 0; j <= size.width - 4; j += 4 )
 403         {
 404             int v0 = src[j], v1 = src[j+1];
 405             h[v0]++; h[v1]++;
 406             v0 = src[j+2]; v1 = src[j+3];
 407             h[v0]++; h[v1]++;
 408         }
 409         for( ; j < size.width; j++ )
 410             h[src[j]]++;
 411     }
 412
 413     double mu = 0, scale = 1./(size.width*size.height);
 414     for( i = 0; i < N; i++ )
 415         mu += i*h[i];
 416
 417     mu *= scale;
 418     double mu1 = 0, q1 = 0;
 419     double max_sigma = 0, max_val = 0;
 420
 421     for( i = 0; i < N; i++ )
 422     {
 423         double p_i, q2, mu2, sigma;
 424
 425         p_i = h[i]*scale;
 426         mu1 *= q1;
 427         q1 += p_i;
 428         q2 = 1. - q1;
 429
 430         if( std::min(q1,q2) < FLT_EPSILON || std::max(q1,q2) > 1. - FLT_EPSILON )
 431             continue;
 432
 433         mu1 = (mu1 + i*p_i)/q1;
 434         mu2 = (mu - q1*mu1)/q2;
 435         sigma = q1*q2*(mu1 - mu2)*(mu1 - mu2);
 436         if( sigma > max_sigma )
 437         {
 438             max_sigma = sigma;
 439             max_val = i;
 440         }
 441     }
 442
 443     return max_val;
 444 }
 445
 446
 447 double threshold( const Mat& _src, Mat& _dst, double thresh, double maxval, int type )
 448 {
 449     bool use_otsu = (type & THRESH_OTSU) != 0;
 450     type &= THRESH_MASK;
 451
 452     if( use_otsu )
 453     {
 454         CV_Assert( _src.type() == CV_8UC1 );
 455         thresh = getThreshVal_Otsu_8u(_src);
 456     }
 457
 458     _dst.create( _src.size(), _src.type() );
 459     if( _src.depth() == CV_8U )
 460     {
 461         int ithresh = cvFloor(thresh);
 462         thresh = ithresh;
 463         int imaxval = cvRound(maxval);
 464         if( type == THRESH_TRUNC )
 465             imaxval = ithresh;
 466         imaxval = saturate_cast<uchar>(imaxval);
 467
 468         if( ithresh < 0 || ithresh >= 255 )
 469         {
 470             if( type == THRESH_BINARY || type == THRESH_BINARY_INV ||
 471                 ((type == THRESH_TRUNC || type == THRESH_TOZERO_INV) && ithresh < 0) ||
 472                 (type == THRESH_TOZERO && ithresh >= 255) )
 473             {
 474                 int v = type == THRESH_BINARY ? (ithresh >= 255 ? 0 : imaxval) :
 475                         type == THRESH_BINARY_INV ? (ithresh >= 255 ? imaxval : 0) :
 476                         type == THRESH_TRUNC ? imaxval : 0;
 477                 _dst = Scalar::all(v);
 478             }
 479             else
 480                 _src.copyTo(_dst);
 481         }
 482         else
 483             thresh_8u( _src, _dst, (uchar)ithresh, (uchar)imaxval, type );
 484     }
 485     else if( _src.depth() == CV_32F )
 486         thresh_32f( _src, _dst, (float)thresh, (float)maxval, type );
 487     else
 488         CV_Error( CV_StsUnsupportedFormat, "" );
 489
 490     return thresh;
 491 }
 492
 493
 494 void adaptiveThreshold( const Mat& _src, Mat& _dst, double maxValue,
 495                         int method, int type, int blockSize, double delta )
 496 {
 497     CV_Assert( _src.type() == CV_8UC1 );
 498     CV_Assert( blockSize % 2 == 1 && blockSize > 1 );
 499     Size size = _src.size();
 500
 501     _dst.create( size, _src.type() );
 502
 503     if( maxValue < 0 )
 504     {
 505         _dst = Scalar(0);
 506         return;
 507     }
 508
 509     Mat _mean;
 510
 511     if( _src.data != _dst.data )
 512         _mean = _dst;
 513
 514     if( method == ADAPTIVE_THRESH_MEAN_C )
 515         boxFilter( _src, _mean, _src.type(), Size(blockSize, blockSize),
 516                    Point(-1,-1), true, BORDER_REPLICATE );
 517     else if( method == ADAPTIVE_THRESH_GAUSSIAN_C )
 518         GaussianBlur( _src, _mean, Size(blockSize, blockSize), 0, 0, BORDER_REPLICATE );
 519     else
 520         CV_Error( CV_StsBadFlag, "Unknown/unsupported adaptive threshold method" );
 521
 522     int i, j;
 523     uchar imaxval = saturate_cast<uchar>(maxValue);
 524     int idelta = type == THRESH_BINARY ? cvCeil(delta) : cvFloor(delta);
 525     uchar tab[768];
 526
 527     if( type == CV_THRESH_BINARY )
 528         for( i = 0; i < 768; i++ )
 529             tab[i] = (uchar)(i - 255 > -idelta ? imaxval : 0);
 530     else if( type == CV_THRESH_BINARY_INV )
 531         for( i = 0; i < 768; i++ )
 532             tab[i] = (uchar)(i - 255 <= -idelta ? imaxval : 0);
 533     else
 534         CV_Error( CV_StsBadFlag, "Unknown/unsupported threshold type" );
 535
 536     if( _src.isContinuous() && _mean.isContinuous() && _dst.isContinuous() )
 537     {
 538         size.width *= size.height;
 539         size.height = 1;
 540     }
 541
 542     for( i = 0; i < size.height; i++ )
 543     {
 544         const uchar* src = _src.data + _src.step*i;
 545         const uchar* mean = _mean.data + _mean.step*i;
 546         uchar* dst = _dst.data + _dst.step*i;
 547
 548         for( j = 0; j < size.width; j++ )
 549             dst[j] = tab[src[j] - mean[j] + 255];
 550     }
 551 }
 552
 553 }
 554
 555 CV_IMPL double
 556 cvThreshold( const void* srcarr, void* dstarr, double thresh, double maxval, int type )
 557 {
 558     cv::Mat src = cv::cvarrToMat(srcarr), dst = cv::cvarrToMat(dstarr), dst0 = dst;
 559
 560     CV_Assert( src.size() == dst.size() && src.channels() == dst.channels() &&
 561         (src.depth() == dst.depth() || dst.depth() == CV_8U));
 562
 563     thresh = cv::threshold( src, dst, thresh, maxval, type );
 564     if( dst0.data != dst.data )
 565         dst.convertTo( dst0, dst0.depth() );
 566     return thresh;
 567 }
 568
 569
 570 CV_IMPL void
 571 cvAdaptiveThreshold( const void *srcIm, void *dstIm, double maxValue,
 572                      int method, int type, int blockSize, double delta )
 573 {
 574     cv::Mat src = cv::cvarrToMat(srcIm), dst = cv::cvarrToMat(dstIm);
 575     CV_Assert( src.size() == dst.size() && src.type() == dst.type() );
 576     cv::adaptiveThreshold( src, dst, maxValue, method, type, blockSize, delta );
 577 }
 578
 579 /* End of file. */