vnl_rnpoly_solve.h

InsightToolkit-1.4.0(有大量的优化算法程序)

// This is core/vnl/algo/vnl_rnpoly_solve.h
#ifndef vnl_rnpoly_solve_h_
#define vnl_rnpoly_solve_h_
#ifdef VCL_NEEDS_PRAGMA_INTERFACE
#pragma interface
#endif
//:
// \file
// \brief Solves for roots of system of real polynomials
// \author Marc Pollefeys, ESAT-VISICS, K.U.Leuven
// \date   12-08-97
//
// \verbatim
// Modifications
//  Peter Vanroose, 20 Oct 1999: implementation simplified through "cmplx" class for doing complex arithmetic.
//  dac (Manchester) 28/03/2001: tidied up documentation
//   Feb.2002 - Peter Vanroose - brief doxygen comment placed on single line
//   May.2002 - Peter Vanroose - added operator*=(cmplx) and operator/=(cmplx)
// \endverbatim

#include <vnl/vnl_vector.h>
#include <vnl/vnl_real_npolynomial.h>
#include <vcl_vector.h>

#ifdef static
# error "grr!!"
#endif

//: Solves for roots of system of real polynomials
//  Calculates all the roots of a system of N polynomials in N variables
//  through continuation.
//  Adapted from the  PARALLEL CONTINUATION algorithm , written by Darrell
//  Stam, 1991, and further improved by  Kriegman and Ponce, 1992.

class vnl_rnpoly_solve
{
 public:
#ifndef VCL_WIN32
  static const unsigned int M = 11;   // Maximum dimension of problem
  static const unsigned int T = 2500; // Max. number of terms in a polynomial
#else
  enum { M = 11 };   // Maximum dimension of problem
  enum { T = 2500 }; // Maximum number of terms in a polynomial
#endif

  // Constructor---------------------------------------------------------------

  //: The constructor already does all the calculations
  inline vnl_rnpoly_solve(vcl_vector<vnl_real_npolynomial*> const& ps)
    : ps_(ps) { 
    this->vnl_rnpoly_solve::compute(); 
  }

  // Destructor----------------------------------------------------------------

 ~vnl_rnpoly_solve();

  // Operations----------------------------------------------------------------

  //: Array of real parts of roots
  inline vcl_vector<vnl_vector<double>*> real() { return r_; }

  //: Array of imaginary parts of roots
  inline vcl_vector<vnl_vector<double>*> imag() { return i_; }

  //: Return real roots only.
  //  Roots are real if the absolute value of their imaginary part is less than
  //  the optional argument tol, which defaults to 1e-12 [untested]
  vcl_vector<vnl_vector<double>*> realroots(double tol = 1e-12);

  // Computations--------------------------------------------------------------

 private:
  //: Compute roots using continuation algorithm.
public:
  bool compute();
  bool compute2();

  int Read_Input(int ideg[M], int terms[M],
                 int polyn[M][T][M], double coeff[M][T]);

  // Data Members--------------------------------------------------------------
  vcl_vector<vnl_real_npolynomial*> ps_;   // the input
  vcl_vector<vnl_vector<double>*> r_; // the output (real part)
  vcl_vector<vnl_vector<double>*> i_; // the output (imaginary part)
};

#endif // vnl_rnpoly_solve_h_