// hough.c              Section: 10.3
//      Copyright 1994, Center for Computer Vision and Visualization,
//      University of Florida.  All rights reserved.


#include "hough.h"
#include "PSIter.h"
#include "DDTempl.h"
#include "UcharConv.h"


IA_Image<IA_Point<int>,int>
Hough( const IA_Image<IA_Point<int>,IA_Bit> b,
       unsigned int positive_rows_requested,
       unsigned int columns )
{
    HoughTemplate ht( b, positive_rows_requested, columns );

    int rows = 2*positive_rows_requested + 1;
    return
	sum( b.domain(),
	     IA_DDTemplate<IA_Image<IA_Point<int>,u_char> >( ht ),
	     IA_boxy_pset( IA_Point<int>(      0,         0 ),
			   IA_Point<int>( rows-1, columns-1 ) ) );
}


HoughTemplate::
HoughTemplate( const IA_Image<IA_Point<int>,IA_Bit> & b_,
	       unsigned int positive_rows_requested,
	       unsigned int cols_requested )
: IA_ClosureDT<IA_Image<IA_Point<int>,u_char> >( b_.domain() ),
  b( b_ ), b_inf( b_.domain().inf() ),
  rows( 2*positive_rows_requested + 1 ),
  cols( cols_requested ),
  R( (b_.domain().sup() - b_.domain().inf()).enorm() ),

  ipa( new IA_Point<int>[ cols_requested ] ),
  cosines( new double[ cols_requested ] ),
  sines(   new double[ cols_requested ] )
{
    if ( b_.domain().dim() != 2 ) {
	ia_throw( IA::INVALID_OPERATION, __FILE__, __LINE__ );
	*((int *) &R) = 0;  // its const by here
	return;
    }

    // The  ipa  array allows us to only allocate the points once
    //    defined by the Hough line  x*cos(theta) + y*sin(theta)
    // The first co-ordinate value is a function of sinusoidal
    //    expresion defined above, the second encodes theta,
    //    so it can be filled in now.
    // cosines[] and sines[] store the precomputed cos(theta) and
    //    sin(theta) values for each encoded theta.
    for( int c=0; c<cols; c++ ) {
	ipa[c] = IA_Point<int>( 0, c );
	double theta = M_PI * double( c ) / cols;
	cosines[ c ] = cos( theta );
	sines[   c ] = sin( theta );
    }
}


// Copy constructor, required because we allocate memory

HoughTemplate::
HoughTemplate( const HoughTemplate & ht )
: IA_ClosureDT<IA_Image<IA_Point<int>,u_char> >( ht.domain() ),
  b( ht.b ), b_inf( ht.b_inf ),
  rows( ht.rows ), cols( ht.cols ), R( ht.R ),
  ipa( new IA_Point<int>[ int(ht.cols) ] ),
  cosines( new double[ int(cols) ] ),
  sines(   new double[ int(cols) ] )
{
    // Initialize our internals
    for( int c=0; c<cols; c++ ) {
	ipa[ c ]     = ht.ipa[c];
	cosines[ c ] = ht.cosines[c];
	sines[   c ] = ht.sines[c];
    }
}


HoughTemplate::
~HoughTemplate() {
    delete[] ipa;
    delete[] cosines;
    delete[] sines;
}


IA_Image<IA_Point<int>,u_char>
HoughTemplate:: operator() ( const IA_Point<int> & xy ) const {

    if ( b( xy ) == IA_Bit( 0 ) ) {
	return
	    IA_Image<IA_Point<int>,u_char>();  // Empty image
    }

    // The image we will return has 1s everywhere it is defined,
    //    but it will only be defined over the sinusoidal curve
    //    the point xy defines in the accumulator array.

    double
	x = xy( 0 ) - b_inf(0),  // Compensate for non <0 0>
	y = xy( 1 ) - b_inf(1);  // based images

    for( int c=0; c<cols; c++ ) {
	// Map c to theta, calculate the corresponding r
	//   and assign that as the first co-ordiante value
	//   for this point.

	double  rho = x*cosines[c] + y*sines[c];

	int r = int( rows/R/2.0*( rho + R ) );
	ipa[c][0] = r;  // 1st co-ordinate of work array of points
    }

    return
	IA_Image<IA_Point<int>,u_char>(
	    IA_Set<IA_Point<int> >( (unsigned int) 2,        // 2D
				    ipa,                     // vector
				    (unsigned int) (cols) ), // constructor
	    u_char( 1 ) );
}


IA_ClosureDT<IA_Image<IA_Point<int>,u_char> > *
HoughTemplate::clone_self() const {
    return new HoughTemplate( *this );
}
