1 #include "mrilib.h"
2
mri_spatial_concentration(MRI_IMAGE * im)3 float mri_spatial_concentration( MRI_IMAGE *im )
4 {
5 float xcen, ycen, zcen , wsum,ww , gval,hval , dxyz ;
6 float *imar ;
7 int nx,ny,nz,nxy , nvox , ii,jj,kk,qq , nww ;
8
9 if( im == NULL ) return 0.0f ;
10 imar = MRI_FLOAT_PTR(im) ; if( imar == NULL ) return 0.0f ;
11 nx = im->nx ; ny = im->ny ; nz = im->nz ; nvox = nx*ny*nz ; nxy = nx*ny ;
12
13 xcen = ycen = zcen = wsum = 0.0f ;
14 for( nww=qq=0 ; qq < nvox ; qq++ ){
15 ww = fabsf(imar[qq]) ; if( ww == 0.0f ) continue ;
16 ii = qq % nx ; kk = qq / nxy ; jj = (qq%nxy)/nx ;
17 wsum += ww ; xcen += ww*ii ; ycen += ww*jj ; zcen += ww*kk ; nww++ ;
18 }
19 if( wsum <= 0.0f || nww < 2 ) return 0.0f ;
20 xcen /= wsum ; ycen /= wsum ; zcen /= wsum ;
21
22 gval = hval = 0.0f ;
23 for( qq=0 ; qq < nvox ; qq++ ){
24 ww = fabsf(imar[qq]) ; if( ww == 0.0f ) continue ;
25 ii = qq % nx ; kk = qq / nxy ; jj = (qq%nxy)/nx ;
26 #if 1
27 dxyz = (xcen-ii)*(xcen-ii) + (ycen-jj)*(ycen-jj) + (zcen-kk)*(zcen-kk) ;
28 #else
29 dxyz = fabsf(xcen-ii)+fabsf(ycen-jj)+fabsf(zcen-kk) ;
30 #endif
31 gval += ww*dxyz ; hval += dxyz ;
32 }
33 gval /= wsum ; hval /= nww ;
34
35 return (hval/gval) ;
36 }
37