git.maemo.org Git - opencv/blob - cv/src/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 //                        Intel License Agreement
  11 //                For Open Source Computer Vision Library
  12 //
  13 // Copyright (C) 2000, Intel Corporation, all rights reserved.
  14 // Third party copyrights are property of their respective owners.
  15 //
  16 // Redistribution and use in source and binary forms, with or without modification,
  17 // are permitted provided that the following conditions are met:
  18 //
  19 //   * Redistribution's of source code must retain the above copyright notice,
  20 //     this list of conditions and the following disclaimer.
  21 //
  22 //   * Redistribution's in binary form must reproduce the above copyright notice,
  23 //     this list of conditions and the following disclaimer in the documentation
  24 //     and/or other materials provided with the distribution.
  25 //
  26 //   * The name of Intel Corporation may not be used to endorse or promote products
  27 //     derived from this software without specific prior written permission.
  28 //
  29 // This software is provided by the copyright holders and contributors "as is" and
  30 // any express or implied warranties, including, but not limited to, the implied
  31 // warranties of merchantability and fitness for a particular purpose are disclaimed.
  32 // In no event shall the Intel Corporation or contributors be liable for any direct,
  33 // indirect, incidental, special, exemplary, or consequential damages
  34 // (including, but not limited to, procurement of substitute goods or services;
  35 // loss of use, data, or profits; or business interruption) however caused
  36 // and on any theory of liability, whether in contract, strict liability,
  37 // or tort (including negligence or otherwise) arising in any way out of
  38 // the use of this software, even if advised of the possibility of such damage.
  39 //
  40 //M*/
  41
  42 #include "_cv.h"
  43
  44 static CvStatus CV_STDCALL
  45 icvThresh_8u_C1R( const uchar* src, int src_step, uchar* dst, int dst_step,
  46                   CvSize roi, uchar thresh, uchar maxval, int type )
  47 {
  48     int i, j;
  49     uchar tab[256];
  50
  51     switch( type )
  52     {
  53     case CV_THRESH_BINARY:
  54         for( i = 0; i <= thresh; i++ )
  55             tab[i] = 0;
  56         for( ; i < 256; i++ )
  57             tab[i] = maxval;
  58         break;
  59     case CV_THRESH_BINARY_INV:
  60         for( i = 0; i <= thresh; i++ )
  61             tab[i] = maxval;
  62         for( ; i < 256; i++ )
  63             tab[i] = 0;
  64         break;
  65     case CV_THRESH_TRUNC:
  66         for( i = 0; i <= thresh; i++ )
  67             tab[i] = (uchar)i;
  68         for( ; i < 256; i++ )
  69             tab[i] = thresh;
  70         break;
  71     case CV_THRESH_TOZERO:
  72         for( i = 0; i <= thresh; i++ )
  73             tab[i] = 0;
  74         for( ; i < 256; i++ )
  75             tab[i] = (uchar)i;
  76         break;
  77     case CV_THRESH_TOZERO_INV:
  78         for( i = 0; i <= thresh; i++ )
  79             tab[i] = (uchar)i;
  80         for( ; i < 256; i++ )
  81             tab[i] = 0;
  82         break;
  83     default:
  84         return CV_BADFLAG_ERR;
  85     }
  86
  87     for( i = 0; i < roi.height; i++, src += src_step, dst += dst_step )
  88     {
  89         for( j = 0; j <= roi.width - 4; j += 4 )
  90         {
  91             uchar t0 = tab[src[j]];
  92             uchar t1 = tab[src[j+1]];
  93
  94             dst[j] = t0;
  95             dst[j+1] = t1;
  96
  97             t0 = tab[src[j+2]];
  98             t1 = tab[src[j+3]];
  99
 100             dst[j+2] = t0;
 101             dst[j+3] = t1;
 102         }
 103
 104         for( ; j < roi.width; j++ )
 105             dst[j] = tab[src[j]];
 106     }
 107
 108     return CV_NO_ERR;
 109 }
 110
 111
 112 static CvStatus CV_STDCALL
 113 icvThresh_32f_C1R( const float *src, int src_step, float *dst, int dst_step,
 114                    CvSize roi, float thresh, float maxval, int type )
 115 {
 116     int i, j;
 117     const int* isrc = (const int*)src;
 118     int* idst = (int*)dst;
 119     Cv32suf v;
 120     int iThresh, iMax;
 121
 122     v.f = thresh; iThresh = CV_TOGGLE_FLT(v.i);
 123     v.f = maxval; iMax = v.i;
 124
 125     src_step /= sizeof(src[0]);
 126     dst_step /= sizeof(dst[0]);
 127
 128     switch( type )
 129     {
 130     case CV_THRESH_BINARY:
 131         for( i = 0; i < roi.height; i++, isrc += src_step, idst += dst_step )
 132         {
 133             for( j = 0; j < roi.width; j++ )
 134             {
 135                 int temp = isrc[j];
 136                 idst[j] = ((CV_TOGGLE_FLT(temp) <= iThresh) - 1) & iMax;
 137             }
 138         }
 139         break;
 140
 141     case CV_THRESH_BINARY_INV:
 142         for( i = 0; i < roi.height; i++, isrc += src_step, idst += dst_step )
 143         {
 144             for( j = 0; j < roi.width; j++ )
 145             {
 146                 int temp = isrc[j];
 147                 idst[j] = ((CV_TOGGLE_FLT(temp) > iThresh) - 1) & iMax;
 148             }
 149         }
 150         break;
 151
 152     case CV_THRESH_TRUNC:
 153         for( i = 0; i < roi.height; i++, src += src_step, dst += dst_step )
 154         {
 155             for( j = 0; j < roi.width; j++ )
 156             {
 157                 float temp = src[j];
 158
 159                 if( temp > thresh )
 160                     temp = thresh;
 161                 dst[j] = temp;
 162             }
 163         }
 164         break;
 165
 166     case CV_THRESH_TOZERO:
 167         for( i = 0; i < roi.height; i++, isrc += src_step, idst += dst_step )
 168         {
 169             for( j = 0; j < roi.width; j++ )
 170             {
 171                 int temp = isrc[j];
 172                 idst[j] = ((CV_TOGGLE_FLT( temp ) <= iThresh) - 1) & temp;
 173             }
 174         }
 175         break;
 176
 177     case CV_THRESH_TOZERO_INV:
 178         for( i = 0; i < roi.height; i++, isrc += src_step, idst += dst_step )
 179         {
 180             for( j = 0; j < roi.width; j++ )
 181             {
 182                 int temp = isrc[j];
 183                 idst[j] = ((CV_TOGGLE_FLT( temp ) > iThresh) - 1) & temp;
 184             }
 185         }
 186         break;
 187
 188     default:
 189         return CV_BADFLAG_ERR;
 190     }
 191
 192     return CV_OK;
 193 }
 194
 195
 196 static double
 197 icvGetThreshVal_Otsu( const CvHistogram* hist )
 198 {
 199     double max_val = 0;
 200
 201     CV_FUNCNAME( "icvGetThreshVal_Otsu" );
 202
 203     __BEGIN__;
 204
 205     int i, count;
 206     const float* h;
 207     double sum = 0, mu = 0;
 208     bool uniform = false;
 209     double low = 0, high = 0, delta = 0;
 210     float* nu_thresh = 0;
 211     double mu1 = 0, q1 = 0;
 212     double max_sigma = 0;
 213
 214     if( !CV_IS_HIST(hist) || CV_IS_SPARSE_HIST(hist) || hist->mat.dims != 1 )
 215         CV_ERROR( CV_StsBadArg,
 216         "The histogram in Otsu method must be a valid dense 1D histogram" );
 217
 218     count = hist->mat.dim[0].size;
 219     h = (float*)cvPtr1D( hist->bins, 0 );
 220
 221     if( !CV_HIST_HAS_RANGES(hist) || CV_IS_UNIFORM_HIST(hist) )
 222     {
 223         if( CV_HIST_HAS_RANGES(hist) )
 224         {
 225             low = hist->thresh[0][0];
 226             high = hist->thresh[0][1];
 227         }
 228         else
 229         {
 230             low = 0;
 231             high = count;
 232         }
 233
 234         delta = (high-low)/count;
 235         low += delta*0.5;
 236         uniform = true;
 237     }
 238     else
 239         nu_thresh = hist->thresh2[0];
 240
 241     for( i = 0; i < count; i++ )
 242     {
 243         sum += h[i];
 244         if( uniform )
 245             mu += (i*delta + low)*h[i];
 246         else
 247             mu += (nu_thresh[i*2] + nu_thresh[i*2+1])*0.5*h[i];
 248     }
 249
 250     sum = fabs(sum) > FLT_EPSILON ? 1./sum : 0;
 251     mu *= sum;
 252
 253     mu1 = 0;
 254     q1 = 0;
 255
 256     for( i = 0; i < count; i++ )
 257     {
 258         double p_i, q2, mu2, val_i, sigma;
 259
 260         p_i = h[i]*sum;
 261         mu1 *= q1;
 262         q1 += p_i;
 263         q2 = 1. - q1;
 264
 265         if( MIN(q1,q2) < FLT_EPSILON || MAX(q1,q2) > 1. - FLT_EPSILON )
 266             continue;
 267
 268         if( uniform )
 269             val_i = i*delta + low;
 270         else
 271             val_i = (nu_thresh[i*2] + nu_thresh[i*2+1])*0.5;
 272
 273         mu1 = (mu1 + val_i*p_i)/q1;
 274         mu2 = (mu - q1*mu1)/q2;
 275         sigma = q1*q2*(mu1 - mu2)*(mu1 - mu2);
 276         if( sigma > max_sigma )
 277         {
 278             max_sigma = sigma;
 279             max_val = val_i;
 280         }
 281     }
 282
 283     __END__;
 284
 285     return max_val;
 286 }
 287
 288
 289 icvAndC_8u_C1R_t icvAndC_8u_C1R_p = 0;
 290 icvCompareC_8u_C1R_cv_t icvCompareC_8u_C1R_cv_p = 0;
 291 icvThreshold_GTVal_8u_C1R_t icvThreshold_GTVal_8u_C1R_p = 0;
 292 icvThreshold_GTVal_32f_C1R_t icvThreshold_GTVal_32f_C1R_p = 0;
 293 icvThreshold_LTVal_8u_C1R_t icvThreshold_LTVal_8u_C1R_p = 0;
 294 icvThreshold_LTVal_32f_C1R_t icvThreshold_LTVal_32f_C1R_p = 0;
 295
 296 CV_IMPL double
 297 cvThreshold( const void* srcarr, void* dstarr, double thresh, double maxval, int type )
 298 {
 299     CvHistogram* hist = 0;
 300
 301     CV_FUNCNAME( "cvThreshold" );
 302
 303     __BEGIN__;
 304
 305     CvSize roi;
 306     int src_step, dst_step;
 307     CvMat src_stub, *src = (CvMat*)srcarr;
 308     CvMat dst_stub, *dst = (CvMat*)dstarr;
 309     CvMat src0, dst0;
 310     int coi1 = 0, coi2 = 0;
 311     int ithresh, imaxval, cn;
 312     bool use_otsu;
 313
 314     CV_CALL( src = cvGetMat( src, &src_stub, &coi1 ));
 315     CV_CALL( dst = cvGetMat( dst, &dst_stub, &coi2 ));
 316
 317     if( coi1 + coi2 )
 318         CV_ERROR( CV_BadCOI, "COI is not supported by the function" );
 319
 320     if( !CV_ARE_CNS_EQ( src, dst ) )
 321         CV_ERROR( CV_StsUnmatchedFormats, "Both arrays must have equal number of channels" );
 322
 323     cn = CV_MAT_CN(src->type);
 324     if( cn > 1 )
 325     {
 326         src = cvReshape( src, &src0, 1 );
 327         dst = cvReshape( dst, &dst0, 1 );
 328     }
 329
 330     use_otsu = (type & ~CV_THRESH_MASK) == CV_THRESH_OTSU;
 331     type &= CV_THRESH_MASK;
 332
 333     if( use_otsu )
 334     {
 335         float _ranges[] = { 0, 256 };
 336         float* ranges = _ranges;
 337         int hist_size = 256;
 338         void* srcarr0 = src;
 339
 340         if( CV_MAT_TYPE(src->type) != CV_8UC1 )
 341             CV_ERROR( CV_StsNotImplemented, "Otsu method can only be used with 8uC1 images" );
 342
 343         CV_CALL( hist = cvCreateHist( 1, &hist_size, CV_HIST_ARRAY, &ranges ));
 344         cvCalcArrHist( &srcarr0, hist );
 345         thresh = cvFloor(icvGetThreshVal_Otsu( hist ));
 346     }
 347
 348     if( !CV_ARE_DEPTHS_EQ( src, dst ) )
 349     {
 350         if( CV_MAT_TYPE(dst->type) != CV_8UC1 )
 351             CV_ERROR( CV_StsUnsupportedFormat, "In case of different types destination should be 8uC1" );
 352
 353         if( type != CV_THRESH_BINARY && type != CV_THRESH_BINARY_INV )
 354             CV_ERROR( CV_StsBadArg,
 355             "In case of different types only CV_THRESH_BINARY "
 356             "and CV_THRESH_BINARY_INV thresholding types are supported" );
 357
 358         if( maxval < 0 )
 359         {
 360             CV_CALL( cvSetZero( dst ));
 361         }
 362         else
 363         {
 364             CV_CALL( cvCmpS( src, thresh, dst, type == CV_THRESH_BINARY ? CV_CMP_GT : CV_CMP_LE ));
 365             if( maxval < 255 )
 366                 CV_CALL( cvAndS( dst, cvScalarAll( maxval ), dst ));
 367         }
 368         EXIT;
 369     }
 370
 371     if( !CV_ARE_SIZES_EQ( src, dst ) )
 372         CV_ERROR( CV_StsUnmatchedSizes, "" );
 373
 374     roi = cvGetMatSize( src );
 375     if( CV_IS_MAT_CONT( src->type & dst->type ))
 376     {
 377         roi.width *= roi.height;
 378         roi.height = 1;
 379         src_step = dst_step = CV_STUB_STEP;
 380     }
 381     else
 382     {
 383         src_step = src->step;
 384         dst_step = dst->step;
 385     }
 386
 387     switch( CV_MAT_DEPTH(src->type) )
 388     {
 389     case CV_8U:
 390
 391         ithresh = cvFloor(thresh);
 392         imaxval = cvRound(maxval);
 393         if( type == CV_THRESH_TRUNC )
 394             imaxval = ithresh;
 395         imaxval = CV_CAST_8U(imaxval);
 396
 397         if( ithresh < 0 || ithresh >= 255 )
 398         {
 399             if( type == CV_THRESH_BINARY || type == CV_THRESH_BINARY_INV ||
 400                 (type == CV_THRESH_TRUNC || type == CV_THRESH_TOZERO_INV) && ithresh < 0 ||
 401                 type == CV_THRESH_TOZERO && ithresh >= 255 )
 402             {
 403                 int v = type == CV_THRESH_BINARY ? (ithresh >= 255 ? 0 : imaxval) :
 404                         type == CV_THRESH_BINARY_INV ? (ithresh >= 255 ? imaxval : 0) :
 405                         type == CV_THRESH_TRUNC ? imaxval : 0;
 406
 407                 cvSet( dst, cvScalarAll(v) );
 408                 EXIT;
 409             }
 410             else
 411             {
 412                 cvCopy( src, dst );
 413                 EXIT;
 414             }
 415         }
 416
 417         if( type == CV_THRESH_BINARY || type == CV_THRESH_BINARY_INV )
 418         {
 419             if( icvCompareC_8u_C1R_cv_p && icvAndC_8u_C1R_p )
 420             {
 421                 IPPI_CALL( icvCompareC_8u_C1R_cv_p( src->data.ptr, src_step,
 422                     (uchar)ithresh, dst->data.ptr, dst_step, roi,
 423                     type == CV_THRESH_BINARY ? cvCmpGreater : cvCmpLessEq ));
 424
 425                 if( imaxval < 255 )
 426                     IPPI_CALL( icvAndC_8u_C1R_p( dst->data.ptr, dst_step,
 427                     (uchar)imaxval, dst->data.ptr, dst_step, roi ));
 428                 EXIT;
 429             }
 430         }
 431         else if( type == CV_THRESH_TRUNC || type == CV_THRESH_TOZERO_INV )
 432         {
 433             if( icvThreshold_GTVal_8u_C1R_p )
 434             {
 435                 IPPI_CALL( icvThreshold_GTVal_8u_C1R_p( src->data.ptr, src_step,
 436                     dst->data.ptr, dst_step, roi, (uchar)ithresh,
 437                     (uchar)(type == CV_THRESH_TRUNC ? ithresh : 0) ));
 438                 EXIT;
 439             }
 440         }
 441         else
 442         {
 443             assert( type == CV_THRESH_TOZERO );
 444             if( icvThreshold_LTVal_8u_C1R_p )
 445             {
 446                 ithresh = cvFloor(thresh+1.);
 447                 ithresh = CV_CAST_8U(ithresh);
 448                 IPPI_CALL( icvThreshold_LTVal_8u_C1R_p( src->data.ptr, src_step,
 449                     dst->data.ptr, dst_step, roi, (uchar)ithresh, 0 ));
 450                 EXIT;
 451             }
 452         }
 453
 454         icvThresh_8u_C1R( src->data.ptr, src_step,
 455                           dst->data.ptr, dst_step, roi,
 456                           (uchar)ithresh, (uchar)imaxval, type );
 457         break;
 458     case CV_32F:
 459
 460         if( type == CV_THRESH_TRUNC || type == CV_THRESH_TOZERO_INV )
 461         {
 462             if( icvThreshold_GTVal_32f_C1R_p )
 463             {
 464                 IPPI_CALL( icvThreshold_GTVal_32f_C1R_p( src->data.fl, src_step,
 465                     dst->data.fl, dst_step, roi, (float)thresh,
 466                     type == CV_THRESH_TRUNC ? (float)thresh : 0 ));
 467                 EXIT;
 468             }
 469         }
 470         else if( type == CV_THRESH_TOZERO )
 471         {
 472             if( icvThreshold_LTVal_32f_C1R_p )
 473             {
 474                 IPPI_CALL( icvThreshold_LTVal_32f_C1R_p( src->data.fl, src_step,
 475                     dst->data.fl, dst_step, roi, (float)(thresh*(1 + FLT_EPSILON)), 0 ));
 476                 EXIT;
 477             }
 478         }
 479
 480         icvThresh_32f_C1R( src->data.fl, src_step, dst->data.fl, dst_step, roi,
 481                            (float)thresh, (float)maxval, type );
 482         break;
 483     default:
 484         CV_ERROR( CV_BadDepth, cvUnsupportedFormat );
 485     }
 486
 487     __END__;
 488
 489     if( hist )
 490         cvReleaseHist( &hist );
 491
 492     return thresh;
 493 }
 494
 495 /* End of file. */