// mat.c                Section: 4.4  Medial Axis Transform (MAT)
//      Copyright 1994, Center for Computer Vision and Visualization,
//      University of Florida.  All rights reserved.
// Medial Axis Transform

#include "Nbh.h"
#include "mat.h"


IA_Image<IA_Point<int>,int>
medial_axis( const IA_Image<IA_Point<int>,IA_Bit> & aImage,
	     const IA_Set<IA_Point<int> > & ball,
	     unsigned int & count ) {

    IA_Neighborhood<IA_Point<int>,IA_Point<int> >
	m_nbhd( aImage.domain().dim(), ball );

    IA_Image<IA_Point<int>,int>  m( aImage.domain(), 0 );

    IA_Image<IA_Point<int>,IA_Bit>
	a = aImage,
	b( aImage.domain(), IA_Bit(0) );

    count = 0;
    while ( a != IA_Bit( 0 ) ) {
	a = min( a, m_nbhd );

	// Remember a is binary valued image, ~x == chi(0,x)
	b = a & ~( max( min( a, m_nbhd ),
			m_nbhd ) );
	++count;
	m += count * IA_Image<IA_Point<int>,int>( b );
    }

    return m;
}


IA_Image<IA_Point<int>,IA_Bit>
inverse_medial_axis( const IA_Image<IA_Point<int>,int> & m,
		     const IA_Set<IA_Point<int> > & ball )  {

    // Each non-zero value in m corresponds to where a ball of
    //   the appropriate radius value should be placed.
    // However, the ball\'s radius is actually one less than
    //   the image value, since a value of one (1)
    //   corresponds to a single point.
    // Thus, radus value > 0 => k = radius value -1

    IA_Set<IA_Point<int> >
	ball_ps   = extend_to_point( 0, m.domain().dim() ),
	result_ps = IA_empty_ipset( m.domain().dim() );

    int  maxK = max( m );

    for( int k=1; k<=maxK; k++ ) {
	result_ps |= m.restrict( IA_Set<int>( k ) ).domain()
	    + ball_ps;

	// Grow the ball for the next loop through
	ball_ps += ball;  // Minkowski sum of point sets
    }

    // The result is an image over the same domain but with the
    //    corresponding inverse medial axis masked in with 1\'s.
    IA_Image<IA_Point<int>,IA_Bit>
	result_image( m.domain(), IA_Bit(0) ),
	ball_image( result_ps, IA_Bit(1) );

    return
	result_image.restrict_assign( ball_image );
}

