// sis.c                Section: 3.8
//      Copyright 1994, Center for Computer Vision and Visualization,
//      University of Florida.  All rights reserved.
// Image algebra C++ program to determine the threshold
//      level for an image based on SIS
//

#include "sis.h"
#include "DDTempl.h"
//#include "UcharConv.h"  // This needs to be specialized
#include "IntConv.h"
#include "FloatConv.h"  // This requires abs(FDI), provided herein

static
IA_Image<IA_Point<int>,float>
abs( const IA_Image<IA_Point<int>,float> & fdi ) {

    return fabs( fdi );
}



template <class P, class T>
T
sis_threshold( const IA_Image<P,T> & a ) {

    // Construct the s & t images used for templates.

    IA_Image<P,T>
	s_img( IA_Set<IA_Point<int> >( IA_Point<int>(0,-1) ) |
	       IA_Set<IA_Point<int> >( IA_Point<int>(0, 1) ),
	       T( 1 ) ),
	t_img( IA_Set<IA_Point<int> >( IA_Point<int>(-1,0) ) |
	       IA_Set<IA_Point<int> >( IA_Point<int>( 1,0) ),
	       T( 1 ) );

    s_img[ IA_Point<int>( 0,1) ] = T( -1 );  // Assign correct value.
    t_img[ IA_Point<int>(-1,0) ] = T( -1 );  // Assign correct value.

    // Construct the s & t templates
    unsigned int a_dimen = a.domain().dim();
    IA_DDTemplate<IA_Image<P,T> >
	s( a_dimen, s_img ),
	t( a_dimen, t_img );

    // Calculate the images we want
    T tau;
    IA_Image<P,T>
	e1 = abs( linear_product(a,s) ),
	e2 = abs( linear_product(a,t) ),
	e = max(e1,e2);

    tau = sum( e * abs(a) )
	  / sum( e );

    return tau;
}
