operations.hpp

00001 /*M///////////////////////////////////////////////////////////////////////////////////////
00002 //
00003 //  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
00004 //
00005 //  By downloading, copying, installing or using the software you agree to this license.
00006 //  If you do not agree to this license, do not download, install,
00007 //  copy or use the software.
00008 //
00009 //
00010 //                           License Agreement
00011 //                For Open Source Computer Vision Library
00012 //
00013 // Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
00014 // Copyright (C) 2009, Willow Garage Inc., all rights reserved.
00015 // Copyright (C) 2013, OpenCV Foundation, all rights reserved.
00016 // Copyright (C) 2015, Itseez Inc., all rights reserved.
00017 // Third party copyrights are property of their respective owners.
00018 //
00019 // Redistribution and use in source and binary forms, with or without modification,
00020 // are permitted provided that the following conditions are met:
00021 //
00022 //   * Redistribution's of source code must retain the above copyright notice,
00023 //     this list of conditions and the following disclaimer.
00024 //
00025 //   * Redistribution's in binary form must reproduce the above copyright notice,
00026 //     this list of conditions and the following disclaimer in the documentation
00027 //     and/or other materials provided with the distribution.
00028 //
00029 //   * The name of the copyright holders may not be used to endorse or promote products
00030 //     derived from this software without specific prior written permission.
00031 //
00032 // This software is provided by the copyright holders and contributors "as is" and
00033 // any express or implied warranties, including, but not limited to, the implied
00034 // warranties of merchantability and fitness for a particular purpose are disclaimed.
00035 // In no event shall the Intel Corporation or contributors be liable for any direct,
00036 // indirect, incidental, special, exemplary, or consequential damages
00037 // (including, but not limited to, procurement of substitute goods or services;
00038 // loss of use, data, or profits; or business interruption) however caused
00039 // and on any theory of liability, whether in contract, strict liability,
00040 // or tort (including negligence or otherwise) arising in any way out of
00041 // the use of this software, even if advised of the possibility of such damage.
00042 //
00043 //M*/
00044 
00045 #ifndef __OPENCV_CORE_OPERATIONS_HPP__
00046 #define __OPENCV_CORE_OPERATIONS_HPP__
00047 
00048 #ifndef __cplusplus
00049 #  error operations.hpp header must be compiled as C++
00050 #endif
00051 
00052 #include <cstdio>
00053 
00054 //! @cond IGNORED
00055 
00056 namespace cv
00057 {
00058 
00059 ////////////////////////////// Matx methods depending on core API /////////////////////////////
00060 
00061 namespace internal
00062 {
00063 
00064 template<typename _Tp, int m> struct Matx_FastInvOp
00065 {
00066     bool operator()(const Matx<_Tp, m, m>& a, Matx<_Tp, m, m>& b, int method) const
00067     {
00068         Matx<_Tp, m, m> temp = a;
00069 
00070         // assume that b is all 0's on input => make it a unity matrix
00071         for( int i = 0; i < m; i++ )
00072             b(i, i) = (_Tp)1;
00073 
00074         if( method == DECOMP_CHOLESKY )
00075             return Cholesky(temp.val, m*sizeof(_Tp), m, b.val, m*sizeof(_Tp), m);
00076 
00077         return LU(temp.val, m*sizeof(_Tp), m, b.val, m*sizeof(_Tp), m) != 0;
00078     }
00079 };
00080 
00081 template<typename _Tp> struct Matx_FastInvOp<_Tp, 2>
00082 {
00083     bool operator()(const Matx<_Tp, 2, 2>& a, Matx<_Tp, 2, 2>& b, int) const
00084     {
00085         _Tp d = determinant(a);
00086         if( d == 0 )
00087             return false;
00088         d = 1/d;
00089         b(1,1) = a(0,0)*d;
00090         b(0,0) = a(1,1)*d;
00091         b(0,1) = -a(0,1)*d;
00092         b(1,0) = -a(1,0)*d;
00093         return true;
00094     }
00095 };
00096 
00097 template<typename _Tp> struct Matx_FastInvOp<_Tp, 3>
00098 {
00099     bool operator()(const Matx<_Tp, 3, 3>& a, Matx<_Tp, 3, 3>& b, int) const
00100     {
00101         _Tp d = (_Tp)determinant(a);
00102         if( d == 0 )
00103             return false;
00104         d = 1/d;
00105         b(0,0) = (a(1,1) * a(2,2) - a(1,2) * a(2,1)) * d;
00106         b(0,1) = (a(0,2) * a(2,1) - a(0,1) * a(2,2)) * d;
00107         b(0,2) = (a(0,1) * a(1,2) - a(0,2) * a(1,1)) * d;
00108 
00109         b(1,0) = (a(1,2) * a(2,0) - a(1,0) * a(2,2)) * d;
00110         b(1,1) = (a(0,0) * a(2,2) - a(0,2) * a(2,0)) * d;
00111         b(1,2) = (a(0,2) * a(1,0) - a(0,0) * a(1,2)) * d;
00112 
00113         b(2,0) = (a(1,0) * a(2,1) - a(1,1) * a(2,0)) * d;
00114         b(2,1) = (a(0,1) * a(2,0) - a(0,0) * a(2,1)) * d;
00115         b(2,2) = (a(0,0) * a(1,1) - a(0,1) * a(1,0)) * d;
00116         return true;
00117     }
00118 };
00119 
00120 
00121 template<typename _Tp, int m, int n> struct Matx_FastSolveOp
00122 {
00123     bool operator()(const Matx<_Tp, m, m>& a, const Matx<_Tp, m, n>& b,
00124                     Matx<_Tp, m, n>& x, int method) const
00125     {
00126         Matx<_Tp, m, m> temp = a;
00127         x = b;
00128         if( method == DECOMP_CHOLESKY )
00129             return Cholesky(temp.val, m*sizeof(_Tp), m, x.val, n*sizeof(_Tp), n);
00130 
00131         return LU(temp.val, m*sizeof(_Tp), m, x.val, n*sizeof(_Tp), n) != 0;
00132     }
00133 };
00134 
00135 template<typename _Tp> struct Matx_FastSolveOp<_Tp, 2, 1>
00136 {
00137     bool operator()(const Matx<_Tp, 2, 2>& a, const Matx<_Tp, 2, 1>& b,
00138                     Matx<_Tp, 2, 1>& x, int) const
00139     {
00140         _Tp d = determinant(a);
00141         if( d == 0 )
00142             return false;
00143         d = 1/d;
00144         x(0) = (b(0)*a(1,1) - b(1)*a(0,1))*d;
00145         x(1) = (b(1)*a(0,0) - b(0)*a(1,0))*d;
00146         return true;
00147     }
00148 };
00149 
00150 template<typename _Tp> struct Matx_FastSolveOp<_Tp, 3, 1>
00151 {
00152     bool operator()(const Matx<_Tp, 3, 3>& a, const Matx<_Tp, 3, 1>& b,
00153                     Matx<_Tp, 3, 1>& x, int) const
00154     {
00155         _Tp d = (_Tp)determinant(a);
00156         if( d == 0 )
00157             return false;
00158         d = 1/d;
00159         x(0) = d*(b(0)*(a(1,1)*a(2,2) - a(1,2)*a(2,1)) -
00160                 a(0,1)*(b(1)*a(2,2) - a(1,2)*b(2)) +
00161                 a(0,2)*(b(1)*a(2,1) - a(1,1)*b(2)));
00162 
00163         x(1) = d*(a(0,0)*(b(1)*a(2,2) - a(1,2)*b(2)) -
00164                 b(0)*(a(1,0)*a(2,2) - a(1,2)*a(2,0)) +
00165                 a(0,2)*(a(1,0)*b(2) - b(1)*a(2,0)));
00166 
00167         x(2) = d*(a(0,0)*(a(1,1)*b(2) - b(1)*a(2,1)) -
00168                 a(0,1)*(a(1,0)*b(2) - b(1)*a(2,0)) +
00169                 b(0)*(a(1,0)*a(2,1) - a(1,1)*a(2,0)));
00170         return true;
00171     }
00172 };
00173 
00174 } // internal
00175 
00176 template<typename _Tp, int m, int n> inline
00177 Matx<_Tp,m,n> Matx<_Tp,m,n>::randu(_Tp a, _Tp b)
00178 {
00179     Matx<_Tp,m,n> M;
00180     cv::randu(M, Scalar(a), Scalar(b));
00181     return M;
00182 }
00183 
00184 template<typename _Tp, int m, int n> inline
00185 Matx<_Tp,m,n> Matx<_Tp,m,n>::randn(_Tp a, _Tp b)
00186 {
00187     Matx<_Tp,m,n> M;
00188     cv::randn(M, Scalar(a), Scalar(b));
00189     return M;
00190 }
00191 
00192 template<typename _Tp, int m, int n> inline
00193 Matx<_Tp, n, m> Matx<_Tp, m, n>::inv(int method, bool *p_is_ok /*= NULL*/) const
00194 {
00195     Matx<_Tp, n, m> b;
00196     bool ok;
00197     if( method == DECOMP_LU || method == DECOMP_CHOLESKY )
00198         ok = cv::internal::Matx_FastInvOp<_Tp, m>()(*this, b, method);
00199     else
00200     {
00201         Mat A(*this, false), B(b, false);
00202         ok = (invert(A, B, method) != 0);
00203     }
00204     if( NULL != p_is_ok ) { *p_is_ok = ok; }
00205     return ok ? b : Matx<_Tp, n, m>::zeros();
00206 }
00207 
00208 template<typename _Tp, int m, int n> template<int l> inline
00209 Matx<_Tp, n, l> Matx<_Tp, m, n>::solve(const Matx<_Tp, m, l>& rhs, int method) const
00210 {
00211     Matx<_Tp, n, l> x;
00212     bool ok;
00213     if( method == DECOMP_LU || method == DECOMP_CHOLESKY )
00214         ok = cv::internal::Matx_FastSolveOp<_Tp, m, l>()(*this, rhs, x, method);
00215     else
00216     {
00217         Mat A(*this, false), B(rhs, false), X(x, false);
00218         ok = cv::solve(A, B, X, method);
00219     }
00220 
00221     return ok ? x : Matx<_Tp, n, l>::zeros();
00222 }
00223 
00224 
00225 
00226 ////////////////////////// Augmenting algebraic & logical operations //////////////////////////
00227 
00228 #define CV_MAT_AUG_OPERATOR1(op, cvop, A, B) \
00229     static inline A& operator op (A& a, const B& b) { cvop; return a; }
00230 
00231 #define CV_MAT_AUG_OPERATOR(op, cvop, A, B)   \
00232     CV_MAT_AUG_OPERATOR1(op, cvop, A, B)      \
00233     CV_MAT_AUG_OPERATOR1(op, cvop, const A, B)
00234 
00235 #define CV_MAT_AUG_OPERATOR_T(op, cvop, A, B)                   \
00236     template<typename _Tp> CV_MAT_AUG_OPERATOR1(op, cvop, A, B) \
00237     template<typename _Tp> CV_MAT_AUG_OPERATOR1(op, cvop, const A, B)
00238 
00239 CV_MAT_AUG_OPERATOR  (+=, cv::add(a,b,a), Mat, Mat)
00240 CV_MAT_AUG_OPERATOR  (+=, cv::add(a,b,a), Mat, Scalar)
00241 CV_MAT_AUG_OPERATOR_T(+=, cv::add(a,b,a), Mat_<_Tp>, Mat)
00242 CV_MAT_AUG_OPERATOR_T(+=, cv::add(a,b,a), Mat_<_Tp>, Scalar)
00243 CV_MAT_AUG_OPERATOR_T(+=, cv::add(a,b,a), Mat_<_Tp>, Mat_<_Tp>)
00244 
00245 CV_MAT_AUG_OPERATOR  (-=, cv::subtract(a,b,a), Mat, Mat)
00246 CV_MAT_AUG_OPERATOR  (-=, cv::subtract(a,b,a), Mat, Scalar)
00247 CV_MAT_AUG_OPERATOR_T(-=, cv::subtract(a,b,a), Mat_<_Tp>, Mat)
00248 CV_MAT_AUG_OPERATOR_T(-=, cv::subtract(a,b,a), Mat_<_Tp>, Scalar)
00249 CV_MAT_AUG_OPERATOR_T(-=, cv::subtract(a,b,a), Mat_<_Tp>, Mat_<_Tp>)
00250 
00251 CV_MAT_AUG_OPERATOR  (*=, cv::gemm(a, b, 1, Mat(), 0, a, 0), Mat, Mat)
00252 CV_MAT_AUG_OPERATOR_T(*=, cv::gemm(a, b, 1, Mat(), 0, a, 0), Mat_<_Tp>, Mat)
00253 CV_MAT_AUG_OPERATOR_T(*=, cv::gemm(a, b, 1, Mat(), 0, a, 0), Mat_<_Tp>, Mat_<_Tp>)
00254 CV_MAT_AUG_OPERATOR  (*=, a.convertTo(a, -1, b), Mat, double)
00255 CV_MAT_AUG_OPERATOR_T(*=, a.convertTo(a, -1, b), Mat_<_Tp>, double)
00256 
00257 CV_MAT_AUG_OPERATOR  (/=, cv::divide(a,b,a), Mat, Mat)
00258 CV_MAT_AUG_OPERATOR_T(/=, cv::divide(a,b,a), Mat_<_Tp>, Mat)
00259 CV_MAT_AUG_OPERATOR_T(/=, cv::divide(a,b,a), Mat_<_Tp>, Mat_<_Tp>)
00260 CV_MAT_AUG_OPERATOR  (/=, a.convertTo((Mat&)a, -1, 1./b), Mat, double)
00261 CV_MAT_AUG_OPERATOR_T(/=, a.convertTo((Mat&)a, -1, 1./b), Mat_<_Tp>, double)
00262 
00263 CV_MAT_AUG_OPERATOR  (&=, cv::bitwise_and(a,b,a), Mat, Mat)
00264 CV_MAT_AUG_OPERATOR  (&=, cv::bitwise_and(a,b,a), Mat, Scalar)
00265 CV_MAT_AUG_OPERATOR_T(&=, cv::bitwise_and(a,b,a), Mat_<_Tp>, Mat)
00266 CV_MAT_AUG_OPERATOR_T(&=, cv::bitwise_and(a,b,a), Mat_<_Tp>, Scalar)
00267 CV_MAT_AUG_OPERATOR_T(&=, cv::bitwise_and(a,b,a), Mat_<_Tp>, Mat_<_Tp>)
00268 
00269 CV_MAT_AUG_OPERATOR  (|=, cv::bitwise_or(a,b,a), Mat, Mat)
00270 CV_MAT_AUG_OPERATOR  (|=, cv::bitwise_or(a,b,a), Mat, Scalar)
00271 CV_MAT_AUG_OPERATOR_T(|=, cv::bitwise_or(a,b,a), Mat_<_Tp>, Mat)
00272 CV_MAT_AUG_OPERATOR_T(|=, cv::bitwise_or(a,b,a), Mat_<_Tp>, Scalar)
00273 CV_MAT_AUG_OPERATOR_T(|=, cv::bitwise_or(a,b,a), Mat_<_Tp>, Mat_<_Tp>)
00274 
00275 CV_MAT_AUG_OPERATOR  (^=, cv::bitwise_xor(a,b,a), Mat, Mat)
00276 CV_MAT_AUG_OPERATOR  (^=, cv::bitwise_xor(a,b,a), Mat, Scalar)
00277 CV_MAT_AUG_OPERATOR_T(^=, cv::bitwise_xor(a,b,a), Mat_<_Tp>, Mat)
00278 CV_MAT_AUG_OPERATOR_T(^=, cv::bitwise_xor(a,b,a), Mat_<_Tp>, Scalar)
00279 CV_MAT_AUG_OPERATOR_T(^=, cv::bitwise_xor(a,b,a), Mat_<_Tp>, Mat_<_Tp>)
00280 
00281 #undef CV_MAT_AUG_OPERATOR_T
00282 #undef CV_MAT_AUG_OPERATOR
00283 #undef CV_MAT_AUG_OPERATOR1
00284 
00285 
00286 
00287 ///////////////////////////////////////////// SVD /////////////////////////////////////////////
00288 
00289 inline SVD::SVD() {}
00290 inline SVD::SVD( InputArray m, int flags ) { operator ()(m, flags); }
00291 inline void SVD::solveZ( InputArray m, OutputArray _dst )
00292 {
00293     Mat mtx = m.getMat();
00294     SVD svd(mtx, (mtx.rows >= mtx.cols ? 0 : SVD::FULL_UV));
00295     _dst.create(svd.vt.cols, 1, svd.vt.type());
00296     Mat dst = _dst.getMat();
00297     svd.vt.row(svd.vt.rows-1).reshape(1,svd.vt.cols).copyTo(dst);
00298 }
00299 
00300 template<typename _Tp, int m, int n, int nm> inline void
00301     SVD::compute( const Matx<_Tp, m, n>& a, Matx<_Tp, nm, 1>& w, Matx<_Tp, m, nm>& u, Matx<_Tp, n, nm>& vt )
00302 {
00303     CV_StaticAssert( nm == MIN(m, n), "Invalid size of output vector.");
00304     Mat _a(a, false), _u(u, false), _w(w, false), _vt(vt, false);
00305     SVD::compute(_a, _w, _u, _vt);
00306     CV_Assert(_w.data == (uchar*)&w.val[0] && _u.data == (uchar*)&u.val[0] && _vt.data == (uchar*)&vt.val[0]);
00307 }
00308 
00309 template<typename _Tp, int m, int n, int nm> inline void
00310 SVD::compute( const Matx<_Tp, m, n>& a, Matx<_Tp, nm, 1>& w )
00311 {
00312     CV_StaticAssert( nm == MIN(m, n), "Invalid size of output vector.");
00313     Mat _a(a, false), _w(w, false);
00314     SVD::compute(_a, _w);
00315     CV_Assert(_w.data == (uchar*)&w.val[0]);
00316 }
00317 
00318 template<typename _Tp, int m, int n, int nm, int nb> inline void
00319 SVD::backSubst( const Matx<_Tp, nm, 1>& w, const Matx<_Tp, m, nm>& u,
00320                 const Matx<_Tp, n, nm>& vt, const Matx<_Tp, m, nb>& rhs,
00321                 Matx<_Tp, n, nb>& dst )
00322 {
00323     CV_StaticAssert( nm == MIN(m, n), "Invalid size of output vector.");
00324     Mat _u(u, false), _w(w, false), _vt(vt, false), _rhs(rhs, false), _dst(dst, false);
00325     SVD::backSubst(_w, _u, _vt, _rhs, _dst);
00326     CV_Assert(_dst.data == (uchar*)&dst.val[0]);
00327 }
00328 
00329 
00330 
00331 /////////////////////////////////// Multiply-with-Carry RNG ///////////////////////////////////
00332 
00333 inline RNG::RNG()              { state = 0xffffffff; }
00334 inline RNG::RNG(uint64 _state) { state = _state ? _state : 0xffffffff; }
00335 
00336 inline RNG::operator uchar()    { return (uchar)next(); }
00337 inline RNG::operator schar ()    { return (schar)next(); }
00338 inline RNG::operator ushort ()   { return (ushort)next(); }
00339 inline RNG::operator short ()    { return (short)next(); }
00340 inline RNG::operator int ()      { return (int)next(); }
00341 inline RNG::operator unsigned () { return next(); }
00342 inline RNG::operator float ()    { return next()*2.3283064365386962890625e-10f; }
00343 inline RNG::operator double ()   { unsigned t = next(); return (((uint64)t << 32) | next()) * 5.4210108624275221700372640043497e-20; }
00344 
00345 inline unsigned RNG::operator ()(unsigned N) { return (unsigned)uniform(0,N); }
00346 inline unsigned RNG::operator ()()           { return next(); }
00347 
00348 inline int    RNG::uniform(int a, int b)       { return a == b ? a : (int)(next() % (b - a) + a); }
00349 inline float  RNG::uniform(float a, float b)   { return ((float)*this)*(b - a) + a; }
00350 inline double RNG::uniform(double a, double b) { return ((double)*this)*(b - a) + a; }
00351 
00352 inline unsigned RNG::next()
00353 {
00354     state = (uint64)(unsigned)state* /*CV_RNG_COEFF*/ 4164903690U + (unsigned)(state >> 32);
00355     return (unsigned)state;
00356 }
00357 
00358 //! returns the next unifomly-distributed random number of the specified type
00359 template<typename _Tp> static inline _Tp randu()
00360 {
00361   return (_Tp)theRNG();
00362 }
00363 
00364 ///////////////////////////////// Formatted string generation /////////////////////////////////
00365 
00366 CV_EXPORTS String format( const char* fmt, ... );
00367 
00368 ///////////////////////////////// Formatted output of cv::Mat /////////////////////////////////
00369 
00370 static inline
00371 Ptr<Formatted> format(InputArray mtx, int fmt)
00372 {
00373     return Formatter::get(fmt)->format(mtx.getMat());
00374 }
00375 
00376 static inline
00377 int print(Ptr<Formatted> fmtd, FILE* stream = stdout)
00378 {
00379     int written = 0;
00380     fmtd->reset();
00381     for(const char* str = fmtd->next(); str; str = fmtd->next())
00382         written += fputs(str, stream);
00383 
00384     return written;
00385 }
00386 
00387 static inline
00388 int print(const Mat& mtx, FILE* stream = stdout)
00389 {
00390     return print(Formatter::get()->format(mtx), stream);
00391 }
00392 
00393 static inline
00394 int print(const UMat& mtx, FILE* stream = stdout)
00395 {
00396     return print(Formatter::get()->format(mtx.getMat(ACCESS_READ)), stream);
00397 }
00398 
00399 template<typename _Tp> static inline
00400 int print(const std::vector<Point_<_Tp> >& vec, FILE* stream = stdout)
00401 {
00402     return print(Formatter::get()->format(Mat(vec)), stream);
00403 }
00404 
00405 template<typename _Tp> static inline
00406 int print(const std::vector<Point3_<_Tp> >& vec, FILE* stream = stdout)
00407 {
00408     return print(Formatter::get()->format(Mat(vec)), stream);
00409 }
00410 
00411 template<typename _Tp, int m, int n> static inline
00412 int print(const Matx<_Tp, m, n>& matx, FILE* stream = stdout)
00413 {
00414     return print(Formatter::get()->format(cv::Mat(matx)), stream);
00415 }
00416 
00417 //! @endcond
00418 
00419 /****************************************************************************************\
00420 *                                  Auxiliary algorithms                                  *
00421 \****************************************************************************************/
00422 
00423 /** @brief Splits an element set into equivalency classes.
00424 
00425 The generic function partition implements an \f$O(N^2)\f$ algorithm for splitting a set of \f$N\f$ elements
00426 into one or more equivalency classes, as described in
00427 <http://en.wikipedia.org/wiki/Disjoint-set_data_structure> . The function returns the number of
00428 equivalency classes.
00429 @param _vec Set of elements stored as a vector.
00430 @param labels Output vector of labels. It contains as many elements as vec. Each label labels[i] is
00431 a 0-based cluster index of `vec[i]`.
00432 @param predicate Equivalence predicate (pointer to a boolean function of two arguments or an
00433 instance of the class that has the method bool operator()(const _Tp& a, const _Tp& b) ). The
00434 predicate returns true when the elements are certainly in the same class, and returns false if they
00435 may or may not be in the same class.
00436 @ingroup core_cluster
00437 */
00438 template<typename _Tp, class _EqPredicate> int
00439 partition( const std::vector<_Tp>& _vec, std::vector<int>& labels,
00440           _EqPredicate predicate=_EqPredicate())
00441 {
00442     int i, j, N = (int)_vec.size();
00443     const _Tp* vec = &_vec[0];
00444 
00445     const int PARENT=0;
00446     const int RANK=1;
00447 
00448     std::vector<int> _nodes(N*2);
00449     int (*nodes)[2] = (int(*)[2])&_nodes[0];
00450 
00451     // The first O(N) pass: create N single-vertex trees
00452     for(i = 0; i < N; i++)
00453     {
00454         nodes[i][PARENT]=-1;
00455         nodes[i][RANK] = 0;
00456     }
00457 
00458     // The main O(N^2) pass: merge connected components
00459     for( i = 0; i < N; i++ )
00460     {
00461         int root = i;
00462 
00463         // find root
00464         while( nodes[root][PARENT] >= 0 )
00465             root = nodes[root][PARENT];
00466 
00467         for( j = 0; j < N; j++ )
00468         {
00469             if( i == j || !predicate(vec[i], vec[j]))
00470                 continue;
00471             int root2 = j;
00472 
00473             while( nodes[root2][PARENT] >= 0 )
00474                 root2 = nodes[root2][PARENT];
00475 
00476             if( root2 != root )
00477             {
00478                 // unite both trees
00479                 int rank = nodes[root][RANK], rank2 = nodes[root2][RANK];
00480                 if( rank > rank2 )
00481                     nodes[root2][PARENT] = root;
00482                 else
00483                 {
00484                     nodes[root][PARENT] = root2;
00485                     nodes[root2][RANK] += rank == rank2;
00486                     root = root2;
00487                 }
00488                 CV_Assert( nodes[root][PARENT] < 0 );
00489 
00490                 int k = j, parent;
00491 
00492                 // compress the path from node2 to root
00493                 while( (parent = nodes[k][PARENT]) >= 0 )
00494                 {
00495                     nodes[k][PARENT] = root;
00496                     k = parent;
00497                 }
00498 
00499                 // compress the path from node to root
00500                 k = i;
00501                 while( (parent = nodes[k][PARENT]) >= 0 )
00502                 {
00503                     nodes[k][PARENT] = root;
00504                     k = parent;
00505                 }
00506             }
00507         }
00508     }
00509 
00510     // Final O(N) pass: enumerate classes
00511     labels.resize(N);
00512     int nclasses = 0;
00513 
00514     for( i = 0; i < N; i++ )
00515     {
00516         int root = i;
00517         while( nodes[root][PARENT] >= 0 )
00518             root = nodes[root][PARENT];
00519         // re-use the rank as the class label
00520         if( nodes[root][RANK] >= 0 )
00521             nodes[root][RANK] = ~nclasses++;
00522         labels[i] = ~nodes[root][RANK];
00523     }
00524 
00525     return nclasses;
00526 }
00527 
00528 } // cv
00529 
00530 #endif
00531