latest version v1.9 - last update 10 Apr 2010 |
00001 /* 00002 * Copyright (C) 2002, 2003, 2004, 2005, 2006 00003 * Lehrstuhl fuer Technische Informatik, RWTH-Aachen, Germany 00004 * 00005 * This file is part of the LTI-Computer Vision Library (LTI-Lib) 00006 * 00007 * The LTI-Lib is free software; you can redistribute it and/or 00008 * modify it under the terms of the GNU Lesser General Public License (LGPL) 00009 * as published by the Free Software Foundation; either version 2.1 of 00010 * the License, or (at your option) any later version. 00011 * 00012 * The LTI-Lib is distributed in the hope that it will be 00013 * useful, but WITHOUT ANY WARRANTY; without even the implied warranty 00014 * of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00015 * GNU Lesser General Public License for more details. 00016 * 00017 * You should have received a copy of the GNU Lesser General Public 00018 * License along with the LTI-Lib; see the file LICENSE. If 00019 * not, write to the Free Software Foundation, Inc., 59 Temple Place - 00020 * Suite 330, Boston, MA 02111-1307, USA. 00021 */ 00022 00023 00024 /*-------------------------------------------------------------------- 00025 * project ....: LTI-Lib: Image Processing and Computer Vision Library 00026 * file .......: ltiBlobEM.h 00027 * authors ....: Suat Akyol 00028 * organization: LTI, RWTH Aachen 00029 * creation ...: 27.5.2002 00030 * revisions ..: $Id: ltiBlobEM.h,v 1.7 2006/02/07 18:31:33 ltilib Exp $ 00031 */ 00032 00033 #ifndef _LTI_BLOB_EM_H_ 00034 #define _LTI_BLOB_EM_H_ 00035 00036 // include only those files which are needed in this header!! 00037 #include <vector> 00038 00039 #include "ltiImage.h" 00040 #include "ltiVector.h" 00041 #include "ltiMatrix.h" 00042 00043 // include to parent class: 00044 #include "ltiFunctor.h" 00045 00046 namespace lti { 00047 00048 /** 00049 * Estimates the position of M overlapping blobs by applying the 00050 * EM-algorithm and estimating the parameters of a gaussian mixture 00051 * model that fits the blobs. As input a distance transformed 00052 * image of the blobs or similar is expected and a vector of ellipses 00053 * for representing a gaussian distribution component. 00054 * 00055 * For a description of the EM-algorithm see Jeff Bilmes "A Gentle 00056 * Tutorial of the EM Algorithm and its Application to Parameter 00057 * Estimation for Gaussian Mixture and Hidden Markov Models", 00058 * Technical Report TR-97-021 by the International Computer 00059 * Science Institute, Berkeley, Ca. 00060 * 00061 * \htmlonly<a name="HINTS"></a>\endhtmlonly 00062 * 00063 * For usage of this functor please see the following examples and 00064 * hints: 00065 * 00066 * \li <a href="#ELLIPSES">Preparing the input ellipses</a> 00067 * \li <a href="#APPLYING">Applying this functor</a> 00068 * \li <a href="#ERRORS">Possible errors on usage</a> 00069 * 00070 * \par Example for preparing ellipses 00071 * \htmlonly 00072 * <a name="ELLIPSES"></a> 00073 * <a href="#HINTS">back</a> 00074 * \endhtmlonly 00075 * 00076 * \code 00077 * 00078 * // Typedef for better readability 00079 * typedef lti::blobEM::gaussEllipse TEllipse; 00080 * 00081 * // Vector of ellipses 00082 * std::vector<TEllipse>& allEllipses; 00083 * 00084 * // Create one ellipse with center(10,10), eigenvalues 100 and 20 00085 * // and 45 degrees orientation 00086 * TEllipse oneEllipse(lti::tpoint<double>(10,10),100,20,lti::Pi/4); 00087 * 00088 * // Constrain ellipse properties 00089 * oneEllipse.constrainShape = true; 00090 * 00091 * oneEllipse.constrainArea = true; 00092 * oneEllipse.areaTolerance = 0.2; 00093 * 00094 * oneEllipse.constrainAngle = true; 00095 * oneEllipse.angleTolerance = lti::Pi/8; 00096 * 00097 * // Create another ellipse with the same properties as the first 00098 * // but center(50,50). This ellipse will be unconstrained (default) 00099 * TEllipse anotherEllipse; 00100 * 00101 * anotherEllipse.fromEllipse(oneEllipse); // is not copy operator! 00102 * // copies only properties 00103 * // angle, center, lambda1,2 00104 * anotherEllipse.center = lti::tpoint<double>(50,50); 00105 * 00106 * allEllipses.push_back(oneEllipse); 00107 * allEllipses.push_back(anotherEllipse); 00108 * 00109 * \endcode 00110 * 00111 * 00112 * \par Example for using this functor 00113 * \htmlonly 00114 * <a name="APPLYING"></a> 00115 * <a href="#HINTS">back</a> 00116 * \endhtmlonly 00117 * 00118 * \code 00119 * 00120 * // prepare functor 00121 * lti::blobEM em; 00122 * lti::blobEM::parameters emParam; 00123 * emParam.maxIterations = 10; // number of maximum allowed iterations 00124 * emParam.convergenceThreshold = -1; // force maximum iterations 00125 * em.setParameters(emParam); 00126 * 00127 * ////////////////////////////////////////////////////////////////////////// 00128 * // Apply with allEllipses but make sure that initialization is good, 00129 * // since result stongly depends on initialization. 00130 * // Result will be delivered back in allEllipses. aChnl8 should be a 00131 * // distance transform of a binary image (see lti::distanceTransform) 00132 * em.apply(aChnl8,allEllipses); 00133 * 00134 * ////////////////////////////////////////////////////////////////////////// 00135 * // Alternatively you can apply without knowing how to initialize by 00136 * // merely specifying the number of ellipses (here 2). These will be 00137 * // set randomly then. 00138 * // You must explicitly request the resulting ellipses afterwards 00139 * // by calling computeEllipses(...) 00140 * em.apply(aChnl8,2); 00141 * em.computeEllipses(allEllipses); 00142 * 00143 * ////////////////////////////////////////////////////////////////////////// 00144 * // Another possibility is to control the iterations manually 00145 * em.initialize(aChnl8,allEllipses); 00146 * // OR 00147 * em.initialize(aChnl8,2); 00148 * 00149 * // perfom iterations 00150 * em.iterate(); 00151 * em.iterate(); 00152 * while (em.getQDiff()>0) { 00153 * // the Q value is maximized by EM. QDiff will be > 0 00154 * // As long as Q increases. To get valid QDiff, we have 00155 * // to perform two iterations before this loop. 00156 * em.iterate(); 00157 * } 00158 * 00159 * em.computeEllipses(allEllipses); 00160 * 00161 * // If you are interested in the internal state of the functor use 00162 * // the getXXX(...) methods, which deliver const references (i.e. 00163 * // you don't have to copy them, if you only want to see them) 00164 * 00165 * std::vector< lti::matrix<double> >& covariances = em.getCovariances(); 00166 * std::vector< lti::vector<double> >& centers = em.getCenters(); 00167 * std::vector< double >& alphas = em.getAlphas(); 00168 * 00169 * // For a description of their meaning refer to the given bibliography 00170 * 00171 * \endcode 00172 * 00173 * 00174 * \par Beware of possible errors 00175 * \htmlonly 00176 * <a name="ERRORS"></a> 00177 * <a href="#HINTS">back</a> 00178 * \endhtmlonly 00179 * 00180 * \li <b>Constraints are relative to given ellipse</b>. If you use the 00181 * returned ellipses for initializing another run of the functor, 00182 * then please consider that their properties could have been changed 00183 * by the first run. 00184 * Example: You have an ellipse that is supposed to stay upright 00185 * with a tolerance of 5°. You call apply once and the ellipse is indeed 00186 * rotated 5°. You use the output for the next call of apply. And again 00187 * the ellipse is rotated 5°, resulting in a total of 10° which is 00188 * correct but possibly not what you expected. 00189 * 00190 */ 00191 class blobEM : public functor { 00192 00193 public: 00194 /** 00195 * the parameters for the class blobEM 00196 */ 00197 class parameters : public functor::parameters { 00198 public: 00199 /** 00200 * default constructor 00201 */ 00202 parameters(); 00203 00204 /** 00205 * copy constructor 00206 * @param other the parameters object to be copied 00207 */ 00208 parameters(const parameters& other); 00209 00210 /** 00211 * destructor 00212 */ 00213 ~parameters(); 00214 00215 /** 00216 * returns name of this type 00217 */ 00218 const char* getTypeName() const; 00219 00220 /** 00221 * copy the contents of a parameters object 00222 * @param other the parameters object to be copied 00223 * @return a reference to this parameters object 00224 */ 00225 parameters& copy(const parameters& other); 00226 00227 /** 00228 * copy the contents of a parameters object 00229 * @param other the parameters object to be copied 00230 * @return a reference to this parameters object 00231 */ 00232 parameters& operator=(const parameters& other); 00233 00234 00235 /** 00236 * returns a pointer to a clone of the parameters 00237 */ 00238 virtual functor::parameters* clone() const; 00239 00240 /** 00241 * write the parameters in the given ioHandler 00242 * @param handler the ioHandler to be used 00243 * @param complete if true (the default) the enclosing begin/end will 00244 * be also written, otherwise only the data block will be written. 00245 * @return true if write was successful 00246 */ 00247 virtual bool write(ioHandler& handler,const bool complete=true) const; 00248 00249 /** 00250 * read the parameters from the given ioHandler 00251 * @param handler the ioHandler to be used 00252 * @param complete if true (the default) the enclosing begin/end will 00253 * be also written, otherwise only the data block will be written. 00254 * @return true if write was successful 00255 */ 00256 virtual bool read(ioHandler& handler,const bool complete=true); 00257 00258 # ifdef _LTI_MSC_6 00259 /** 00260 * this function is required by MSVC only, as a workaround for a 00261 * very awful bug, which exists since MSVC V.4.0, and still by 00262 * V.6.0 with all bugfixes (so called "service packs") remains 00263 * there... This method is also public due to another bug, so please 00264 * NEVER EVER call this method directly: use read() instead 00265 */ 00266 bool readMS(ioHandler& handler,const bool complete=true); 00267 00268 /** 00269 * this function is required by MSVC only, as a workaround for a 00270 * very awful bug, which exists since MSVC V.4.0, and still by 00271 * V.6.0 with all bugfixes (so called "service packs") remains 00272 * there... This method is also public due to another bug, so please 00273 * NEVER EVER call this method directly: use write() instead 00274 */ 00275 bool writeMS(ioHandler& handler,const bool complete=true) const; 00276 # endif 00277 00278 // ------------------------------------------------ 00279 // the parameters 00280 // ------------------------------------------------ 00281 00282 // If you add more parameters manually, do not forget to do following: 00283 // 1. indicate in the default constructor the default values 00284 // 2. make sure that the copy member also copy your new parameters 00285 // 3. make sure that the read and write members also read and 00286 // write your parameters 00287 00288 /** Maximum number of iterations. Default = 100 00289 */ 00290 int maxIterations; 00291 00292 /** If the difference of the Q value between two successive 00293 * iterations (relative to the last Q value) falls below 00294 * this threshold, then the iterations are aborted. 00295 * A value < zero enforces the maximum number of iterations. 00296 * Default 1e-6 00297 */ 00298 double convergenceThreshold; 00299 00300 }; 00301 00302 00303 /** 00304 * An internal class of lti::blobEM for handling 2D gaussian ellipses. 00305 * The properties can be set manually. 00306 * 00307 * Another possibility is to obtain the properties from 00308 * a 2x2 covariance matrix (except center) or another 00309 * ellipse, either considering constraints or not. 00310 */ 00311 class gaussEllipse { 00312 public: 00313 00314 gaussEllipse(); 00315 00316 gaussEllipse(lti::tpoint<double> cen, 00317 double l1, 00318 double l2, 00319 double ang); 00320 00321 ~gaussEllipse(); 00322 00323 00324 /** @name Ellipse properties 00325 */ 00326 //@{ 00327 00328 /** Center of ellipse 00329 */ 00330 lti::tpoint<double> center; 00331 00332 /** Variance along main axis (eigenvalue 1). 00333 * Ensure lambda1 to be > lambda2 for plausibility. 00334 * 00335 * \b Note: This is not the standard deviation! 00336 * For visualization of an ellipse use std dev. 00337 */ 00338 double lambda1; 00339 00340 /** Variance along second axis (eigenvalue 2) 00341 * 00342 * \b Note: This is not the standard deviation! 00343 * For visualization of an ellipse use std dev. 00344 */ 00345 double lambda2; 00346 00347 /** Angle of main axis (-Pi/2 to +Pi/2). Second axis is 00348 * perpendicular to main axis. 00349 */ 00350 double angle; 00351 00352 //@} 00353 00354 00355 /** @name Ellipse constraints 00356 * Use these to indicate constraints for the properties of an ellipse. 00357 * \b NOTE: constraining properties can make it impossible 00358 * for the EM-Algorithm to converge. Use with 00359 * caution, preferably only for one component 00360 * amongst several! 00361 */ 00362 //@{ 00363 00364 /** keep center within centerTolerance (default = false) 00365 */ 00366 bool constrainCenter; 00367 00368 /** keep shape (lambda1/lambda2) within shapeTolerance 00369 * (default = false) 00370 */ 00371 bool constrainShape; 00372 00373 /** keep area size (lambda1*lambda2) within areaTolerance (default = false) 00374 */ 00375 bool constrainArea; 00376 00377 /** keep orientation angle within angleTolerance (default = false) 00378 */ 00379 bool constrainAngle; 00380 00381 //@} 00382 00383 00384 /** @name Constraint tolerances 00385 * Setting these to values > 0 (default) relaxes the 00386 * constraint conditions. 00387 */ 00388 //@{ 00389 00390 /** The center tolerance as maximum Euclidean distance in pixel. 00391 * (default = 0) 00392 * 00393 * \f$ dist(newCenter - center) < centerTolerance \f$ 00394 */ 00395 double centerTolerance; 00396 00397 /** The relative shape tolerance (default = 0). 00398 * 00399 * <b> NOT CONSIDERED YET!!! Development in Progress! </b> 00400 */ 00401 double shapeTolerance; 00402 00403 /** The relative size tolerance (default = 0), which means 00404 * Ellipse area may change in range 00405 * 00406 * \f$ area/a <= newArea <= area*a \f$ with \f$ a = (1+areaTolerance)\f$ 00407 */ 00408 double areaTolerance; 00409 00410 /** The (+/-) angle tolerance (default = 0). Reasonable range 00411 * is from 0 to Pi/2, since this is Pi/2 cyclic. 00412 * 00413 * \f$ angle - angleTolerance <= newAngle <= angle + angleTolerance \f$ 00414 */ 00415 double angleTolerance; 00416 00417 //@} 00418 00419 00420 /** Get Ellipse properties lambda and angle 00421 * from a 2x2 covariance matrix considering 00422 * shape, size, and/or angle constraints, if 00423 * desired. 00424 * 00425 * \b Note: center can NOT be obtained from 00426 * covariance. 00427 * 00428 * Returns false, if not possible 00429 */ 00430 bool from2x2Covariance(const lti::matrix<double>& cov, 00431 const bool& constraints = true); 00432 00433 /** Get all ellipse properties (center, lambda, angle) 00434 * from another ellipse considering center, shape, size, 00435 * and/or angle constraints, if desired. 00436 * THIS IS NOT A COPY METHOD!!! 00437 * 00438 * Returns false, if not possible. 00439 */ 00440 bool fromEllipse(const gaussEllipse& other, 00441 const bool& constraints = true); 00442 00443 /** Generate 2x2 covariance matrix with ellipse params 00444 * lambda and angle. Returns false, if not possible 00445 */ 00446 bool to2x2Covariance(lti::matrix<double>& cov) const; 00447 00448 private: 00449 00450 /** Perform pca on 2x2 covariance matrix. 00451 * Lambdas are the eigenvalues and e are the 00452 * eigenvectors (ordered by magnitude). 00453 * Returns false, if not possible. 00454 */ 00455 static bool pca2x2(const lti::matrix<double>& cov, 00456 double& lambda1, 00457 double& lambda2, 00458 lti::vector<double>& e1, 00459 lti::vector<double>& e2); 00460 }; 00461 00462 /** 00463 * default constructor 00464 */ 00465 blobEM(); 00466 00467 /** 00468 * copy constructor 00469 * @param other the object to be copied 00470 */ 00471 blobEM(const blobEM& other); 00472 00473 /** 00474 * destructor 00475 */ 00476 virtual ~blobEM(); 00477 00478 /** 00479 * returns the name of this type ("blobEM") 00480 */ 00481 virtual const char* getTypeName() const; 00482 00483 /** 00484 * Calls initialize() once, then iterate() several times until 00485 * maximum number of iterations or abortion criterion reached 00486 * (percental increase of Q value below given threshold) 00487 * @param src channel8 with the source data. Should be distance 00488 * transform of a binary image. 00489 * @param elem vector of blobEM::gaussEllipse with initial distribution 00490 * parameters, represented by 2D-ellipses. Result of EM 00491 * iteration is also left here. 00492 * @return true if apply successful or false otherwise. 00493 */ 00494 bool apply(const channel8& src, std::vector<gaussEllipse>& elem); 00495 00496 /** 00497 * Calls initialize() once, then iterate() several times until 00498 * maximum number of iterations or abortion criterion reached 00499 * (percental increase of Q value below given threshold) 00500 * @param src channel8 with the source data. Should be distance 00501 * transform of a binary image. 00502 * @param components number of mixture components. Will be initialized 00503 * randomly. For deterministic initialization use other 00504 * apply method. 00505 * @return true if apply successful or false otherwise. 00506 */ 00507 bool apply(const channel8& src, const int& components); 00508 00509 /** 00510 * copy data of "other" functor. 00511 * @param other the functor to be copied 00512 * @return a reference to this functor object 00513 */ 00514 blobEM& copy(const blobEM& other); 00515 00516 /** 00517 * alias for copy member 00518 * @param other the functor to be copied 00519 * @return a reference to this functor object 00520 */ 00521 blobEM& operator=(const blobEM& other); 00522 00523 /** 00524 * returns a pointer to a clone of this functor. 00525 */ 00526 virtual functor* clone() const; 00527 00528 /** 00529 * returns used parameters 00530 */ 00531 const parameters& getParameters() const; 00532 00533 // If you add more attributes manually, do not forget to do following: 00534 // 1. indicate in the default constructor the default values 00535 // 2. make sure that the copy member also copy your new attributes, or 00536 // to ensure there, that these attributes are properly initialized. 00537 00538 /** @name Manual control 00539 * These methods are required for the manual control of the EM 00540 */ 00541 //@{ 00542 00543 /** Explicitly initialize this functor, for example if you intend to 00544 * perform each iteration singly. Initializes with the given ellipses. 00545 */ 00546 bool initialize(const channel8& src, const std::vector<gaussEllipse>& elem); 00547 00548 /** Explicitly initialize this functor, for example if you intend to 00549 * perform each iteration singly. Initializes with M random ellipses, 00550 * with M = components 00551 */ 00552 bool initialize(const channel8& src, const int& components); 00553 00554 /** perform one single EM iteration. returns the total Q value, that has 00555 * to be maximized by EM. \b Note: Call initialize() first !!! 00556 */ 00557 double iterate(); 00558 00559 /** compute and return ellipses based on internal covariances and centers 00560 */ 00561 bool computeEllipses(std::vector<gaussEllipse>& ellipses); 00562 00563 /** reset this functors internal members and the working data 00564 */ 00565 void reset(); 00566 00567 //@} 00568 00569 /** @name State access 00570 * Use these to access the internal state of this functor. 00571 */ 00572 //@{ 00573 00574 /** Access the number of iterations so far 00575 */ 00576 const int& getIterations() const; 00577 00578 /** Access the number of components 00579 */ 00580 const int& getM() const; 00581 00582 /** Access the number of data elements considered 00583 */ 00584 const int& getN() const; 00585 00586 /** Access the vector of the M mixture component weights 00587 */ 00588 const std::vector<double>& getAlphas() const; 00589 00590 /** Access the centers of the M center components 00591 */ 00592 const std::vector<tpoint<double> >& getCenters() const; 00593 00594 /** Access the covariance matrices of the M center components 00595 */ 00596 const std::vector<matrix<double> >& getCovariances() const; 00597 00598 /** Access the Q value, that is maximized by EM 00599 */ 00600 const double& getQ() const; 00601 00602 /** Access the difference of current and last Q value 00603 */ 00604 const double& getQDiff() const; 00605 00606 //@} 00607 00608 private: 00609 /** @name Internal state 00610 * These are the internal state members. 00611 */ 00612 //@{ 00613 00614 /** The weights (a-priori probs) of the M mixture components 00615 */ 00616 std::vector<double> m_alphas; 00617 00618 /** The centers of the M mixture components 00619 */ 00620 std::vector< tpoint<double> > m_centers; 00621 00622 /** The covariance matrices of the M center components 00623 */ 00624 std::vector< matrix<double> > m_covariances; 00625 00626 //@} 00627 00628 /** @name Working data 00629 */ 00630 //@{ 00631 00632 /** The number of components 00633 */ 00634 int m_iterations; 00635 00636 /** The number of components 00637 */ 00638 int m_M; 00639 00640 /** The number of data elements (sum of Ni) 00641 */ 00642 int m_N; 00643 00644 /** The x coordinates 00645 */ 00646 std::vector<int> m_xi; 00647 00648 /** The y coordinates 00649 */ 00650 std::vector<int> m_yi; 00651 00652 /** The frequency of the i-th x,y pair 00653 */ 00654 std::vector<int> m_Ni; 00655 00656 /** The ellipses used for initializing 00657 */ 00658 std::vector<gaussEllipse> m_ellipses; 00659 00660 /** The Q value, that is maximized by EM 00661 */ 00662 double m_Q; 00663 00664 /** The difference of current and last Q value 00665 */ 00666 double m_QDiff; 00667 00668 //@} 00669 00670 00671 /** Performs a single EM step to estimate the new parameters, 00672 * and returns the likelihood 00673 */ 00674 double singleEMStep(); 00675 00676 /** 00677 * An internal helper class for 2D gaussian PDF p(x,y). 00678 * Use only with 2x2 covariance matrices (symmetric)! 00679 */ 00680 class gauss2DPDF { 00681 public: 00682 00683 gauss2DPDF() { 00684 reset(); 00685 } 00686 00687 ~gauss2DPDF() {} 00688 00689 /** evaluates the pdf with the given x,y pair and center and 00690 * the internally stored covariance matrix. 00691 */ 00692 double evaluate(int x, int y, const tpoint<double>& cen) const; 00693 00694 /** sets the covariance matrix, computes its determinant, 00695 * and its inverse. returns false and unsets covariance 00696 * and inverse, on error. 00697 */ 00698 bool set2x2CovarianceMatrix(const lti::matrix<double>& covar); 00699 00700 /** resets the internal state 00701 */ 00702 void reset(); 00703 00704 private: 00705 lti::matrix<double> m_covar; 00706 lti::matrix<double> m_covarInverse; 00707 double m_det; 00708 00709 static const lti::point m_validSize; 00710 }; 00711 00712 // covariance regularization constant (see Ueda, Nakano et al. 00713 // "SMEM Algorithm for Mixture Models", resp. Ormoneit & Tresp 00714 // "Improved Gaussian Mixture Density Estimates Using Bayesian 00715 // Penalty Terms and Network Averaging". 00716 static const double LAMBDA; 00717 00718 // regularization constant for alphas 00719 static const double ALPHAEPSILON; 00720 }; 00721 00722 /** 00723 * write a blobEM::gaussEllipse into the given stream handler 00724 * @param handler ioHandler where the ellipse should be written. 00725 * @param el ellipse 00726 * @param complete if true (default), begin and end tokens will be written 00727 * around the object. 00728 * @return true if successful, false otherwise 00729 */ 00730 bool write(ioHandler& handler, 00731 const blobEM::gaussEllipse& el, 00732 const bool complete=true); 00733 00734 /** 00735 * read a blobEM::gaussEllipse from the given stream handler 00736 * @param handler ioHandler where the ellipse should be read from. 00737 * @param el ellipse will be left here 00738 * @param complete if true (default), begin and end tokens will be read too. 00739 * @return true if successful, false otherwise 00740 */ 00741 bool read(ioHandler& handler, 00742 blobEM::gaussEllipse& el, 00743 const bool complete=true); 00744 } 00745 00746 #endif