1 #include "mrilib.h"
2
main(int argc,char * argv[])3 int main( int argc , char * argv[] )
4 {
5 int iarg , nvox=0 , iv,ii,cnum , verb=1 ;
6 int automask=1 ; /* allow masks as input 14 Jul 2010 [rickr] */
7 THD_3dim_dataset *aset , *bset ;
8 byte *amm , *bmm ; int naa , nbb , nabu,nabi , naout , nbout ;
9 float paout , pbout , xrat,yrat,zrat ;
10 float_triple axyz , bxyz ;
11
12 /*-- read command line arguments --*/
13
14 if( argc < 2 || strncmp(argv[1],"-help",5) == 0 ){
15 printf(
16 "Usage: 3dABoverlap [options] A B\n"
17 "Output (to screen) is a count of various things about how\n"
18 "the automasks of datasets A and B overlap or don't overlap.\n"
19 "\n"
20 "* Dataset B will be resampled to match dataset A, if necessary,\n"
21 " which will be slow if A is high resolution. In such a case,\n"
22 " you should only use one sub-brick from dataset B.\n"
23 " ++ The resampling of B is done before the automask is generated.\n"
24 "* The values output are labeled thusly:\n"
25 " #A = number of voxels in the A mask\n"
26 " #B = number of voxels in the B mask\n"
27 " #(A uni B) = number of voxels in the either or both masks (set union)\n"
28 " #(A int B) = number of voxels present in BOTH masks (set intesection)\n"
29 " #(A \\ B) = number of voxels in A mask that aren't in B mask\n"
30 " #(B \\ A) = number of voxels in B mask that arent' in A mask\n"
31 " %%(A \\ B) = percentage of voxels from A mask that aren't in B mask\n"
32 " %%(B \\ A) = percentage of voxels from B mask that aren't in A mask\n"
33 " Rx(B/A) = radius of gyration of B mask / A mask, in x direction\n"
34 " Ry(B/A) = radius of gyration of B mask / A mask, in y direction\n"
35 " Rz(B/A) = radius of gyration of B mask / A mask, in z direction\n"
36 "* If B is an EPI dataset sub-brick, and A is a skull stripped anatomical\n"
37 " dataset, then %%(B \\ A) might be useful for assessing if the EPI\n"
38 " brick B is grossly misaligned with respect to the anatomical brick A.\n"
39 "* The radius of gyration ratios might be useful for determining if one\n"
40 " dataset is grossly larger or smaller than the other.\n"
41 "\n"
42 "OPTIONS\n"
43 "-------\n"
44 " -no_automask = consider input datasets as masks\n"
45 " (automask does not work on mask datasets)\n"
46 " -quiet = be as quiet as possible (without being entirely mute)\n"
47 " -verb = print out some progress reports (to stderr)\n"
48 "\n"
49 "NOTES\n"
50 "-----\n"
51 " * If an input dataset is comprised of bytes and contains only one\n"
52 " sub-brick, then this program assumes it is already an automask-\n"
53 " generated dataset and the automask operation will be skipped.\n"
54 ) ;
55 PRINT_COMPILE_DATE ; exit(0) ;
56 }
57
58 iarg = 1 ;
59
60 /*-- 20 Apr 2001: addto the arglist, if user wants to [RWCox] --*/
61
62 /* check options */
63
64 while( iarg < argc && argv[iarg][0] == '-' ){
65
66 if( strcmp(argv[iarg],"-verb") == 0 ){
67 verb++ ; iarg++ ; continue ;
68 }
69
70 if( strcmp(argv[iarg],"-quiet") == 0 ){
71 verb = 0 ; iarg++ ; continue ;
72 }
73
74 if( strcmp(argv[iarg],"-no_automask") == 0 ){
75 automask = 0 ; iarg++ ; continue ;
76 }
77
78 ERROR_exit("Illegal option: %s",argv[iarg]) ;
79 }
80
81 mainENTRY("3dABoverlap main") ; machdep() ;
82 if( verb ) PRINT_VERSION("3dABoverlap") ;
83 AFNI_logger("3dABoverlap",argc,argv) ;
84
85 /* input datasets */
86
87 if( iarg+1 >= argc ) ERROR_exit("Need 2 input datasets on command line") ;
88
89 aset = THD_open_dataset(argv[iarg]) ; CHECK_OPEN_ERROR(aset,argv[iarg]) ; iarg++ ;
90 bset = THD_open_dataset(argv[iarg]) ; CHECK_OPEN_ERROR(bset,argv[iarg]) ; iarg++ ;
91
92 nvox = DSET_NVOX(aset) ;
93
94 if( !EQUIV_GRIDS(aset,bset) ){ /** must resample **/
95 THD_3dim_dataset *cset ;
96 if( verb ) INFO_message("resampling dataset B to match dataset A") ;
97
98 cset = r_new_resam_dset( bset, aset, 0.0,0.0,0.0,NULL,
99 MRI_BILINEAR, NULL, 1, 0 ) ;
100 DSET_delete(bset) ; bset = cset ;
101 }
102
103 if( iarg < argc ) WARNING_message("Extra arguments?") ;
104
105 DSET_load(aset); CHECK_LOAD_ERROR(aset);
106 DSET_load(bset); CHECK_LOAD_ERROR(bset);
107
108 /* 10 Aug 2009: keep input datasets without automask, if appropriate */
109
110 if( DSET_NVALS(aset) > 1 || DSET_BRICK_TYPE(aset,0) != MRI_byte ){
111 /* allow masks as input (via -no_automask) 14 Jul 2010 [rickr] */
112 if( automask ) amm = THD_automask(aset);
113 else amm = THD_makemask(aset, 0, 1, 0); /* use any non-zero */
114 DSET_unload(aset);
115 } else {
116 amm = DSET_BRICK_ARRAY(aset,0) ;
117 }
118
119 if( DSET_NVALS(bset) > 1 || DSET_BRICK_TYPE(bset,0) != MRI_byte ){
120 /* allow masks as input (via -no_automask) 14 Jul 2010 [rickr] */
121 if( automask ) bmm = THD_automask(bset);
122 else bmm = THD_makemask(bset, 0, 1, 0); /* use any non-zero */
123 DSET_unload(bset);
124 } else {
125 bmm = DSET_BRICK_ARRAY(bset,0) ;
126 }
127
128 naa = mask_count ( nvox , amm ) ;
129 nbb = mask_count ( nvox , bmm ) ;
130 nabi = mask_intersect_count( nvox , amm , bmm ) ;
131 nabu = mask_union_count ( nvox , amm , bmm ) ;
132 naout = naa - nabi ;
133 nbout = nbb - nabi ;
134 paout = (naa > 0) ? naout/(float)naa : 0.0f ;
135 pbout = (nbb > 0) ? nbout/(float)nbb : 0.0f ;
136
137 axyz = mask_rgyrate( DSET_NX(aset),DSET_NY(aset),DSET_NZ(aset) , amm ) ;
138 bxyz = mask_rgyrate( DSET_NX(bset),DSET_NY(bset),DSET_NZ(bset) , bmm ) ;
139
140 xrat = (axyz.a > 0.0f && bxyz.a > 0.0f) ? bxyz.a / axyz.a : 0.0f ;
141 yrat = (axyz.b > 0.0f && bxyz.b > 0.0f) ? bxyz.b / axyz.b : 0.0f ;
142 zrat = (axyz.c > 0.0f && bxyz.c > 0.0f) ? bxyz.c / axyz.c : 0.0f ;
143
144 if( verb )
145 printf("#A=%s B=%s\n",DSET_BRIKNAME(aset),DSET_BRIKNAME(bset)) ;
146 if( verb )
147 printf("#A #B #(A uni B) #(A int B) "
148 "#(A \\ B) #(B \\ A) %%(A \\ B) %%(B \\ A) "
149 "Rx(B/A) Ry(B/A) Rz(B/A)\n") ;
150 printf("%-12d %-12d %-12d %-12d %-12d %-12d %7.4f %7.4f %7.4f %7.4f %7.4f\n",
151 naa , nbb , nabu, nabi, naout,nbout,100.0f*paout,100.0f*pbout,
152 xrat,yrat,zrat ) ;
153
154 exit(0) ;
155 }
156