1 /*--------------------------------------------------------------------*/ 2 /*----- 02 Mar 2001 - RWCox - Used for AUTOGZIP decision -------------*/ 3 4 #include "mrilib.h" 5 6 #undef SNUM 7 #define SNUM 65536 /* one for every possible 16 bit pattern */ 8 9 static int64_t *scount=NULL ; /* holds count of unsigned shorts */ 10 static int64_t snum=0 ; /* total number of unsigned shorts processed */ 11 12 static int do_zskip = 0 ; /* skip zero? */ 13 static int do_perbin = 0 ; 14 static int do_permax = 0 ; 15 16 /*-----------------------------------------------------------------------*/ 17 ENTROPY_setup(void)18void ENTROPY_setup(void) 19 { 20 if( scount == NULL ) scount = (int64_t *) malloc(sizeof(int64_t)*SNUM) ; 21 memset( scount , 0 , sizeof(int64_t)*SNUM ) ; 22 snum = 0 ; 23 } 24 25 /*-----------------------------------------------------------------------*/ 26 ENTROPY_setdown(void)27void ENTROPY_setdown(void) 28 { 29 if( scount != NULL ){ free(scount); scount = NULL; snum = 0; } 30 } 31 32 /*-----------------------------------------------------------------------*/ 33 ENTROPY_accumulate(int64_t nbytes,void * var)34void ENTROPY_accumulate( int64_t nbytes , void * var ) 35 { 36 int64_t nn = nbytes/2 , ii ; 37 unsigned short * sar = (unsigned short *) var ; 38 39 if( scount == NULL ) ENTROPY_setup() ; 40 41 if( do_zskip ) 42 for( ii=0 ; ii < nn ; ii++ ){ 43 if( sar[ii] != 0 ){ scount[sar[ii]]++ ; snum++ ; } 44 } else { 45 for( ii=0 ; ii < nn ; ii++ ) scount[sar[ii]]++ ; 46 snum += nn ; 47 } 48 return ; 49 } 50 51 /*----------------------------------------------------------------------- 52 Value returned is in bits per 16 bit unsigned short 53 -------------------------------------------------------------------------*/ 54 ENTROPY_compute(void)55double ENTROPY_compute(void) 56 { 57 int64_t ii ; 58 double sum , bfac=1.0 ; 59 60 if( scount == NULL || snum == 0 ) return 0.0 ; 61 62 sum = 0.0 ; 63 for( ii=0 ; ii < SNUM ; ii++ ) 64 if( scount[ii] > 0 ) sum += scount[ii] * log((double)scount[ii]) ; 65 66 if( do_perbin ){ 67 int64_t nbin = 0 ; 68 for( ii=0 ; ii < SNUM ; ii++ ){ if( scount[ii] > 0 ) nbin++ ; } 69 if( nbin > 1 ) bfac = 1.0 / log((double)nbin) ; 70 } else if( do_permax ){ 71 int64_t nmax = 0 ; 72 for( ii=0 ; ii < SNUM ; ii++ ){ if( scount[ii] > 0 ) nmax = ii ; } 73 if( nmax > 1 ) bfac = 1.0 / log((double)nmax) ; 74 } 75 76 sum = -bfac * (sum - snum*log((double)snum)) / ( log(2.0) * snum ) ; 77 return sum ; 78 } 79 80 /*-----------------------------------------------------------------------*/ 81 ENTROPY_dataset(THD_3dim_dataset * dset)82double ENTROPY_dataset( THD_3dim_dataset *dset ) 83 { 84 if( !ISVALID_DSET(dset) ) return(0.0) ; 85 DSET_load(dset) ; 86 if( !DSET_LOADED(dset) ) return(0.0) ; 87 return ENTROPY_datablock( dset->dblk ) ; 88 } 89 90 /*-----------------------------------------------------------------------*/ 91 ENTROPY_datablock(THD_datablock * blk)92double ENTROPY_datablock( THD_datablock *blk ) 93 { 94 int iv ; 95 double sum ; 96 97 ENTRY("ENTROPY_datablock") ; 98 99 ENTROPY_setup() ; 100 101 for( iv=0 ; iv < blk->nvals ; iv++ ) 102 ENTROPY_accumulate( DBLK_BRICK_BYTES(blk,iv) , DBLK_ARRAY(blk,iv) ) ; 103 104 sum = ENTROPY_compute() ; 105 ENTROPY_setdown() ; 106 RETURN(sum) ; 107 } 108