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
44static CvStatus CV_STDCALL
45icvThresh_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
112static CvStatus CV_STDCALL
113icvThresh_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
196static double
197icvGetThreshVal_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
289icvAndC_8u_C1R_t icvAndC_8u_C1R_p = 0;
290icvCompareC_8u_C1R_cv_t icvCompareC_8u_C1R_cv_p = 0;
291icvThreshold_GTVal_8u_C1R_t icvThreshold_GTVal_8u_C1R_p = 0;
292icvThreshold_GTVal_32f_C1R_t icvThreshold_GTVal_32f_C1R_p = 0;
293icvThreshold_LTVal_8u_C1R_t icvThreshold_LTVal_8u_C1R_p = 0;
294icvThreshold_LTVal_32f_C1R_t icvThreshold_LTVal_32f_C1R_p = 0;
295
296CV_IMPL double
297cvThreshold( 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. */
496