core/vnl/algo/vnl_rnpoly_solve.h
Go to the documentation of this file.
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_