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


1.8
date	95.01.31.07.20.28;	author rjj;	state Exp;
branches;
next	1.7;

1.7
date	95.01.08.22.11.54;	author rjj;	state Exp;
branches;
next	1.6;

1.6
date	94.08.29.02.45.26;	author rjj;	state Exp;
branches;
next	1.5;

1.5
date	94.07.23.02.51.05;	author rjj;	state Exp;
branches;
next	1.4;

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

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

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

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


desc
@First cut at the Hough line detection algorithm.
@


1.8
log
@Some simple formatting and comment changes.
Also changed to use rounding rather than trancation.
@
text
@// hough.c              Section: 10.3, Line Detection using the Hough Transform
//      Copyright 1994, Center for Computer Vision and Visualization,
//      University of Florida.  All rights reserved.

#include "PSIter.h"
#include "DDTempl.h"
#include "UcharConv.h"
#include "hough.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-1);
	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   = nint( (rows-1)/R/2.0*( rho + R ) ); // round to nearest
	ipa[c][0]   = r;  // 1st co-ord of c-th element of the point arrary
    }

    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 );
}
@


1.7
log
@Added standard header.
Reformatted some lines to shrinken the length of the file,
   or for readibility with surrounding lines.
@
text
@d53 2
a54 2
	ipa[c]       = IA_Point<int>( 0, c );
	double theta = M_PI * double( c ) / cols;
d73 1
a73 1
	ipa[ c ]     = ht.ipa[c];
d108 2
a109 2
	int     r   = int( rows/R/2.0*( rho + R ) );
	ipa[c][0]   = r;  // 1st co-ordinate of work array of points
@


1.6
log
@Modified include stmts to be iac++ 2.0 compliant
@
text
@d1 1
a1 1
// hough.c              Section: 10.3
a4 2

#include "hough.h"
d8 1
a8 1

d32 1
a32 2
  rows( 2*positive_rows_requested + 1 ),
  cols( cols_requested ),
d53 1
a53 1
	ipa[c] = IA_Point<int>( 0, c );
a61 1

d80 1
a80 2
HoughTemplate::
~HoughTemplate() {
d91 1
a91 2
	return
	    IA_Image<IA_Point<int>,u_char>();  // Empty image
d108 2
a109 3

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


1.5
log
@Checking in this version before adding additional functionality that
  will calculate the 'reverse' hough, to see what lines were found.
@
text
@d9 1
a9 1
#include "UcharConvolutions.h"
@


1.4
log
@Added speedup enhancement of co/sines[].
@
text
@d15 2
a16 3
       unsigned int columns,
       double &R ) {

a18 2
    R = ht.R;

a27 15
void init_co_sines( double *cosines, double *sines, unsigned int N ) {
    if ( (cosines == 0) || (sines == 0) ) {
	ia_throw( IA::INVALID_OPERATION, __FILE__, __LINE__ );
	// This is really bad, because we can not recover from this;
	exit( 1 );
    }
    // OK, assuming we are here, fill in the values
    for ( int i=0; i<N; i++ ) {
	double theta = M_PI * double( i ) / double( N );
	cosines[i] = cos( theta );
	sines[i]   = sin( theta );
    }
}


d33 1
a33 1
  b( b_ ),
d39 2
a40 1
  b_inf( b_.domain().inf() )
d48 2
a49 2
    // This array servers as a convient point array for the filling
    //    in of the Hough line  x*cos(theta) + y*sin(theta)
d51 4
a54 3
    //    expresion defined above,
    //    but the second co-ordinate will represent theta, so
    //    it can be filled in now.
d57 3
a60 3
    cosines = new double[ cols_requested ];
    sines   = new double[ cols_requested ];
    init_co_sines( cosines, sines, cols_requested );
d69 5
a73 2
  b( ht.b ), rows( ht.rows ), cols( ht.cols ), R( ht.R ),
  ipa( new IA_Point<int>[ int(ht.cols) ] ), b_inf( ht.b_inf )
d77 3
a79 1
	ipa[c] = ht.ipa[c];
a80 3
    cosines = new double[ int(ht.cols) ];
    sines   = new double[ int(ht.cols) ];
    init_co_sines( cosines, sines, (unsigned int)( cols ) );
@


1.3
log
@Corrected column limit in operator()
@
text
@a11 1

d31 14
d74 3
d86 9
a94 6
  ipa( new IA_Point<int>[ int(ht.cols) ] ), b_inf( ht.b_inf ) {

      // Copy the points over...
      for( int c=0; c<cols; c++ ) {
	  ipa[c] = ht.ipa[c];
      }
d101 2
d127 1
a127 2
	double  theta = M_PI * double( c ) / cols;
	double  rho   = x*cos( theta ) + y*sin( theta );
a139 1

@


1.2
log
@This version seems to work, but is not real fast.
  There are still some debug stmts left, but I'm putting
  the file under rcs before working with it further.
@
text
@d106 1
a106 1
	double  theta = M_PI * double( c ) / (cols+1);
@


1.1
log
@Initial revision
@
text
@d11 18
a28 11
static
int R_max( int i0, int i1, int s0, int s1 ) {
    int m[4];
    m[0] = i0 < 0 ? -i0 : i0;    m[1] = i1 < 0 ? -i1 : i1;
    m[2] = s0 < 0 ? -s0 : s0;    m[3] = s1 < 0 ? -s1 : s1;

    int r_max = m[0];
    for( int i=1; i<4; i++ ) {
	if (m[i] > r_max) r_max = m[i];
    }
    return r_max;
d32 1
d35 2
a36 2
	       unsigned int rows_, unsigned int cols_,
	       double epsilon_ )
d38 7
a44 1
  b( b_ ), rows( rows_ ), cols( cols_ ), epsilon( epsilon_ )
d48 1
a48 1
	R = 0;
d52 9
a60 2
    IA_Point<int> i( b_.domain().inf() ), s( b_.domain().sup() );
    R = R_max( i(0), i(1), s(0), s(1) );
d64 1
d66 11
a76 2
IA_Image<IA_Point<int>,u_char>
HoughTemplate:: operator() ( const IA_Point<int> & xy ) const {
a77 3
    IA_Set<IA_Point<int> >
	ips( IA_boxy_pset( IA_Point<int>( 0, 0 ),
			   IA_Point<int>( rows-1, cols-1 ) ) );
d79 4
a82 2
    IA_Image<IA_Point<int>,u_char>
	rTheta( ips, u_char( 0 ) );
a83 3
    if ( b( xy ) == 0 ) {
	return  rTheta;
    }
d85 2
d88 3
a90 18
    // We know the domain of the returned image, so we can
    //    iterate over that domain and fill in the 1s.
    // Calculate the image values at the rho, theta points
    IA_PSIter<IA_Point<int> >
	psIter( ips );

    IA_Point<int> rt;
    double  x = xy( 0 ), y = xy( 1 );
    while( psIter( rt ) ){
	double
	    i = rt( 0 ), j = rt( 1 ),
	    rho   = 2*R*( i - rows/2 ) / rows,
	    theta = M_PI * j / cols;

	if ( abs( rho - ( x*cos( theta ) + y*sin( theta ) ) )
	     < epsilon ) {
	    rTheta[ rt ] = 1;
	}
d92 13
a104 2
    return rTheta;
}
d106 2
d109 3
d113 7
a120 3
IA_Image<IA_Point<int>,int>
Hough( const IA_Image<IA_Point<int>,IA_Bit> b,
       unsigned int rows, unsigned int columns, double epsilon ) {
a121 3
    IA_Set<IA_Point<int> >
	dest_ps( IA_boxy_pset( IA_Point<int>( 0, 0 ),
			       IA_Point<int>( rows-1, columns-1 ) ) );
d123 3
a125 5
    return
	sum( b.domain(),
	     IA_DDTemplate<IA_Image<IA_Point<int>,u_char> >(
		 HoughTemplate( b, rows, columns, epsilon ) ),
	     dest_ps );
@
