/************************************************************************
 *                                                                      *
 *  Program package         T O O L D I A G                             *
 *                                                                      *
 *  Version   1.5                                                       *
 *  Date:     Tue Feb  8 13:39:06 1994                                  *
 *                                                                      *
 *  NOTE: This program package is copyrighted in the sense that it      *
 *  may be used for scientific purposes. The package as a whole, or     *
 *  parts thereof, cannot be included or used in any commercial         *
 *  application without written permission granted by the author.       *
 *  No programs contained in this package may be copied for commercial  *
 *  distribution.                                                       *
 *                                                                      *
 *  All comments  concerning this program package may be sent to the    *
 *  e-mail address 'tr@fct.unl.pt'.                                     *
 *                                                                      *
 ************************************************************************/



#include <stdio.h>
#include <math.h>
#include "def.h"
#include "featslct.h"
#include "matrix.h"

extern universe *U;
extern MatElem trace();
extern int euclid_dist;

static str80 buf;

#define PEEPMAT(M) printf("M=\n");print_Matrix(&M);DBG;   /**/
/* #define PEEPMAT(M)   /**/
#define SOON {printf("Soon at the theater near to you...type <Ctrl> c");d=-1.0;DBG;}


float Euclidian_Distance( S1, S2, dim )
FeatVector S1, S2;
int dim;
{
 int i;
 float aux, D = 0.0;

 for( i = 0; i < dim; i++ )
 {
   aux = S1[i] - S2[i];
   D += aux * aux;
 }
 return( (float)sqrt((double)D) );
}


float vector_len( V, dim )
FeatVector V;
int dim;
{
 int i;
 float L = 0.0;

 for( i = 0; i < dim; i++ )
   L += V[i] * V[i];
 return( (float)sqrt((double)L) );
}


void norm_vector( V, dim )
FeatVector V;
int dim;
{
 int i; float len;

 len = vector_len( V, dim );
 if( len == 0.0 )
 {
   printf("Error in norm_vector: vector has 0 length - exit ...\n");
   exit(1);
 }
 for( i = 0; i < dim; i++ )
   V[i] /= len;
}


float dot_product( S1, S2, dim )
FeatVector S1, S2;
int dim;
{
 int i;
 float sum = 0.0;

 for( i = 0; i < dim; i++ )
   sum += S1[i] * S2[i];
 return( (float)sum );
}


float selectINTER_CLASS_DISTANCE_Euclid( FSV, len )
FeatSelectVector FSV;
int len;
{
 Matrix T, SB, SW, mat1, invSW, invT;
 MatElem dd, traceSW, DetSW, DetT;
 float d;
 init_Matrix( &T ); init_Matrix( &SB ); init_Matrix( &SW );
 init_Matrix( &mat1 ); init_Matrix( &invSW ); init_Matrix( &invT );
 scatter_mats_kittler( FSV, len, &T, &SB, &SW );

 switch( euclid_dist )
 {
   case J1 : dd = trace( &T ); break;
   case J2 : traceSW = trace( &SW );
     if( traceSW != 0.0 )
       dd = trace( &SB ) / traceSW;
     else
       {fprintf(stderr,"Had a zero SW scatter matrix - exit..."); exit(1);}
     break;
   case J3 : Copy_Matrix( &SW, &invSW );
     Invert_Matrix( &invSW, &DetSW );
     Mult_Matrix( &invSW, &SB, &mat1 );
     dd = trace( &mat1 );
     break;
   case J4 : Copy_Matrix( &SW, &invSW );
     Invert_Matrix( &invSW, &DetSW );
     Copy_Matrix( &T, &invT );
     Invert_Matrix( &invT, &DetT );
     dd = DetT / DetSW;
     break;
   default : fprintf(stderr,"What is Euclidean distance %d? - exit...\n",
	euclid_dist ); exit(1);
 }
 d = (float)dd;
 /* printf("d=%f", d ); DBG;     /**/
 Matrix_Free( &T ); Matrix_Free( &SB ); Matrix_Free( &SW );
 Matrix_Free( &mat1 ); Matrix_Free( &invSW ); Matrix_Free( &invT );
 return( d );
}