00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022 #include "mrilib.h"
00023 #include <stdio.h>
00024 #include <stdlib.h>
00025 #include <ctype.h>
00026 #include <math.h>
00027
00028 short non_zero[65536];
00029
00030 void Error_Exit(char *message)
00031 {
00032 fprintf(stderr, "\n\nError: %s\n", message);
00033 exit(1);
00034 }
00035
00036
00037 int main(int argc, char *argv[])
00038 {
00039 THD_3dim_dataset *mask_dset = NULL, *input_dset = NULL;
00040 int mask_subbrik = 0;
00041 int sigma = 0, nzmean = 0, nzcount = 0, debug = 0, quiet = 0, summary = 0;
00042 int minmax = 0, nzminmax = 0;
00043 short *mask_data;
00044 int nvox, i, brik;
00045 int num_ROI, ROI, mask_f2s = 0;
00046 int force_num_ROI = 0;
00047 int narg = 1;
00048 double *sum, *sumsq, *nzsum, sig, *sumallbriks;
00049 double *min, *max, *nzmin, *nzmax;
00050 long *voxels, *nzvoxels;
00051 float *input_data;
00052 byte *temp_datab;
00053 short *temp_datas;
00054
00055 if (argc < 3 || strcmp(argv[1], "-help") == 0) {
00056 printf("Usage: 3dROIstats -mask[n] mset [options] datasets\n"
00057 "Options:\n"
00058 " -mask[n] mset Means to use the dataset 'mset' as a mask:\n"
00059 " If n is present, it specifies which sub-brick\n"
00060 " in mset to use a la 3dcalc. Note: do not include\n"
00061 " the brackets if specifing a sub-brick, they are\n"
00062 " there to indicate that they are optional. If not\n"
00063 " present, 0 is assumed\n"
00064 " Voxels with the same nonzero values in 'mset'\n"
00065 " will be statisticized from 'dataset'. This will\n"
00066 " be repeated for all the different values in mset.\n"
00067 " I.e. all of the 1s in mset are one ROI, as are all\n"
00068 " of the 2s, etc.\n"
00069 " Note that the mask dataset and the input dataset\n"
00070 " must have the same number of voxels and that mset\n"
00071 " must be BYTE or SHORT (i.e., float masks won't work\n"
00072 " without the -mask_f2short option).\n"
00073 " \n"
00074 " -mask_f2short Tells the program to convert a float mask to short\n"
00075 " integers, by simple rounding. This option is needed\n"
00076 " when the mask dataset is a 1D file, for instance\n"
00077 " (since 1D files are read as floats).\n"
00078 "\n"
00079 " Be careful with this, it may not be appropriate to do!\n"
00080 "\n"
00081 " -numROI n Forces the assumption that the mask dataset's ROIs are\n"
00082 " denoted by 1 to n inclusive. Normally, the program\n"
00083 " figures out the ROIs on its own. This option is \n"
00084 " useful if a) you are certain that the mask dataset\n"
00085 " has no values outside the range [0 n], b) there may \n"
00086 " be some ROIs missing between [1 n] in the mask data-\n"
00087 " set and c) you want those columns in the output any-\n"
00088 " way so the output lines up with the output from other\n"
00089 " invocations of 3dROIstats. Confused? Then don't use\n"
00090 " this option!\n"
00091 "\n"
00092 " -debug Print out debugging information\n"
00093 " -quiet Do not print out labels for columns or rows\n"
00094 "\n"
00095 "The following options specify what stats are computed. By default\n"
00096 "the mean is always computed.\n"
00097 "\n"
00098 " -nzmean Compute the mean using only non_zero voxels. Implies\n"
00099 " the oppisite for the normal mean computed\n"
00100 " -nzvoxels Compute the number of non_zero voxels\n"
00101 " -minmax Compute the min/max of all voxels\n"
00102 " -nzminmax Compute the min/max of non_zero voxels\n"
00103 " -sigma Means to compute the standard deviation as well\n"
00104 " as the mean.\n"
00105 " -summary Only output a summary line with the grand mean across all briks\n"
00106 " in the input dataset. \n"
00107 "\n"
00108 "The output is printed to stdout (the terminal), and can be\n"
00109 "saved to a file using the usual redirection operation '>'.\n"
00110 "\n"
00111 "N.B.: The input datasets and the mask dataset can use sub-brick\n"
00112 " selectors, as detailed in the output of 3dcalc -help.\n"
00113 );
00114
00115 printf("\n" MASTER_SHORTHELP_STRING);
00116
00117 exit(0);
00118 }
00119
00120
00121
00122 machdep() ;
00123 { int new_argc ; char ** new_argv ;
00124 addto_args( argc , argv , &new_argc , &new_argv ) ;
00125 if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; }
00126 }
00127
00128
00129
00130 while (narg < argc && argv[narg][0] == '-') {
00131
00132 if (strncmp(argv[narg], "-mask_f2short", 9) == 0) {
00133 mask_f2s = 1;
00134 narg++;
00135 continue;
00136 }
00137 if (strncmp(argv[narg], "-mask", 5) == 0) {
00138 if (mask_dset != NULL)
00139 Error_Exit("Cannot have two -mask options!");
00140
00141 if (narg + 1 >= argc)
00142 Error_Exit("-mask option requires a following argument!");
00143
00144 mask_dset = THD_open_dataset(argv[++narg]);
00145
00146 if (mask_dset == NULL)
00147 Error_Exit("Cannot open mask dataset!");
00148
00149 if (DSET_BRICK_TYPE(mask_dset, 0) == MRI_complex)
00150 Error_Exit("Cannot deal with complex-valued mask dataset!");
00151
00152
00153
00154 if (isdigit(argv[narg - 1][5])) {
00155 mask_subbrik = (int) atol(argv[narg - 1] + 5);
00156 if ((mask_subbrik < 0) || (mask_subbrik >= (DSET_NVALS(mask_dset))))
00157 Error_Exit("Illegal sub-brisk following the -mask option");
00158 }
00159 narg++;
00160 continue;
00161 }
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172 if (strncmp(argv[narg], "-sigma", 5) == 0) {
00173 sigma = 1;
00174 narg++;
00175 continue;
00176 }
00177 if (strncmp(argv[narg], "-nzmean", 5) == 0) {
00178 nzmean = 1;
00179 narg++;
00180 continue;
00181 }
00182 if (strncmp(argv[narg], "-minmax", 5) == 0) {
00183 minmax = 1;
00184 narg++;
00185 continue;
00186 }
00187 if (strncmp(argv[narg], "-nzminmax", 5) == 0) {
00188 nzminmax = 1;
00189 narg++;
00190 continue;
00191 }
00192 if (strncmp(argv[narg], "-debug", 5) == 0) {
00193 debug = 1;
00194 narg++;
00195 continue;
00196 }
00197 if (strncmp(argv[narg], "-quiet", 5) == 0) {
00198 quiet = 1;
00199 narg++;
00200 continue;
00201 }
00202 if (strncmp(argv[narg], "-summary", 5) == 0) {
00203 summary = 1;
00204 narg++;
00205 continue;
00206 }
00207 if (strncmp(argv[narg], "-nzvoxels", 5) == 0) {
00208 nzcount = 1;
00209 narg++;
00210 continue;
00211 }
00212
00213 if (strncmp(argv[narg], "-numROI", 5) == 0) {
00214 force_num_ROI = (int) atol(argv[++narg]);
00215 narg++;
00216 continue;
00217 }
00218 Error_Exit("Unknown option");
00219 }
00220
00221
00222
00223 if (narg >= argc)
00224 Error_Exit("No input datasets!?\n");
00225
00226 if (mask_dset == NULL){
00227 fprintf(stderr,"Error: missing '-mask MASK_DATASET'\n");
00228 return 1;
00229 }
00230
00231
00232 if (DSET_BRICK_TYPE(mask_dset, 0) == MRI_float && mask_f2s == 0 ){
00233 fprintf(stderr,"\nError: cannot deal with float-valued mask dataset\n");
00234 fprintf(stderr,"(consider the -mask_f2short option)\n");
00235 return 1;
00236 }
00237
00238
00239
00240 for (i = 0; i < 65536; non_zero[i++] = 0);
00241
00242 DSET_load(mask_dset);
00243 if (DSET_ARRAY(mask_dset, mask_subbrik) == NULL)
00244 Error_Exit("Cannot read in mask dataset BRIK!");
00245
00246 nvox = DSET_NVOX(mask_dset);
00247 if (debug)
00248 fprintf(stderr, "The number of voxels in the mask dataset is %d\n", nvox);
00249
00250 switch (DSET_BRICK_TYPE(mask_dset, mask_subbrik)) {
00251 default:
00252 Error_Exit("Cannot deal with mask dataset datum!");
00253
00254 case MRI_short:{
00255 mask_data = (short *) DSET_ARRAY(mask_dset, mask_subbrik);
00256 for (i = 0; i < nvox; i++)
00257 if (mask_data[i]) {
00258 if (debug)
00259 fprintf(stderr, "Nonzero mask voxel %d value is %d\n", i, mask_data[i]);
00260 non_zero[mask_data[i] + 32768] = 1;
00261 }
00262 break;
00263 }
00264
00265 case MRI_byte:{
00266 byte *temp_byte = (byte *) DSET_ARRAY(mask_dset, mask_subbrik);
00267 if ((mask_data = (short *) malloc(nvox * sizeof(short))) == NULL)
00268 Error_Exit("Memory allocation error");
00269
00270 for (i = 0; i < nvox; i++) {
00271 mask_data[i] = (short) temp_byte[i];
00272 if (mask_data[i]) {
00273 if (debug)
00274 fprintf(stderr, "Nonzero mask voxel %d value is %d\n", i, mask_data[i]);
00275 non_zero[mask_data[i] + 32768] = 1;
00276 }
00277 }
00278 break;
00279 }
00280
00281 case MRI_float:{
00282 float *temp_float = (float *) DSET_ARRAY(mask_dset, mask_subbrik);
00283 if ((mask_data = (short *) malloc(nvox * sizeof(short))) == NULL)
00284 Error_Exit("Memory allocation error");
00285
00286 for (i = 0; i < nvox; i++) {
00287 mask_data[i] = (short) temp_float[i];
00288 if (mask_data[i]) {
00289 if (debug)
00290 fprintf(stderr, "Nonzero mask voxel %d value is %d\n", i, mask_data[i]);
00291 non_zero[mask_data[i] + 32768] = 1;
00292 }
00293 }
00294 break;
00295 }
00296 }
00297
00298
00299 if (!quiet && !summary)
00300 fprintf(stdout, "File\tSub-brick");
00301
00302 for (i = 0, num_ROI = 0; i < 65536; i++)
00303 if (non_zero[i]) {
00304 if (force_num_ROI && (((i - 32768) < 0) || ((i - 32768) > force_num_ROI)))
00305 Error_Exit("You used the numROI option, yet in the mask there was a\n"
00306 "value in the mask outside the range [1 n].\nMaybe you shouldn't use\n"
00307 "that option\n");
00308 non_zero[i] = num_ROI;
00309 num_ROI++;
00310 if (!quiet && !summary && !force_num_ROI) {
00311 fprintf(stdout, "\tMean_%d ", i - 32768);
00312 if (nzmean)
00313 fprintf(stdout, "\tNZMean_%d", i - 32768);
00314 if (nzcount)
00315 fprintf(stdout, "\tNZcount_%d", i - 32768);
00316 if (sigma)
00317 fprintf(stdout, "\tSigma_%d", i - 32768);
00318 if (minmax) {
00319 fprintf(stdout, "\tMin_%d ", i - 32768);
00320 fprintf(stdout, "\tMax_%d ", i - 32768);
00321 }
00322 if (nzminmax) {
00323 fprintf(stdout, "\tNZMin_%d ", i - 32768);
00324 fprintf(stdout, "\tNZMax_%d ", i - 32768);
00325 }
00326 }
00327 }
00328
00329 if (force_num_ROI) {
00330 for (i = 1; i <= force_num_ROI; i++) {
00331 non_zero[i + 32768] = (i - 1);
00332 if (!quiet && !summary) {
00333 fprintf(stdout, "\tMean_%d ", i );
00334 if (nzmean)
00335 fprintf(stdout, "\tNZMean_%d", i );
00336 if (nzcount)
00337 fprintf(stdout, "\tNZcount_%d", i );
00338 if (sigma)
00339 fprintf(stdout, "\tSigma_%d", i );
00340 if (minmax) {
00341 fprintf(stdout, "\tMin_%d ", i );
00342 fprintf(stdout, "\tMax_%d ", i );
00343 }
00344 if (nzminmax) {
00345 fprintf(stdout, "\tNZMin_%d ", i );
00346 fprintf(stdout, "\tNZMax_%d ", i );
00347 }
00348 }
00349 }
00350 num_ROI = force_num_ROI;
00351 }
00352 if (!quiet && !summary)
00353 fprintf(stdout, "\n");
00354
00355 if (debug)
00356 fprintf(stderr, "Total Number of ROIs are %d\n", num_ROI);
00357
00358
00359
00360
00361
00362 if ((sum = (double *) malloc(num_ROI * sizeof(double))) == NULL)
00363 Error_Exit("Memory allocation error");
00364 if ((voxels = (long *) malloc(num_ROI * sizeof(long))) == NULL)
00365 Error_Exit("Memory allocation error");
00366 if (nzmean || nzcount || nzminmax) {
00367 if ((nzsum = (double *) malloc(num_ROI * sizeof(double))) == NULL)
00368 Error_Exit("Memory allocation error");
00369 if ((nzvoxels = (long *) malloc(num_ROI * sizeof(long))) == NULL)
00370 Error_Exit("Memory allocation error");
00371 if ((nzmin = (double *) malloc(num_ROI * sizeof(double))) == NULL)
00372 Error_Exit("Memory allocation error");
00373 if ((nzmax = (double *) malloc(num_ROI * sizeof(double))) == NULL)
00374 Error_Exit("Memory allocation error");
00375 }
00376 if ( minmax ) {
00377 if ((min = (double *) malloc(num_ROI * sizeof(double))) == NULL)
00378 Error_Exit("Memory allocation error - min");
00379 if ((max = (double *) malloc(num_ROI * sizeof(double))) == NULL)
00380 Error_Exit("Memory allocation error - max");
00381 }
00382 if (sigma)
00383 if ((sumsq = (double *) malloc(num_ROI * sizeof(double))) == NULL)
00384 Error_Exit("Memory allocation error");
00385
00386 if (summary)
00387 if ((sumallbriks = (double *) malloc(num_ROI * sizeof(double))) == NULL)
00388 Error_Exit("Memory allocation error");
00389
00390
00391
00392 for (; narg < argc; narg++) {
00393
00394 input_dset = THD_open_dataset(argv[narg]);
00395
00396 if (input_dset == NULL) {
00397 fprintf(stderr, "Warning: Cannot open input dataset %s\n", argv[narg]);
00398 continue;
00399 }
00400 if (DSET_NVOX(input_dset) != nvox) {
00401 fprintf(stderr, "Warning: Input dataset %s is a different size than the mask\n", argv[narg]);
00402 continue;
00403 }
00404 DSET_load(input_dset);
00405
00406 if (summary)
00407 for (i = 0; i < num_ROI; i++)
00408 sumallbriks[i] = 0;
00409
00410 for (brik = 0; brik < DSET_NVALS(input_dset); brik++) {
00411
00412 for (i = 0; i < num_ROI; i++) {
00413 sum[i] = 0;
00414 voxels[i] = 0;
00415 if (nzmean || nzcount) {
00416 nzsum[i] = 0;
00417 nzvoxels[i] = 0;
00418 }
00419 if (sigma)
00420 sumsq[i] = 0;
00421 }
00422
00423 switch (DSET_BRICK_TYPE(input_dset, brik)) {
00424
00425 case MRI_byte:{
00426 float fac = DSET_BRICK_FACTOR(input_dset, brik);
00427 if (fac == 0)
00428 fac = 1.0;
00429 temp_datab = (byte *) DSET_ARRAY(input_dset, brik);
00430 if ((input_data = (float *) malloc(nvox * sizeof(float))) == NULL)
00431 Error_Exit("Memory allocation error");
00432 for (i = 0; i < nvox; i++)
00433 input_data[i] = fac * (float) temp_datab[i];
00434 break;
00435 }
00436
00437 case MRI_short:{
00438 float fac = DSET_BRICK_FACTOR(input_dset, brik);
00439 if (fac == 0)
00440 fac = 1.0;
00441 temp_datas = (short *) DSET_ARRAY(input_dset, brik);
00442 if ((input_data = (float *) malloc(nvox * sizeof(float))) == NULL)
00443 Error_Exit("Memory allocation error");
00444 for (i = 0; i < nvox; i++)
00445 input_data[i] = fac * (float) temp_datas[i];
00446 break;
00447 }
00448
00449 case MRI_float:{
00450 float fac = DSET_BRICK_FACTOR(input_dset, brik);
00451 input_data = (float *) DSET_ARRAY(input_dset, brik);
00452 if (fac == 0)
00453 fac = 1.0;
00454 else
00455 for (i = 0; i < nvox; input_data[i++] *= fac);
00456 break;
00457 }
00458
00459 default:{
00460 fprintf(stderr, "Cannot use sub-brick %d for file %s. Is it complex?\n", brik, argv[narg]);
00461 continue;
00462 }
00463
00464 }
00465
00466
00467 for (i = 0; i < num_ROI; i++) {
00468 if ( minmax ) {
00469 min[i] = 1e30;
00470 max[i] = -1e30;
00471 }
00472 if ( nzminmax ) {
00473 nzmin[i] = 1e30;
00474 nzmax[i] = -1e30;
00475 }
00476 }
00477
00478
00479
00480 for (i = 0; i < nvox; i++) {
00481 if (mask_data[i]) {
00482 ROI = non_zero[mask_data[i] + 32768];
00483 if ((ROI < 0) || (ROI >= num_ROI))
00484 Error_Exit("Somehow I boned computing how many ROIs existed");
00485
00486 if ( minmax ) {
00487 if (input_data[i] < min[ROI] ) min[ROI] = input_data[i];
00488 if (input_data[i] > max[ROI] ) max[ROI] = input_data[i];
00489 }
00490
00491 sum[ROI] += (double) input_data[i];
00492 voxels[ROI]++;
00493 if (nzmean || nzcount || nzminmax) {
00494 if (input_data[i] != 0.0) {
00495 nzsum[ROI] += (double) input_data[i];
00496 nzvoxels[ROI]++;
00497 if (input_data[i] < nzmin[ROI] )
00498 nzmin[ROI] = input_data[i];
00499 if (input_data[i] > nzmax[ROI] )
00500 nzmax[ROI] = input_data[i];
00501 }
00502 }
00503 if (sigma)
00504 sumsq[ROI] += input_data[i] * input_data[i];
00505 }
00506 }
00507
00508
00509 if (!quiet && !summary)
00510 fprintf(stdout, "%s\t%d", argv[narg], brik);
00511 if (!summary) {
00512 for (i = 0; i < num_ROI; i++) {
00513 if (voxels[i]) {
00514 fprintf(stdout, "\t%f", (float) (sum[i] / (double) voxels[i]));
00515 if (nzmean)
00516 fprintf(stdout, "\t%f", nzvoxels[i] ? (float) (nzsum[i] / (double) nzvoxels[i]) : 0.0);
00517 if (nzcount)
00518 fprintf(stdout, "\t%ld", nzvoxels[i]);
00519 if (sigma) {
00520 double mean = sum[i] / (double) voxels[i];
00521 sumsq[i] /= (double) voxels[i];
00522 if (voxels[i] == 1)
00523 sig = 1e30;
00524 else
00525 sig = sqrt((voxels[i] / (voxels[i] - 1)) * (sumsq[i] - mean * mean));
00526 fprintf(stdout, "\t%f", (float) sig);
00527 }
00528 if (minmax) {
00529 fprintf(stdout, "\t%f", min[i] );
00530 fprintf(stdout, "\t%f", max[i] );
00531 }
00532 if (nzminmax) {
00533 fprintf(stdout, "\t%f", nzmin[i] );
00534 fprintf(stdout, "\t%f", nzmax[i] );
00535 }
00536 } else {
00537 fprintf(stdout, "\t ");
00538 if (nzmean)
00539 fprintf(stdout, "\t ");
00540 if (nzcount)
00541 fprintf(stdout, "\t ");
00542 if (sigma)
00543 fprintf(stdout, "\t ");
00544 if (minmax) {
00545 fprintf(stdout, "\t ");
00546 fprintf(stdout, "\t ");
00547 }
00548 if (nzminmax) {
00549 fprintf(stdout, "\t ");
00550 fprintf(stdout, "\t ");
00551 }
00552 }
00553 }
00554
00555 fprintf(stdout, "\n");
00556 } else
00557 for (i = 0; i < num_ROI; i++) {
00558 if (nzmean)
00559 sumallbriks[i] += nzvoxels[i] ? (double) (nzsum[i] / (double) nzvoxels[i]) : 0.0;
00560 else
00561 sumallbriks[i] += (double) (sum[i] / (double) voxels[i]);
00562 }
00563
00564
00565 if (DSET_BRICK_TYPE(input_dset, brik) != MRI_float)
00566 free(input_data);
00567
00568 }
00569
00570 if (summary) {
00571 fprintf(stdout, "Average Across %i sub-briks in input dataset %s\n", DSET_NVALS(input_dset), argv[narg]);
00572 if (nzmean)
00573 fprintf(stdout, "Non-zero voxels only!\n");
00574 fprintf(stdout, "ROI\t");
00575 for (i = 0; i < num_ROI; i++)
00576 fprintf(stdout, "%i\t", i);
00577 fprintf(stdout, "\nVoxels\t");
00578 for (i = 0; i < num_ROI; i++)
00579 fprintf(stdout, "%ld\t", nzmean ? nzvoxels[i] : voxels[i]);
00580 fprintf(stdout, "\nValue\t");
00581 for (i = 0; i < num_ROI; i++)
00582 fprintf(stdout, "%f\t", (float) (sumallbriks[i] / (double) DSET_NVALS(input_dset)));
00583 fprintf(stdout, "\n");
00584 }
00585 DSET_unload(input_dset);
00586
00587 }
00588
00589 if (DSET_BRICK_TYPE(mask_dset, mask_subbrik) == MRI_byte)
00590 free(mask_data);
00591 DSET_unload(mask_dset);
00592
00593 free(sum);
00594 free(voxels);
00595 if (nzmean || nzcount) {
00596 free(nzsum);
00597 free(nzvoxels);
00598 }
00599 if (sigma)
00600 free(sumsq);
00601 exit(0);
00602 }
