1 #include "mrilib.h"
2
3 /*------------------------------------------------------------------------*/
4 /*! Computes various metrics betweeen the histograms (joint and marginal)
5 of two images. At this stage, experimental -- RWCox - 25 Apr 2006.
6 --------------------------------------------------------------------------*/
7
mri_metrics(MRI_IMAGE * imp,MRI_IMAGE * imq,float * met)8 void mri_metrics( MRI_IMAGE *imp , MRI_IMAGE *imq , float *met )
9 {
10 int nvox ;
11 int *rst , *pst , *qst ;
12 byte *par, *qar ;
13 float qj,rij , rat,tmp,lrr,rm1,rp1 , fac ;
14 float esum,tsum,hsum , jsum,dsum,xsum , qsum,asum ;
15 register int ii,jj,kk ;
16 MRI_IMAGE *imqq, *impp ;
17
18 if( imp == NULL || imq == NULL ) return ;
19 if( met == NULL || imp->nvox != imq->nvox ) return ;
20
21 nvox = imp->nvox ; fac = 1.0f / nvox ;
22
23 impp = (imp->kind==MRI_byte) ? imp : mri_to_byte(imp) ;
24 imqq = (imq->kind==MRI_byte) ? imq : mri_to_byte(imq) ;
25 par = MRI_BYTE_PTR(impp) ;
26 qar = MRI_BYTE_PTR(imqq) ;
27
28 pst = (int *)calloc(256 ,sizeof(int)) ;
29 qst = (int *)calloc(256 ,sizeof(int)) ;
30 rst = (int *)calloc(256*256,sizeof(int)) ;
31
32 for( kk=0 ; kk < nvox ; kk++ ){
33 ii = par[kk] ; jj = qar[kk] ;
34 pst[ii]++ ; qst[jj]++ ; rst[ii+256*jj]++ ;
35 }
36
37 esum = tsum = hsum = 0.0f ;
38 jsum = dsum = xsum = 0.0f ;
39 qsum = asum = 0.0f ;
40 for( jj=0 ; jj < 256 ; jj++ ){
41 qj = (float)qst[jj] ;
42 if( qj > 0.0f ){
43 kk = 256*jj ; qj *= fac ;
44 for( ii=0 ; ii < 256 ; ii++ ){
45 rij = (float)rst[ii+kk] ;
46 if( rij > 0.0f ){
47 rat = (qj *(float)pst[ii]) / rij ;
48 lrr = logf(rat) ; rm1 = rat-1.0f ; rp1 = rat+1.0f ;
49
50 esum += rij * lrr ;
51 jsum += rij * lrr*rm1 ;
52 tsum += rij * rm1*rm1/rp1 ;
53 dsum += rij * (rat*lrr - rp1*logf(0.5*rp1)) ;
54 xsum += rij * rm1*rm1*rp1/rat ;
55 #if 0
56 tmp = sqrtf(rat)-1.0f ; hsum += rij * tmp*tmp ;
57 #else
58 tmp = 1.0f/sqrtf(rat)-1.0f ; hsum += rij * tmp ;
59 #endif
60
61 qsum += rij * (1.0f/rat-1.0f) ;
62 tmp = 0.5*rp1 ; asum += rij * tmp * logf(tmp/sqrtf(rat)) ;
63 }
64 }
65 }
66 }
67
68 free((void*)rst); free((void*)qst); free((void*)pst);
69 if( impp != imp ) mri_free(impp);
70 if( imqq != imq ) mri_free(imqq);
71
72 met[METRIC_KULL] = -fac * esum ;
73 met[METRIC_HELL] = fac * hsum ;
74 met[METRIC_TRIA] = fac * tsum ;
75 met[METRIC_JDIV] = fac * jsum * 0.50f ;
76 met[METRIC_JSDV] = fac * dsum * 2.00f ;
77 met[METRIC_XISQ] = fac * xsum * 0.25f ;
78 met[METRIC_XXSQ] = fac * xsum * 0.50f ;
79 met[METRIC_AGDV] = fac * asum ;
80 return ;
81 }
82