1 #include "mrilib.h"
2
main(int argc,char * argv[])3 int main( int argc , char *argv[] )
4 {
5 THD_3dim_dataset *iset , *oset ;
6 int iarg=1 , kk ;
7 int doz=0 ; /* 20 Aug 2019 */
8 int doq=0 , nqbad=0 ; /* 01 Feb 2020 */
9 int dolog2=0 ; /* 29 Jun 2021 */
10 int dolog10=0 ; /* ditto */
11 char *prefix="Pval" ;
12 MRI_IMAGE *iim , *oim ;
13
14 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
15 printf(
16 "\n"
17 "Usage: 3dPval [options] dataset\n"
18 "\n"
19 "* Converts a dataset's statistical sub-bricks to p-values.\n"
20 "\n"
21 "* Sub-bricks not internally marked as statistical volumes are unchanged.\n"
22 "\n"
23 "* However, all output volumes will be converted to float format!\n"
24 "\n"
25 "* If you wish to convert only sub-brick #3 (say) of a dataset, then\n"
26 " something like this command should do the job:\n"
27 " 3dPval -prefix Zork.nii InputDataset.nii'[3]'\n"
28 "\n"
29 "* Note that sub-bricks being marked as statistical volumes, and\n"
30 " having value-to-FDR conversion curves attached, are AFNI-only\n"
31 " ideas, and are not part of any standard, NIfTI or otherwise!\n"
32 " In other words, this program will be useless for a random dataset\n"
33 " which you download from some random non-AFNI-centric site :(\n"
34 "\n"
35 "* Also note that SMALLER p- and q-values are more 'significant', but\n"
36 " that the AFNI GUI provides interactive thresholding for values\n"
37 " ABOVE a user-chosen level, so using the GUI to threshold on a\n"
38 " p-value or q-value volume will have the opposite result to what\n"
39 " you might wish for.\n"
40 "\n"
41 "* Although the program now allows conversion of statistic values\n"
42 " to z-scores or FDR q-values, instead of p-values, you can only\n"
43 " do one type of conversion per run of 3dPval. If you want p-values\n"
44 " AND q-values, you'll have to run this program twice.\n"
45 "\n"
46 "* Finally, 'sub-brick' is AFNI jargon for a single 3D volume inside\n"
47 " a multi-volume dataset.\n"
48 "\n"
49 "Options:\n"
50 "=======\n"
51 " -zscore = Convert statistic to a z-score instead, an N(0,1) deviate\n"
52 " that represents the same p-value.\n"
53 "\n"
54 " -log2 = Convert statistic to -log2(p)\n"
55 " -log10 = Convert statistic to -log10(p)\n"
56 "\n"
57 " -qval = Convert statistic to a q-value (FDR) instead:\n"
58 " + This option only works with datasets that have\n"
59 " FDR curves inserted in their headers, which most\n"
60 " AFNI statistics programs will do. The program\n"
61 " 3drefit can also do this, with the -addFDR option.\n"
62 "\n"
63 " -prefix p = Prefix name for output file (default name is 'Pval')\n"
64 "\n"
65 "\n"
66 "AUTHOR: The Man With The Golden p < 0.000001\n"
67 ) ;
68 PRINT_COMPILE_DATE ; exit(0) ;
69 }
70
71 mainENTRY("3dPval") ; machdep() ;
72
73 while( iarg < argc && argv[iarg][0] == '-' ){
74
75 if( strcasecmp(argv[iarg],"-zscore") == 0 || strcasecmp(argv[iarg],"-zstat") == 0 ){
76 dolog10 = dolog2 = doq = 0 ; doz++ ; iarg++ ; continue ; /* 20 Aug 2019 */
77 }
78
79 if( strcasecmp(argv[iarg],"-qval") == 0 ){
80 doq++ ; dolog10 = dolog2 = doz = 0 ; iarg++ ; continue ; /* 01 Feb 2020 */
81 }
82
83 if( strcasecmp(argv[iarg],"-log2") == 0 ){
84 dolog2++ ; dolog10 = doz = doq = 0 ; iarg++ ; continue ; /* 29 Jun 2021 */
85 }
86
87 if( strcasecmp(argv[iarg],"-log10") == 0 ){
88 dolog10++ ; dolog2 = doz = doq = 0 ; iarg++ ; continue ; /* 29 Jun 2021 */
89 }
90
91 if( strcmp(argv[iarg],"-prefix") == 0 ){
92 prefix = argv[++iarg] ;
93 if( !THD_filename_ok(prefix) ) ERROR_exit("bad -prefix value") ;
94 iarg++ ; continue ;
95 }
96
97 WARNING_message("Skipping unknown option %s",argv[iarg]) ;
98 }
99
100 iset = THD_open_dataset( argv[iarg] ) ; CHECK_OPEN_ERROR(iset,argv[iarg]) ;
101 DSET_load(iset) ; ; CHECK_LOAD_ERROR(iset) ;
102
103 oset = EDIT_empty_copy(iset) ;
104
105 EDIT_dset_items( oset ,
106 ADN_prefix , prefix ,
107 ADN_datum_all , MRI_float ,
108 ADN_brick_fac , NULL ,
109 ADN_none ) ;
110
111 fprintf(stderr,"++ Processing:") ;
112 for( kk=0 ; kk < DSET_NVALS(iset) ; kk++ ){
113 iim = THD_extract_float_brick(kk,iset) ;
114 if( iim == NULL )
115 ERROR_exit("Can't get sub-brick %d of input dataset :-(",kk) ;
116 DSET_unload_one(iset,kk) ;
117
118 if( doz ){
119 oim = mri_to_zscore( iim , DSET_BRICK_STATCODE(iset,kk) ,
120 DSET_BRICK_STATAUX (iset,kk) ) ;
121 } else if( doq ){
122 floatvec *fv = DSET_BRICK_FDRCURVE(iset,kk) ;
123 oim = mri_to_qval( iim , fv ) ;
124 if( fv == NULL && FUNC_IS_STAT(DSET_BRICK_STATCODE(iset,kk)) ){
125 WARNING_message(
126 "volume %d is marked as a stat, but does not have an FDR curve :(",kk) ;
127 nqbad++ ;
128 }
129 } else {
130 oim = mri_to_pval( iim , DSET_BRICK_STATCODE(iset,kk) ,
131 DSET_BRICK_STATAUX (iset,kk) ) ;
132
133 if( oim != NULL && (dolog2 || dolog10) ){ /* 29 Jun 2021 */
134 int ii , nvox=oim->nvox ; float val , *oar=MRI_FLOAT_PTR(oim) ;
135 for( ii=0 ; ii < nvox ; ii++ ){
136 val = oar[ii] ;
137 if( val > 0.0f ){
138 oar[ii] = (dolog2) ? -log2f (val)
139 : -log10f(val) ;
140 }
141 }
142 }
143 }
144
145 if( oim == NULL ){
146 oim = iim ; fprintf(stderr,"-") ; /* new data = old data */
147 } else {
148 char *olab , nlab[128] ;
149 mri_free(iim) ; fprintf(stderr,"+") ;
150 olab = DSET_BRICK_LABEL(iset,kk) ;
151 if( doz ) sprintf(nlab,"%.116s_zstat" ,olab) ;
152 else if( doq ) sprintf(nlab,"%.116s_qval" ,olab) ;
153 else if( dolog2 ) sprintf(nlab,"%.116s_log2p" ,olab) ;
154 else if( dolog10 ) sprintf(nlab,"%.116s_log10p",olab) ;
155 else sprintf(nlab,"%.116s_pval" ,olab) ; /* default */
156 EDIT_BRICK_LABEL(oset,kk,nlab) ;
157 if( doz ) EDIT_BRICK_TO_FIZT(oset,kk) ; /* change stat code */
158 else EDIT_BRICK_TO_NOSTAT(oset,kk) ; /* erase stat code */
159 }
160 /* shove new data into dataset */
161 EDIT_substitute_brick( oset , kk , MRI_float , MRI_FLOAT_PTR(oim) ) ;
162 mri_clear_and_free(oim) ;
163 }
164 fprintf(stderr,"\n") ;
165
166 if( nqbad > 0 ){
167 WARNING_message(
168 "lack of FDR curves can be supplied by using 3drefit -addFDR") ;
169 }
170
171 DSET_write(oset) ; WROTE_DSET(oset) ;
172 exit(0) ;
173 }
174