1 #include "mrilib.h"
2 #include <string.h>
3
main(int argc,char * argv[])4 int main( int argc , char * argv[] )
5 {
6 MRI_IMAGE * imbar , * imsdev , * imwt ;
7 int ii , nvox , nsum ;
8 float * bar , * sdev , * wt ;
9 float bmax , smax , bcut , scl , bsum ;
10 float hist[11] ;
11
12 if( argc < 4 ){
13 printf("Usage: %s mean_image sdev_image wt_image\n",argv[0]) ;
14 exit(0) ;
15 }
16
17 imbar = mri_read_just_one( argv[1] ) ; if( imbar == NULL ) exit(1) ;
18 imsdev = mri_read_just_one( argv[2] ) ; if( imsdev== NULL ) exit(1) ;
19
20 if( imbar->kind != MRI_float ){
21 imwt = mri_to_float(imbar) ; mri_free(imbar) ; imbar = imwt ;
22 }
23
24 if( imsdev->kind != MRI_float ){
25 imwt = mri_to_float(imsdev) ; mri_free(imsdev) ; imsdev = imwt ;
26 }
27
28 nvox = imbar->nvox ;
29 imwt = mri_new( imbar->nx , imbar->ny , MRI_float ) ;
30
31 bar = MRI_FLOAT_PTR(imbar) ;
32 sdev = MRI_FLOAT_PTR(imsdev) ;
33 wt = MRI_FLOAT_PTR(imwt) ;
34
35 bcut = 0.05 * mri_maxabs(imbar) ; bsum = 0.0 ; nsum = 0 ;
36 for( ii=0 ; ii < nvox ; ii++ ){
37 if( bar[ii] > bcut ){ bsum += bar[ii] ; nsum++ ; }
38 }
39 bcut = 0.25 * bsum / nsum ;
40
41 smax = mri_maxabs(imsdev) ; scl = 0.25 * bcut*bcut ;
42
43 printf("cutoff value = %f\n",bcut) ;
44
45 for( ii=0 ; ii < nvox ; ii++ ){
46 if( bar[ii] < bcut || sdev[ii] <= 0.0 ) wt[ii] = 0.0 ;
47 else wt[ii] = scl / SQR(sdev[ii]) ;
48 }
49
50 mri_write( argv[3] , imwt ) ;
51 exit(0) ;
52 }
53