head	1.4;
access;
symbols;
locks; strict;
comment	@// @;


1.4
date	95.01.06.05.16.52;	author rjj;	state Exp;
branches;
next	1.3;

1.3
date	94.12.02.06.22.04;	author rjj;	state Exp;
branches;
next	1.2;

1.2
date	94.08.29.00.55.56;	author rjj;	state Exp;
branches;
next	1.1;

1.1
date	94.06.20.07.10.39;	author rjj;	state Exp;
branches;
next	;


desc
@Initial version for the Hueckel algorithm suite.
@


1.4
log
@Added support for initalization of NUM_BASIS_FUNCTIONS.
Figured out pointers to member functions and moved code to use this
  approach rather than the switch stmts use perviously.
Added implementation of destructor for HueckelBasis.
@
text
@// 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 );
}
@


1.3
log
@Deleted unreachable stmts that were generating warning msgs.
@
text
@d1 1
a1 1
// hueckel-basis.c    Section: 6.15  Hueckel Edge Operator
d16 1
a16 1
// static const int HueckelBasis::NUM_BASIS_FUNCTIONS = 8;
a20 4
//
// NOTE: The H() basis functions have been modified to return 0
//       when the  r()  value is greater than 1.
//       This simplifies the defintion of the HueckelTemplate
d28 5
a32 19
    switch (index) {
    case 0:
	return this->H0(ip);
    case 1:
	return this->H1(ip);
    case 2:
	return this->H2(ip);
    case 3:
	return this->H3(ip);
    case 4:
	return this->H4(ip);
    case 5:
	return this->H5(ip);
    case 6:
	return this->H6(ip);
    case 7:
	return this->H7(ip);
    default:
	return this->basis_error(ip);
d34 2
d95 4
a98 19
    switch (index) {
    case 0:
	return this->h0(ip);
    case 1:
	return this->h1(ip);
    case 2:
	return this->h2(ip);
    case 3:
	return this->h3(ip);
    case 4:
	return this->h4(ip);
    case 5:
	return this->h5(ip);
    case 6:
	return this->h6(ip);
    case 7:
	return this->h7(ip);
    default:
	return this->basis_error(ip);
d100 2
d147 2
a148 1
    return p(0)-1;
d158 1
a158 1
				  const IA_Point<int> & origin )
d164 7
a170 7
    this->basis_H[0] = &HueckelBasis::H1;
    this->basis_H[0] = &HueckelBasis::H2;
    this->basis_H[0] = &HueckelBasis::H3;
    this->basis_H[0] = &HueckelBasis::H4;
    this->basis_H[0] = &HueckelBasis::H5;
    this->basis_H[0] = &HueckelBasis::H6;
    this->basis_H[0] = &HueckelBasis::H7;
d175 7
a181 7
    this->basis_h[0] = &HueckelBasis::h1;
    this->basis_h[0] = &HueckelBasis::h2;
    this->basis_h[0] = &HueckelBasis::h3;
    this->basis_h[0] = &HueckelBasis::h4;
    this->basis_h[0] = &HueckelBasis::h5;
    this->basis_h[0] = &HueckelBasis::h6;
    this->basis_h[0] = &HueckelBasis::h7;
d184 4
d191 1
a191 2
//    if (index < HueckelBasis::NUM_BASIS_FUNCTIONS) {
    if (index < 8) {
d196 1
a196 1
//    IA::ia_throw(IA::INVALID_OPERATION, __FILE__, __LINE__);
@


1.2
log
@Deleted include of <math.h> stmt.
Modified enorm() call to use the non-member version.
@
text
@d1 1
a1 1
// hueckel-basis.c    Section: 6.15
a51 1
    return this->basis_error(ip);
a130 1
    return this->basis_error(ip);
@


1.1
log
@Initial revision
@
text
@a5 1
#include <math.h>
a10 1

d231 1
a231 1
	(this->p0 - ip).enorm() / R;
@
