00001 // This is core/vnl/algo/vnl_rnpoly_solve.h 00002 #ifndef vnl_rnpoly_solve_h_ 00003 #define vnl_rnpoly_solve_h_ 00004 #ifdef VCL_NEEDS_PRAGMA_INTERFACE 00005 #pragma interface 00006 #endif 00007 //: 00008 // \file 00009 // \brief Solves for roots of system of real polynomials 00010 // \author Marc Pollefeys, ESAT-VISICS, K.U.Leuven 00011 // \date 12-Aug-1997 00012 // 00013 // \verbatim 00014 // Modifications 00015 // Oct.1999 - Peter Vanroose - implementation simplified through "cmplx" class for doing complex arithmetic. 00016 // May.2002 - Peter Vanroose - added operator*=(cmplx) and operator/=(cmplx) 00017 // Mar.2003 - Peter Vanroose - renamed M to M_, T to T_ 00018 // Feb.2004 - Peter Vanroose - removed hard limits on dimensionality; this gets rid of M_ and T_; 00019 // now using std::vector throughout instead of C arrays of fixed size 00020 // \endverbatim 00021 00022 #include <vnl/vnl_vector.h> 00023 #include <vnl/vnl_real_npolynomial.h> 00024 #include <vcl_vector.h> 00025 00026 //: Solves for roots of system of real polynomials 00027 // Calculates all the roots of a system of N polynomials in N variables 00028 // through continuation. 00029 // Adapted from the PARALLEL CONTINUATION algorithm, written by Darrell 00030 // Stam, 1991, and further improved by Kriegman and Ponce, 1992. 00031 00032 class vnl_rnpoly_solve 00033 { 00034 // Data Members-------------------------------------------------------------- 00035 vcl_vector<vnl_real_npolynomial*> ps_; // the input 00036 vcl_vector<vnl_vector<double>*> r_; // the output (real part) 00037 vcl_vector<vnl_vector<double>*> i_; // the output (imaginary part) 00038 00039 public: 00040 00041 // Constructor--------------------------------------------------------------- 00042 00043 //: The constructor already does all the calculations 00044 inline vnl_rnpoly_solve(vcl_vector<vnl_real_npolynomial*> const& ps) 00045 : ps_(ps) { compute(); } 00046 00047 // Destructor---------------------------------------------------------------- 00048 00049 ~vnl_rnpoly_solve(); 00050 00051 // Operations---------------------------------------------------------------- 00052 00053 //: Array of real parts of roots 00054 inline vcl_vector<vnl_vector<double>*> real() { return r_; } 00055 00056 //: Array of imaginary parts of roots 00057 inline vcl_vector<vnl_vector<double>*> imag() { return i_; } 00058 00059 //: Return real roots only. 00060 // Roots are real if the absolute value of their imaginary part is less than 00061 // the optional argument tol, which defaults to 1e-12 [untested] 00062 vcl_vector<vnl_vector<double>*> realroots(double tol = 1e-12); 00063 00064 // Computations-------------------------------------------------------------- 00065 00066 private: 00067 //: Compute roots using continuation algorithm. 00068 bool compute(); 00069 00070 void Read_Input(vcl_vector<unsigned int>& ideg, 00071 vcl_vector<unsigned int>& terms, 00072 vcl_vector<int>& polyn, 00073 vcl_vector<double>& coeff); 00074 }; 00075 00076 #endif // vnl_rnpoly_solve_h_