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
icvThresh_8u_C1R(const uchar * src,int src_step,uchar * dst,int dst_step,CvSize roi,uchar thresh,uchar maxval,int type)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
icvThresh_32f_C1R(const float * src,int src_step,float * dst,int dst_step,CvSize roi,float thresh,float maxval,int type)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
icvGetThreshVal_Otsu(const CvHistogram * hist)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
cvThreshold(const void * srcarr,void * dstarr,double thresh,double maxval,int type)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. */
496