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

#include "hueckel.h"
#include "DDTempl.h"
#include "FloatConv.h"
#include "hueckel-basis.h"



//==========================================================
// Hueckel algorithm
//==========================================================


IA_Image<IA_Point<int>,float>
hueckel2D( const IA_Image<IA_Point<int>,float> aImage,
           double radius,
           double theta[],        // Array of angles, in radians
           unsigned theta_m )     // count of number of angles
{
    // Declare working images we will need
    IA_Image<IA_Point<int>,float>
	a[8],  e0, e1, e2,  u,  *lambda;

    lambda = new IA_Image<IA_Point<int>,float>[theta_m];


    unsigned int i, j;

    // Since were are only doing the simple part of the Hueckel
    //   algorithm here, then we only need to calculate
    //   a[1..7], we do NOT need to calculate a[0]

    for( i=1; i<8; i++) {
	a[i] =
	    linear_product(
		aImage,
		IA_DDTemplate< IA_Image<IA_Point<int>,float> >(
		    2, HueckelImage2D( i, radius ) )
		);
    }


    double cosj, sinj;
    for( j=0; j<theta_m; j++) {
	cosj = cos(theta[j]);
	sinj = sin(theta[j]);

	e0 = a[2]*cosj + a[3]*sinj;
	e1 = a[4]*cosj + a[5]*sinj;
	e2 = a[1]
	    + a[6]*( cosj*cosj - sinj*sinj ) + 2*a[7]*sinj*cosj;

	// Use pseudo_div() in case divisor has zero values.
	u = pseudo_div( e0, fabs( e0 ) )
	    *sqrt( sqr( e1 ) + sqr( e2 ) );

	lambda[j] = e0 + u;
    }


    // Now for the final stuff
    IA_Image<IA_Point<int>,float>
	Lambda( fabs( lambda[0] ) ),
	k, e;

    for( j=1; j<theta_m; j++) {
	Lambda = max( Lambda, fabs( lambda[j] ) );
    }

    k = pseudo_div( Lambda,
		    sqrt( 6.0*sqr(a[1]) +
			  2.0*(sqr(a[2])+sqr(a[3])+sqr(a[4])+sqr(a[5])) +
			  3.0*(sqr(a[6])+sqr(a[7]))
			  )
		    );

    delete[] lambda;  // delete this memory after we are done with it.

    return k;
}



IA_Image<IA_Point<int>,u_char>
hueckel2D( const IA_Image<IA_Point<int>,float> aImage,
           double radius,
           double theta[],        // Array of angles, in radians
           unsigned theta_m,      // count of number of angles
           float threshold ) {

    IA_Image<IA_Point<int>,float>
	e = chi_gt( hueckel2D( aImage, radius, theta, theta_m ),
		    threshold );

    return to_uchar( e );
}
