#ifndef POLYNOMIAL_HH
#define POLYNOMIAL_HH
////////////////////////////////////////////////////////////////////////////
//! file="amma/Geometry/AnGeoN/Polynomi.hh"
//! lib=MagN
//! userlevel=Normal
//! author="Radek Marik"
//! date="20.10.92"
//! docentry="Geometry.N-D;Optimisation.Numerical;Numerical Methods"
//! rcsid="$Id: Polynomi.hh,v 1.8 1999/09/16 13:30:52 ees1cg Exp $"

#include "amma/StdType.hh" //RealT
#include "amma/Index.hh"
#include "amma/Vector.hh"

class ostream;
class istream;
class ByteImageC;
class DoubleImageC;

//=======================================================================
// ------------ PolynomialC ---------------------------------------------
//=======================================================================

//: Algebraic polynomial

class PolynomialC : public VectorC
{
  public:

    PolynomialC(): VectorC(), vars(0), order(0) {}
    PolynomialC(SizeT var, SizeT ord);
    PolynomialC(SizeT var, SizeT ord, const VectorC & vec);
    PolynomialC(const PolynomialC & pol);

    PolynomialC(const ByteImageC & image, SizeT ord);
    // Computes the algebraic polynomial of the order 'ord' which is
    // a least square fit of the image function, the coordinates are
    // relative in the range <-1,1> x <-1,1>.

    PolynomialC(const ByteImageC & image, SizeT ord, const ByteImageC & mask);
    // Computes the algebraic polynomial of the order 'ord' which is
    // a least square fit of the image function, the coordinates are
    // relative in the range <-1,1> x <-1,1>. The polynomial is
    // computed using the pixel values which are marked in the image 'mask'
    // by the label different from 'BLACK'.

    PolynomialC(const DoubleImageC & image, SizeT ord);
    // compute the algebraic polynomial of the order 'ord' which is
    // a least square fit of the image function, the coordinates are
    // relative in the range <-1,1> x <-1,1>

    PolynomialC & Polynomial();
    const PolynomialC & Polynomial() const;
    virtual ~PolynomialC();

    SizeT Dim(SizeT var, SizeT order) const;
    SizeT Dim() const;
    SizeT VarsNumber() const;
    SizeT Order() const;
    PolynomialC Copy() const;
    PolynomialC & Set(const VectorC & vec);
    PolynomialC & Set(const PolynomialC & pol);
    const PolynomialC & operator=(const PolynomialC & pol);

    // algebraic polynomials
    //======================
    VectorC VarBase() const;
    VectorC VarBase(RealT x) const;
    VectorC VarBase(RealT x, RealT y) const;
    VectorC VarBase(RealT x, RealT y, RealT z) const;

    RealT Value() const;
    RealT Value(RealT x) const;
    RealT Value(RealT x, RealT y) const;
    RealT Value(RealT x, RealT y, RealT z) const;

    RealT DifValue(RealT f) const;
    RealT DifValue(RealT x, RealT f) const;
    RealT DifValue(RealT x, RealT y, RealT f) const;
    RealT DifValue(RealT x, RealT y, RealT z, RealT f) const;

    PolynomialC Dif(IndexT count) const;
    PolynomialC DifX(IndexT count) const;
    PolynomialC DifY(IndexT count) const;

    PolynomialC Dif() const   {return Dif(1);}
    PolynomialC DifX() const  {return DifX(1);}
    PolynomialC DifY() const  {return DifY(1);}
    PolynomialC DifXY() const {return DifX(1).DifY(1);}

    // Chebyshev polynomials
    //======================
    VectorC VarBaseC() const;
    VectorC VarBaseC(RealT x) const;
    VectorC VarBaseC(RealT x, RealT y) const;
    VectorC VarBaseC(RealT x, RealT y, RealT z) const;

    RealT ValueC() const;
    RealT ValueC(RealT x) const;
    RealT ValueC(RealT x, RealT y) const;
    RealT ValueC(RealT x, RealT y, RealT z) const;

  private:
    SizeT vars;  // number of variables
    SizeT order; // polynomial order
 
  friend istream & operator>>(istream & s, PolynomialC & pol);
};

ostream & operator<<(ostream & s, const PolynomialC & pol);
istream & operator>>(istream & s, PolynomialC & pol);

#endif

// IAPS - Image analysis program system. 
// End of include file Polynomial.hh