// hueckel-basis.c    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

#include "PSIter.h"
#include "hueckel-basis.h"
#include "shaped-pointset.h"


//====================================================================
// class HueckelBasis
// Support functions for Hueckel operator
//====================================================================

const int HueckelBasis::NUM_BASIS_FUNCTIONS = 8;


//====================================================================
// Hueckel H() basis functions   
//====================================================================

// Initialize the Hueckel constants
static  HueckelConst  C;


double HueckelBasis::H(unsigned int index, const IA_Point<int> &ip) const {

    if ( index < HueckelBasis::NUM_BASIS_FUNCTIONS ) {
	PFD_IP hbH = this->basis_H[index];
	return (this->*hbH)(ip);
	//    return (*(this->basis_H[index]))(ip);  This does NOT work
    }
    // otherwise
    return this->basis_error(ip);
}


double HueckelBasis::H0 (const IA_Point<int> &p) const {
    double r = this->r(p);
    return
	::C.HC[0] * Q(r)*(1+2*r*r);
}

double HueckelBasis::H1 (const IA_Point<int> &p) const {
    double r = this->r(p);
    return
	::C.HC[1] * Q(r)*(5*r*r - 2);
}

double HueckelBasis::H2 (const IA_Point<int> &p) const { 
    double r = this->r(p);
    return
	::C.HC[2] * Q(r)*(p(0)-p0(0))*(4*r*r-1)/R;
}

double HueckelBasis::H3 (const IA_Point<int> &p) const {
    double r = this->r(p);
    return
	::C.HC[3] * Q(r)*(p(1)-p0(1))*(4*r*r-1)/R;
}

double HueckelBasis::H4 (const IA_Point<int> &p) const {
    double r = this->r(p);
    return
	::C.HC[4] * Q(r)*(p(0)-p0(0))*(3-4*r*r)/R;
}

double HueckelBasis::H5 (const IA_Point<int> &p) const {
    double r = this->r(p);
    return
	::C.HC[5] * Q(r)*(p(1)-p0(1))*(3-4*r*r)/R;
}

double HueckelBasis::H6 (const IA_Point<int> &p) const {
    double r = this->r(p);
    double x = p(0)-p0(0),  y = p(1)-p0(1);
    return
	::C.HC[6] * Q(r)*( (x*x - y*y)/R/R );
}

double HueckelBasis::H7 (const IA_Point<int> &p) const {
    double r = this->r(p);
    double x = p(0)-p0(0),  y = p(1)-p0(1);
    return
	::C.HC[7] * Q(r)*x*y/R/R;
}



//====================================================================
// Hueckel h() basis functions   
//====================================================================
double HueckelBasis::h(unsigned int index, const IA_Point<int> &ip) const {

    if ( index < HueckelBasis::NUM_BASIS_FUNCTIONS ) {
	PFD_IP hbh = this->basis_h[index];
	return (this->*hbh)(ip);
    }
    // otherwise
    return this->basis_error(ip);
}


double HueckelBasis::h0 (const IA_Point<int> &p) const {
    return
	::C.hc[0] * H(0,p);
}

double HueckelBasis::h1 (const IA_Point<int> &p) const {
    return
	::C.hc[1] * H(1,p);
}

double HueckelBasis::h2 (const IA_Point<int> &p) const { 
    return
	H(2,p) + H(4,p);
}

double HueckelBasis::h3 (const IA_Point<int> &p) const {
    return
	H(3,p) + H(5,p);
}

double HueckelBasis::h4 (const IA_Point<int> &p) const {
    return
	::C.hc[4] * H(6,p);
}

double HueckelBasis::h5 (const IA_Point<int> &p) const {
    return
	::C.hc[5] * H(7,p);
}

double HueckelBasis::h6 (const IA_Point<int> &p) const {
    return
	::C.hc[6] * ( H(2,p) - H(4,p) );
}

double HueckelBasis::h7 (const IA_Point<int> &p) const {
    return
	::C.hc[7] * ( H(3,p) - H(5,p) );
}


double HueckelBasis::basis_error (const IA_Point<int> &p) const {
    ia_throw( IA::INVALID_OPERATION, __FILE__, __LINE__ );
    return -1;
};
    


//====================================================================
// class HueckelBasis
//====================================================================

HueckelBasis::HueckelBasis( const double R_,
			    const IA_Point<int> & origin )
: R( R_ ), p0( origin ) {

    basis_H = new HueckelBasis::PFD_IP[8];

    this->basis_H[0] = &HueckelBasis::H0;
    this->basis_H[1] = &HueckelBasis::H1;
    this->basis_H[2] = &HueckelBasis::H2;
    this->basis_H[3] = &HueckelBasis::H3;
    this->basis_H[4] = &HueckelBasis::H4;
    this->basis_H[5] = &HueckelBasis::H5;
    this->basis_H[6] = &HueckelBasis::H6;
    this->basis_H[7] = &HueckelBasis::H7;

    basis_h = new HueckelBasis::PFD_IP[8];

    this->basis_h[0] = &HueckelBasis::h0;
    this->basis_h[1] = &HueckelBasis::h1;
    this->basis_h[2] = &HueckelBasis::h2;
    this->basis_h[3] = &HueckelBasis::h3;
    this->basis_h[4] = &HueckelBasis::h4;
    this->basis_h[5] = &HueckelBasis::h5;
    this->basis_h[6] = &HueckelBasis::h6;
    this->basis_h[7] = &HueckelBasis::h7;
}

HueckelBasis::~HueckelBasis() {
    delete[] basis_H;
    delete[] basis_h;
}

HueckelBasis::PFD_IP 
HueckelBasis::operator() ( const unsigned int index ) const {
    if (index < HueckelBasis::NUM_BASIS_FUNCTIONS) {
	// For computational efficiency we will use the H()
	//    Hueckel functions as the vector space basis
	return this->basis_H[index];
    }
    ia_throw(IA::INVALID_OPERATION, __FILE__, __LINE__);
    return &HueckelBasis::basis_error;
}


double
HueckelBasis::r( const IA_Point<int> &ip ) const {
    return
	enorm(this->p0 - ip) / R;
}


double
HueckelBasis::Q( double r_val ) const {
    return
	sqrt( 1 - r_val*r_val );
}


double
HueckelBasis::Q( const IA_Point<int> &ip ) const {
    double r_val( this->r(ip));
    return
	sqrt( 1 - r_val*r_val );
}



//====================================================================
// 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 ) {

    IA_Point<int>  origin(0,0);

    // The Hueckel basis
    HueckelBasis
	hb( radius, origin );

    // The Hueckel image domain
    IA_Set<IA_Point<int> >
	hi_ps( digital_disc( origin, radius ) );

    // The Hueckel image range
    float *hi_vals = new float[ hi_ps.card() ];

    // Let us fill in the values...
    IA_Point<int>  p;
    IA_PSIter<IA_Point<int> >  hipIter( hi_ps );

    int icount = 0;
    while( hipIter( p ) ) {
	hi_vals[ icount++ ] = hb.H(basis_index,p);
    }

    int ASSUME_STORAGE_MGT = 1;
    return
	IA_Image<IA_Point<int>,float>( hi_ps,
				       hi_vals,
				       hi_ps.card(),
				       ASSUME_STORAGE_MGT );
}
