1 /*M///////////////////////////////////////////////////////////////////////////////////////
3 // IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
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.
10 // Intel License Agreement
11 // For Open Source Computer Vision Library
13 // Copyright (C) 2000, Intel Corporation, all rights reserved.
14 // Third party copyrights are property of their respective owners.
16 // Redistribution and use in source and binary forms, with or without modification,
17 // are permitted provided that the following conditions are met:
19 // * Redistribution's of source code must retain the above copyright notice,
20 // this list of conditions and the following disclaimer.
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.
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.
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.
41 #include "precomp.hpp"
43 #define ICV_DIST_SHIFT 16
44 #define ICV_INIT_DIST0 (INT_MAX >> 2)
47 icvInitTopBottom( int* temp, int tempstep, CvSize size, int border )
50 for( i = 0; i < border; i++ )
52 int* ttop = (int*)(temp + i*tempstep);
53 int* tbottom = (int*)(temp + (size.height + border*2 - i - 1)*tempstep);
55 for( j = 0; j < size.width + border*2; j++ )
57 ttop[j] = ICV_INIT_DIST0;
58 tbottom[j] = ICV_INIT_DIST0;
66 static CvStatus CV_STDCALL
67 icvDistanceTransform_3x3_C1R( const uchar* src, int srcstep, int* temp,
68 int step, float* dist, int dststep, CvSize size, const float* metrics )
72 const int HV_DIST = CV_FLT_TO_FIX( metrics[0], ICV_DIST_SHIFT );
73 const int DIAG_DIST = CV_FLT_TO_FIX( metrics[1], ICV_DIST_SHIFT );
74 const float scale = 1.f/(1 << ICV_DIST_SHIFT);
76 srcstep /= sizeof(src[0]);
77 step /= sizeof(temp[0]);
78 dststep /= sizeof(dist[0]);
80 icvInitTopBottom( temp, step, size, BORDER );
83 for( i = 0; i < size.height; i++ )
85 const uchar* s = src + i*srcstep;
86 int* tmp = (int*)(temp + (i+BORDER)*step) + BORDER;
88 for( j = 0; j < BORDER; j++ )
89 tmp[-j-1] = tmp[size.width + j] = ICV_INIT_DIST0;
91 for( j = 0; j < size.width; j++ )
97 int t0 = tmp[j-step-1] + DIAG_DIST;
98 int t = tmp[j-step] + HV_DIST;
100 t = tmp[j-step+1] + DIAG_DIST;
102 t = tmp[j-1] + HV_DIST;
110 for( i = size.height - 1; i >= 0; i-- )
112 float* d = (float*)(dist + i*dststep);
113 int* tmp = (int*)(temp + (i+BORDER)*step) + BORDER;
115 for( j = size.width - 1; j >= 0; j-- )
120 int t = tmp[j+step+1] + DIAG_DIST;
122 t = tmp[j+step] + HV_DIST;
124 t = tmp[j+step-1] + DIAG_DIST;
126 t = tmp[j+1] + HV_DIST;
130 d[j] = (float)(t0 * scale);
138 static CvStatus CV_STDCALL
139 icvDistanceTransform_5x5_C1R( const uchar* src, int srcstep, int* temp,
140 int step, float* dist, int dststep, CvSize size, const float* metrics )
142 const int BORDER = 2;
144 const int HV_DIST = CV_FLT_TO_FIX( metrics[0], ICV_DIST_SHIFT );
145 const int DIAG_DIST = CV_FLT_TO_FIX( metrics[1], ICV_DIST_SHIFT );
146 const int LONG_DIST = CV_FLT_TO_FIX( metrics[2], ICV_DIST_SHIFT );
147 const float scale = 1.f/(1 << ICV_DIST_SHIFT);
149 srcstep /= sizeof(src[0]);
150 step /= sizeof(temp[0]);
151 dststep /= sizeof(dist[0]);
153 icvInitTopBottom( temp, step, size, BORDER );
156 for( i = 0; i < size.height; i++ )
158 const uchar* s = src + i*srcstep;
159 int* tmp = (int*)(temp + (i+BORDER)*step) + BORDER;
161 for( j = 0; j < BORDER; j++ )
162 tmp[-j-1] = tmp[size.width + j] = ICV_INIT_DIST0;
164 for( j = 0; j < size.width; j++ )
170 int t0 = tmp[j-step*2-1] + LONG_DIST;
171 int t = tmp[j-step*2+1] + LONG_DIST;
173 t = tmp[j-step-2] + LONG_DIST;
175 t = tmp[j-step-1] + DIAG_DIST;
177 t = tmp[j-step] + HV_DIST;
179 t = tmp[j-step+1] + DIAG_DIST;
181 t = tmp[j-step+2] + LONG_DIST;
183 t = tmp[j-1] + HV_DIST;
191 for( i = size.height - 1; i >= 0; i-- )
193 float* d = (float*)(dist + i*dststep);
194 int* tmp = (int*)(temp + (i+BORDER)*step) + BORDER;
196 for( j = size.width - 1; j >= 0; j-- )
201 int t = tmp[j+step*2+1] + LONG_DIST;
203 t = tmp[j+step*2-1] + LONG_DIST;
205 t = tmp[j+step+2] + LONG_DIST;
207 t = tmp[j+step+1] + DIAG_DIST;
209 t = tmp[j+step] + HV_DIST;
211 t = tmp[j+step-1] + DIAG_DIST;
213 t = tmp[j+step-2] + LONG_DIST;
215 t = tmp[j+1] + HV_DIST;
219 d[j] = (float)(t0 * scale);
227 static CvStatus CV_STDCALL
228 icvDistanceTransformEx_5x5_C1R( const uchar* src, int srcstep, int* temp,
229 int step, float* dist, int dststep, int* labels, int lstep,
230 CvSize size, const float* metrics )
232 const int BORDER = 2;
235 const int HV_DIST = CV_FLT_TO_FIX( metrics[0], ICV_DIST_SHIFT );
236 const int DIAG_DIST = CV_FLT_TO_FIX( metrics[1], ICV_DIST_SHIFT );
237 const int LONG_DIST = CV_FLT_TO_FIX( metrics[2], ICV_DIST_SHIFT );
238 const float scale = 1.f/(1 << ICV_DIST_SHIFT);
240 srcstep /= sizeof(src[0]);
241 step /= sizeof(temp[0]);
242 dststep /= sizeof(dist[0]);
243 lstep /= sizeof(labels[0]);
245 icvInitTopBottom( temp, step, size, BORDER );
248 for( i = 0; i < size.height; i++ )
250 const uchar* s = src + i*srcstep;
251 int* tmp = (int*)(temp + (i+BORDER)*step) + BORDER;
252 int* lls = (int*)(labels + i*lstep);
254 for( j = 0; j < BORDER; j++ )
255 tmp[-j-1] = tmp[size.width + j] = ICV_INIT_DIST0;
257 for( j = 0; j < size.width; j++ )
262 //assert( lls[j] != 0 );
266 int t0 = ICV_INIT_DIST0, t;
269 t = tmp[j-step*2-1] + LONG_DIST;
273 l0 = lls[j-lstep*2-1];
275 t = tmp[j-step*2+1] + LONG_DIST;
279 l0 = lls[j-lstep*2+1];
281 t = tmp[j-step-2] + LONG_DIST;
287 t = tmp[j-step-1] + DIAG_DIST;
293 t = tmp[j-step] + HV_DIST;
299 t = tmp[j-step+1] + DIAG_DIST;
305 t = tmp[j-step+2] + LONG_DIST;
311 t = tmp[j-1] + HV_DIST;
325 for( i = size.height - 1; i >= 0; i-- )
327 float* d = (float*)(dist + i*dststep);
328 int* tmp = (int*)(temp + (i+BORDER)*step) + BORDER;
329 int* lls = (int*)(labels + i*lstep);
331 for( j = size.width - 1; j >= 0; j-- )
337 int t = tmp[j+step*2+1] + LONG_DIST;
341 l0 = lls[j+lstep*2+1];
343 t = tmp[j+step*2-1] + LONG_DIST;
347 l0 = lls[j+lstep*2-1];
349 t = tmp[j+step+2] + LONG_DIST;
355 t = tmp[j+step+1] + DIAG_DIST;
361 t = tmp[j+step] + HV_DIST;
367 t = tmp[j+step-1] + DIAG_DIST;
373 t = tmp[j+step-2] + LONG_DIST;
379 t = tmp[j+1] + HV_DIST;
388 d[j] = (float)(t0 * scale);
397 icvGetDistanceTransformMask( int maskType, float *metrics )
400 return CV_NULLPTR_ERR;
416 metrics[1] = 1.3693f;
434 metrics[2] = 2.1969f;
437 return CV_BADRANGE_ERR;
446 struct DTColumnInvoker : ParallelLoopBody
448 DTColumnInvoker( const CvMat* _src, CvMat* _dst, const int* _sat_tab, const float* _sqr_tab)
452 sat_tab = _sat_tab + src->rows*2 + 1;
456 void operator()( const Range& range ) const
458 int i, i1 = range.start, i2 = range.end;
460 size_t sstep = src->step, dstep = dst->step/sizeof(float);
461 AutoBuffer<int> _d(m);
464 for( i = i1; i < i2; i++ )
466 const uchar* sptr = src->data.ptr + i + (m-1)*sstep;
467 float* dptr = dst->data.fl + i;
470 for( j = m-1; j >= 0; j--, sptr -= sstep )
472 dist = (dist + 1) & (sptr[0] == 0 ? 0 : -1);
477 for( j = 0; j < m; j++, dptr += dstep )
479 dist = dist + 1 - sat_tab[dist - d[j]];
481 dptr[0] = sqr_tab[dist];
489 const float* sqr_tab;
493 struct DTRowInvoker : ParallelLoopBody
495 DTRowInvoker( CvMat* _dst, const float* _sqr_tab, const float* _inv_tab )
502 void operator()( const Range& range ) const
504 const float inf = 1e15f;
505 int i, i1 = range.start, i2 = range.end;
507 AutoBuffer<uchar> _buf((n+2)*2*sizeof(float) + (n+2)*sizeof(int));
508 float* f = (float*)(uchar*)_buf;
510 int* v = alignPtr((int*)(z + n + 1), sizeof(int));
512 for( i = i1; i < i2; i++ )
514 float* d = (float*)(dst->data.ptr + i*dst->step);
522 for( q = 1, k = 0; q < n; q++ )
530 float s = (fq + sqr_tab[q] - d[p] - sqr_tab[p])*inv_tab[q - p];
542 for( q = 0, k = 0; q < n; q++ )
547 d[q] = std::sqrt(sqr_tab[std::abs(q - p)] + f[p]);
553 const float* sqr_tab;
554 const float* inv_tab;
560 icvTrueDistTrans( const CvMat* src, CvMat* dst )
562 const float inf = 1e15f;
564 if( !CV_ARE_SIZES_EQ( src, dst ))
565 CV_Error( CV_StsUnmatchedSizes, "" );
567 if( CV_MAT_TYPE(src->type) != CV_8UC1 ||
568 CV_MAT_TYPE(dst->type) != CV_32FC1 )
569 CV_Error( CV_StsUnsupportedFormat,
570 "The input image must have 8uC1 type and the output one must have 32fC1 type" );
572 int i, m = src->rows, n = src->cols;
574 cv::AutoBuffer<uchar> _buf(std::max(m*2*sizeof(float) + (m*3+1)*sizeof(int), n*2*sizeof(float)));
575 // stage 1: compute 1d distance transform of each column
576 float* sqr_tab = (float*)(uchar*)_buf;
577 int* sat_tab = cv::alignPtr((int*)(sqr_tab + m*2), sizeof(int));
580 for( i = 0; i < m; i++ )
581 sqr_tab[i] = (float)(i*i);
582 for( i = m; i < m*2; i++ )
584 for( i = 0; i < shift; i++ )
586 for( ; i <= m*3; i++ )
587 sat_tab[i] = i - shift;
589 cv::parallel_for_(cv::Range(0, n), cv::DTColumnInvoker(src, dst, sat_tab, sqr_tab));
591 // stage 2: compute modified distance transform for each row
592 float* inv_tab = sqr_tab + n;
594 inv_tab[0] = sqr_tab[0] = 0.f;
595 for( i = 1; i < n; i++ )
597 inv_tab[i] = (float)(0.5/i);
598 sqr_tab[i] = (float)(i*i);
601 cv::parallel_for_(cv::Range(0, m), cv::DTRowInvoker(dst, sqr_tab, inv_tab));
605 /*********************************** IPP functions *********************************/
607 typedef CvStatus (CV_STDCALL * CvIPPDistTransFunc)( const uchar* src, int srcstep,
608 void* dst, int dststep,
609 CvSize size, const void* metrics );
611 typedef CvStatus (CV_STDCALL * CvIPPDistTransFunc2)( uchar* src, int srcstep,
612 CvSize size, const int* metrics );
614 /***********************************************************************************/
616 typedef CvStatus (CV_STDCALL * CvDistTransFunc)( const uchar* src, int srcstep,
617 int* temp, int tempstep,
618 float* dst, int dststep,
619 CvSize size, const float* metrics );
622 /****************************************************************************************\
623 Non-inplace and Inplace 8u->8u Distance Transform for CityBlock (a.k.a. L1) metric
624 (C) 2006 by Jay Stavinzky.
625 \****************************************************************************************/
628 /* 8-bit grayscale distance transform function */
630 icvDistanceATS_L1_8u( const CvMat* src, CvMat* dst )
632 int width = src->cols, height = src->rows;
638 const uchar *sbase = src->data.ptr;
639 uchar *dbase = dst->data.ptr;
640 int srcstep = src->step;
641 int dststep = dst->step;
643 CV_Assert( CV_IS_MASK_ARR( src ) && CV_MAT_TYPE( dst->type ) == CV_8UC1 );
644 CV_Assert( CV_ARE_SIZES_EQ( src, dst ));
646 ////////////////////// forward scan ////////////////////////
647 for( x = 0; x < 256; x++ )
648 lut[x] = CV_CAST_8U(x+1);
650 //init first pixel to max (we're going to be skipping it)
651 dbase[0] = (uchar)(sbase[0] == 0 ? 0 : 255);
653 //first row (scan west only, skip first pixel)
654 for( x = 1; x < width; x++ )
655 dbase[x] = (uchar)(sbase[x] == 0 ? 0 : lut[dbase[x-1]]);
657 for( y = 1; y < height; y++ )
662 //for left edge, scan north only
663 a = sbase[0] == 0 ? 0 : lut[dbase[-dststep]];
666 for( x = 1; x < width; x++ )
668 a = sbase[x] == 0 ? 0 : lut[MIN(a, dbase[x - dststep])];
673 ////////////////////// backward scan ///////////////////////
677 // do last row east pixel scan here (skip bottom right pixel)
678 for( x = width - 2; x >= 0; x-- )
681 dbase[x] = (uchar)(CV_CALC_MIN_8U(a, dbase[x]));
684 // right edge is the only error case
685 for( y = height - 2; y >= 0; y-- )
690 a = lut[dbase[width-1+dststep]];
691 dbase[width-1] = (uchar)(MIN(a, dbase[width-1]));
693 for( x = width - 2; x >= 0; x-- )
695 int b = dbase[x+dststep];
697 dbase[x] = (uchar)(MIN(a, dbase[x]));
704 /* Wrapper function for distance transform group */
706 cvDistTransform( const void* srcarr, void* dstarr,
707 int distType, int maskSize,
709 void* labelsarr, int labelType )
711 float _mask[5] = {0};
712 CvMat srcstub, *src = (CvMat*)srcarr;
713 CvMat dststub, *dst = (CvMat*)dstarr;
714 CvMat lstub, *labels = (CvMat*)labelsarr;
716 src = cvGetMat( src, &srcstub );
717 dst = cvGetMat( dst, &dststub );
719 if( !CV_IS_MASK_ARR( src ) || (CV_MAT_TYPE( dst->type ) != CV_32FC1 &&
720 (CV_MAT_TYPE(dst->type) != CV_8UC1 || distType != CV_DIST_L1 || labels)) )
721 CV_Error( CV_StsUnsupportedFormat,
722 "source image must be 8uC1 and the distance map must be 32fC1 "
723 "(or 8uC1 in case of simple L1 distance transform)" );
725 if( !CV_ARE_SIZES_EQ( src, dst ))
726 CV_Error( CV_StsUnmatchedSizes, "the source and the destination images must be of the same size" );
728 if( maskSize != CV_DIST_MASK_3 && maskSize != CV_DIST_MASK_5 && maskSize != CV_DIST_MASK_PRECISE )
729 CV_Error( CV_StsBadSize, "Mask size should be 3 or 5 or 0 (presize)" );
731 if( distType == CV_DIST_C || distType == CV_DIST_L1 )
732 maskSize = !labels ? CV_DIST_MASK_3 : CV_DIST_MASK_5;
733 else if( distType == CV_DIST_L2 && labels )
734 maskSize = CV_DIST_MASK_5;
736 if( maskSize == CV_DIST_MASK_PRECISE )
738 icvTrueDistTrans( src, dst );
744 labels = cvGetMat( labels, &lstub );
745 if( CV_MAT_TYPE( labels->type ) != CV_32SC1 )
746 CV_Error( CV_StsUnsupportedFormat, "the output array of labels must be 32sC1" );
748 if( !CV_ARE_SIZES_EQ( labels, dst ))
749 CV_Error( CV_StsUnmatchedSizes, "the array of labels has a different size" );
751 if( maskSize == CV_DIST_MASK_3 )
752 CV_Error( CV_StsNotImplemented,
753 "3x3 mask can not be used for \"labeled\" distance transform. Use 5x5 mask" );
756 if( distType == CV_DIST_C || distType == CV_DIST_L1 || distType == CV_DIST_L2 )
758 icvGetDistanceTransformMask( (distType == CV_DIST_C ? 0 :
759 distType == CV_DIST_L1 ? 1 : 2) + maskSize*10, _mask );
761 else if( distType == CV_DIST_USER )
764 CV_Error( CV_StsNullPtr, "" );
766 memcpy( _mask, mask, (maskSize/2 + 1)*sizeof(float));
769 CvSize size = cvGetMatSize(src);
771 if( CV_MAT_TYPE(dst->type) == CV_8UC1 )
773 icvDistanceATS_L1_8u( src, dst );
777 int border = maskSize == CV_DIST_MASK_3 ? 1 : 2;
778 cv::Ptr<CvMat> temp = cvCreateMat( size.height + border*2, size.width + border*2, CV_32SC1 );
782 CvDistTransFunc func = maskSize == CV_DIST_MASK_3 ?
783 icvDistanceTransform_3x3_C1R :
784 icvDistanceTransform_5x5_C1R;
786 func( src->data.ptr, src->step, temp->data.i, temp->step,
787 dst->data.fl, dst->step, size, _mask );
793 if( labelType == CV_DIST_LABEL_CCOMP )
796 cv::Ptr<CvMemStorage> st = cvCreateMemStorage();
797 cv::Ptr<CvMat> src_copy = cvCreateMat( size.height+border*2, size.width+border*2, src->type );
798 cvCopyMakeBorder(src, src_copy, cvPoint(border, border), IPL_BORDER_CONSTANT, cvScalarAll(255));
799 cvCmpS( src_copy, 0, src_copy, CV_CMP_EQ );
800 cvFindContours( src_copy, st, &contours, sizeof(CvContour),
801 CV_RETR_CCOMP, CV_CHAIN_APPROX_SIMPLE, cvPoint(-border, -border));
803 for( int label = 1; contours != 0; contours = contours->h_next, label++ )
805 CvScalar area_color = cvScalarAll(label);
806 cvDrawContours( labels, contours, area_color, area_color, -255, -1, 8 );
812 for( int i = 0; i < src->rows; i++ )
814 const uchar* srcptr = src->data.ptr + src->step*i;
815 int* labelptr = (int*)(labels->data.ptr + labels->step*i);
817 for( int j = 0; j < src->cols; j++ )
823 icvDistanceTransformEx_5x5_C1R( src->data.ptr, src->step, temp->data.i, temp->step,
824 dst->data.fl, dst->step, labels->data.i, labels->step, size, _mask );
829 void cv::distanceTransform( InputArray _src, OutputArray _dst, OutputArray _labels,
830 int distanceType, int maskSize, int labelType )
832 Mat src = _src.getMat();
833 _dst.create(src.size(), CV_32F);
834 _labels.create(src.size(), CV_32S);
835 CvMat c_src = src, c_dst = _dst.getMat(), c_labels = _labels.getMat();
836 cvDistTransform(&c_src, &c_dst, distanceType, maskSize, 0, &c_labels, labelType);
839 void cv::distanceTransform( InputArray _src, OutputArray _dst,
840 int distanceType, int maskSize )
842 Mat src = _src.getMat();
843 _dst.create(src.size(), CV_32F);
844 Mat dst = _dst.getMat();
845 CvMat c_src = src, c_dst = _dst.getMat();
846 cvDistTransform(&c_src, &c_dst, distanceType, maskSize, 0, 0, -1);