// hueckel-basis.h    Section: 6.15, Hueckel Edge Operator
//      Copyright 1994, Center for Computer Vision and Visualization,
//      University of Florida.  All rights reserved.
// Hueckel Edge Operator support functions


#ifndef _hueckel_basis_h_
#define _hueckel_basis_h_

#include <math.h>
#include "FloatDI.h"


//====================================================================
// Support functions for Hueckel operator
//====================================================================


// These include ALL constants that are simplied multiplied
//   or divided in the Hueckel basis functions.
class HueckelConst {
public:
    static double HC[ 8 ];
    static double hc[ 8 ];

    HueckelConst() {
	HC[0] = 5.0 * M_PI / 6.0;
	HC[1] = sqrt( 8.0 * M_PI / 27.0 );
	HC[2] = sqrt( 3.0 / M_PI );
	HC[3] = HC[2];
	HC[4] = HC[2];
	HC[5] = HC[2];
	HC[6] = sqrt( 32.0 / 3.0 / M_PI );
	HC[7] = HC[6] * 2.0;

	hc[0] = sqrt( 24.0 * M_PI / 25.0 );
	hc[1] = sqrt( 9.0/2.0 );
	hc[2] = 1;
	hc[3] = 1;
	hc[4] = 3.0/2.0;
	hc[5] = hc[4];
	hc[6] = sqrt( 5.0/4.0 );
	hc[7] = hc[6];
    }
};



class HueckelBasis {
public:
    static const int NUM_BASIS_FUNCTIONS;

    HueckelBasis( const double R_, const IA_Point<int> & origin );
    ~HueckelBasis();

    // The Hueckel H() basis functions
    double H0 (const IA_Point<int> &p) const;
    double H1 (const IA_Point<int> &p) const;
    double H2 (const IA_Point<int> &p) const;
    double H3 (const IA_Point<int> &p) const;
    double H4 (const IA_Point<int> &p) const;
    double H5 (const IA_Point<int> &p) const;
    double H6 (const IA_Point<int> &p) const;
    double H7 (const IA_Point<int> &p) const;

    double H( unsigned int index, const IA_Point<int> & ip ) const;


    // The Hueckel h() basis functions
    double h0 (const IA_Point<int> &p) const;
    double h1 (const IA_Point<int> &p) const;
    double h2 (const IA_Point<int> &p) const;
    double h3 (const IA_Point<int> &p) const;
    double h4 (const IA_Point<int> &p) const;
    double h5 (const IA_Point<int> &p) const;
    double h6 (const IA_Point<int> &p) const;
    double h7 (const IA_Point<int> &p) const;

    double h( unsigned int index, const IA_Point<int> & ip ) const;

    double r( const IA_Point<int> & ) const;

    double Q( double ) const;
    double Q( const IA_Point<int> & ) const;

    typedef  double (HueckelBasis::*PFD_IP) (const IA_Point<int> &) const;

    PFD_IP 
    operator() (const unsigned int) const;

private:
    double R;
    IA_Point<int> p0;

    PFD_IP
    *basis_H, //[ NUM_BASIS_FUNCTIONS ],
    *basis_h; //[ NUM_BASIS_FUNCTIONS ];

    double basis_error (const IA_Point<int> &p) const;
};



//====================================================================
// HueckelImage2D()
// Return the Hueckel image centered at the origin, and
//   parameterized by the radius and basis function
//====================================================================

IA_Image<IA_Point<int>,float>
HueckelImage2D( unsigned int basis_index, double radius );

#endif
