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)18 void 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)27 void 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)34 void 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)55 double 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)82 double 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)92 double 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