00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039 #define PROGRAM_NAME "AlphaSim"
00040 #define PROGRAM_AUTHOR "B. Douglas Ward"
00041 #define PROGRAM_INITIAL "18 June 1997"
00042 #define PROGRAM_LATEST "02 December 2002"
00043
00044
00045
00046 #include "mrilib.h"
00047
00048 #define MAX_NAME_LENGTH THD_MAX_NAME
00049 #define MAX_CLUSTER_SIZE 1000
00050
00051
00052
00053
00054
00055
00056 static char * mask_filename = NULL;
00057 static byte * mask_vol = NULL;
00058 static int mask_ngood = 0;
00059
00060
00061
00062
00063
00064
00065
00066 void display_help_menu()
00067 {
00068 printf
00069 (
00070 "This program performs alpha probability simulations. \n\n"
00071 "Usage: \n"
00072 "AlphaSim \n"
00073 "-nx n1 n1 = number of voxels along x-axis \n"
00074 "-ny n2 n2 = number of voxels along y-axis \n"
00075 "-nz n3 n3 = number of voxels along z-axis \n"
00076 "-dx d1 d1 = voxel size (mm) along x-axis \n"
00077 "-dy d2 d2 = voxel size (mm) along y-axis \n"
00078 "-dz d3 d3 = voxel size (mm) along z-axis \n"
00079 "[-mask mset] Use the 0 sub-brick of dataset 'mset' as a mask \n"
00080 " to indicate which voxels to analyze (a sub-brick \n"
00081 " selector is allowed) [default = use all voxels] \n"
00082 " Note: The -mask command REPLACES the -nx, -ny, -nz,\n"
00083 " -dx, -dy, and -dz commands. \n"
00084 "[-fwhm s] s = Gaussian filter width (FWHM) \n"
00085 "[-fwhmx sx] sx = Gaussian filter width, x-axis (FWHM) \n"
00086 "[-fwhmy sy] sy = Gaussian filter width, y-axis (FWHM) \n"
00087 "[-fwhmz sz] sz = Gaussian filter width, z-axis (FWHM) \n"
00088 "[-sigma s] s = Gaussian filter width (1 sigma) \n"
00089 "[-sigmax sx] sx = Gaussian filter width, x-axis (1 sigma) \n"
00090 "[-sigmay sy] sy = Gaussian filter width, y-axis (1 sigma) \n"
00091 "[-sigmaz sz] sz = Gaussian filter width, z-axis (1 sigma) \n"
00092 "[-power] perform statistical power calculations \n"
00093 "[-ax n1] n1 = extent of active region (in voxels) along x-axis \n"
00094 "[-ay n2] n2 = extent of active region (in voxels) along y-axis \n"
00095 "[-az n3] n3 = extent of active region (in voxels) along z-axis \n"
00096 "[-zsep z] z = z-score separation between signal and noise \n"
00097 "-rmm r r = cluster connection radius (mm) \n"
00098 "-pthr p p = individual voxel threshold probability \n"
00099 "-iter n n = number of Monte Carlo simulations \n"
00100 "[-quiet] suppress screen output \n"
00101 "[-out file] file = name of output file \n"
00102 );
00103
00104 exit(0);
00105 }
00106
00107
00108
00109
00110
00111
00112 void AlphaSim_error (char * message)
00113 {
00114 fprintf (stderr, "%s Error: %s \n", PROGRAM_NAME, message);
00115 exit(1);
00116 }
00117
00118
00119
00120
00121
00122
00123
00124 void initialize_options (
00125 int * nx, int * ny, int * nz,
00126 float * dx, float * dy, float * dz,
00127 int * filter, float * sigmax, float * sigmay, float * sigmaz,
00128 int * egfw,
00129 int * power, int * ax, int * ay, int * az, float * zsep,
00130 float * rmm, float * pthr, int * niter, int * quiet,
00131 char ** outfilename)
00132
00133 {
00134 *nx = 0;
00135 *ny = 0;
00136 *nz = 0;
00137 *dx = 0.0;
00138 *dy = 0.0;
00139 *dz = 0.0;
00140 *filter = 0;
00141 *sigmax = 0.0;
00142 *sigmay = 0.0;
00143 *sigmaz = 0.0;
00144 *egfw = 0;
00145 *power = 0;
00146 *ax = 0;
00147 *ay = 0;
00148 *az = 0;
00149 *zsep = 0.0;
00150 *rmm = 0.0;
00151 *pthr = 0.0;
00152 *niter = 0;
00153 *quiet = 0;
00154 *outfilename = NULL;
00155 }
00156
00157
00158
00159
00160
00161
00162
00163 void get_options (int argc, char ** argv,
00164 int * nx, int * ny, int * nz,
00165 float * dx, float * dy, float * dz,
00166 int * filter, float * sigmax, float * sigmay, float * sigmaz,
00167 int * egfw,
00168 int * power, int * ax, int * ay, int * az, float * zsep,
00169 float * rmm, float * pthr, int * niter, int * quiet,
00170 char ** outfilename)
00171 {
00172 int nopt = 1;
00173 int ival;
00174 float fval;
00175 char message[MAX_NAME_LENGTH];
00176 int mask_nx, mask_ny, mask_nz, mask_nvox;
00177 float mask_dx, mask_dy, mask_dz;
00178
00179
00180
00181 if (argc < 2 || strncmp(argv[1], "-help", 5) == 0) display_help_menu();
00182
00183
00184
00185 AFNI_logger (PROGRAM_NAME,argc,argv);
00186
00187
00188
00189 initialize_options (nx, ny, nz, dx, dy, dz, filter, sigmax, sigmay, sigmaz,
00190 egfw, power, ax, ay, az, zsep, rmm, pthr, niter, quiet,
00191 outfilename);
00192
00193
00194 while (nopt < argc )
00195 {
00196
00197
00198 if (strncmp(argv[nopt], "-nx", 3) == 0)
00199 {
00200 nopt++;
00201 if (nopt >= argc) AlphaSim_error ("need argument after -nx ");
00202 sscanf (argv[nopt], "%d", &ival);
00203 if (ival <= 0)
00204 AlphaSim_error ("illegal argument after -nx ");
00205 *nx = ival;
00206 nopt++;
00207 continue;
00208 }
00209
00210
00211
00212 if (strncmp(argv[nopt], "-ny", 3) == 0)
00213 {
00214 nopt++;
00215 if (nopt >= argc) AlphaSim_error ("need argument after -ny ");
00216 sscanf (argv[nopt], "%d", &ival);
00217 if (ival <= 0)
00218 AlphaSim_error ("illegal argument after -ny ");
00219 *ny = ival;
00220 nopt++;
00221 continue;
00222 }
00223
00224
00225
00226 if (strncmp(argv[nopt], "-nz", 3) == 0)
00227 {
00228 nopt++;
00229 if (nopt >= argc) AlphaSim_error ("need argument after -nz ");
00230 sscanf (argv[nopt], "%d", &ival);
00231 if (ival <= 0)
00232 AlphaSim_error ("illegal argument after -nz ");
00233 *nz = ival;
00234 nopt++;
00235 continue;
00236 }
00237
00238
00239
00240 if (strncmp(argv[nopt], "-dx", 3) == 0)
00241 {
00242 nopt++;
00243 if (nopt >= argc) AlphaSim_error ("need argument after -dx ");
00244 sscanf (argv[nopt], "%f", &fval);
00245 if (fval <= 0.0)
00246 AlphaSim_error ("illegal argument after -dx ");
00247 *dx = fval;
00248 nopt++;
00249 continue;
00250 }
00251
00252
00253
00254 if (strncmp(argv[nopt], "-dy", 3) == 0)
00255 {
00256 nopt++;
00257 if (nopt >= argc) AlphaSim_error ("need argument after -dy ");
00258 sscanf (argv[nopt], "%f", &fval);
00259 if (fval <= 0.0)
00260 AlphaSim_error ("illegal argument after -dy ");
00261 *dy = fval;
00262 nopt++;
00263 continue;
00264 }
00265
00266
00267
00268 if (strncmp(argv[nopt], "-dz", 3) == 0)
00269 {
00270 nopt++;
00271 if (nopt >= argc) AlphaSim_error ("need argument after -dz ");
00272 sscanf (argv[nopt], "%f", &fval);
00273 if (fval <= 0.0)
00274 AlphaSim_error ("illegal argument after -dz ");
00275 *dz = fval;
00276 nopt++;
00277 continue;
00278 }
00279
00280
00281
00282
00283 if (strcmp(argv[nopt],"-mask") == 0 )
00284 {
00285 THD_3dim_dataset * mset ; int ii,mc ;
00286 nopt++ ;
00287 if (nopt >= argc) AlphaSim_error ("need argument after -mask!") ;
00288 mask_filename = (char *) malloc (sizeof(char) * MAX_NAME_LENGTH);
00289 if (mask_filename == NULL)
00290 AlphaSim_error ("unable to allocate memory");
00291 strcpy (mask_filename, argv[nopt]);
00292 mset = THD_open_dataset (mask_filename);
00293 if (mset == NULL) AlphaSim_error ("can't open -mask dataset!") ;
00294
00295 mask_vol = THD_makemask( mset , 0 , 1.0,0.0 ) ;
00296 if (mask_vol == NULL ) AlphaSim_error ("can't use -mask dataset!") ;
00297
00298 mask_nvox = DSET_NVOX(mset) ;
00299 mask_nx = DSET_NX(mset);
00300 mask_ny = DSET_NY(mset);
00301 mask_nz = DSET_NZ(mset);
00302 mask_dx = DSET_DX(mset);
00303 mask_dy = DSET_DY(mset);
00304 mask_dz = DSET_DZ(mset);
00305
00306 DSET_delete(mset) ;
00307
00308 for( ii=mc=0 ; ii < mask_nvox ; ii++ ) if (mask_vol[ii]) mc++ ;
00309 if (mc == 0) AlphaSim_error ("mask is all zeros!") ;
00310 printf("--- %d voxels in mask\n",mc) ;
00311 mask_ngood = mc;
00312 nopt++ ; continue ;
00313 }
00314
00315
00316
00317 if (strncmp(argv[nopt], "-fwhm", 6) == 0)
00318 {
00319 nopt++;
00320 if (nopt >= argc) AlphaSim_error ("need argument after -fwhm ");
00321 sscanf (argv[nopt], "%f", &fval);
00322 if (fval < 0.0)
00323 AlphaSim_error ("illegal argument after -fwhm ");
00324 if (fval > 0.0) *filter = 1;
00325 *sigmax = FWHM_TO_SIGMA(fval);
00326 *sigmay = FWHM_TO_SIGMA(fval);
00327 *sigmaz = FWHM_TO_SIGMA(fval);
00328 nopt++;
00329 continue;
00330 }
00331
00332
00333
00334 if (strncmp(argv[nopt], "-fwhmx", 6) == 0)
00335 {
00336 nopt++;
00337 if (nopt >= argc) AlphaSim_error ("need argument after -fwhmx ");
00338 sscanf (argv[nopt], "%f", &fval);
00339 if (fval < 0.0)
00340 AlphaSim_error ("illegal argument after -fwhmx ");
00341 if (fval > 0.0) *filter = 1;
00342 *sigmax = FWHM_TO_SIGMA(fval);
00343 nopt++;
00344 continue;
00345 }
00346
00347
00348
00349 if (strncmp(argv[nopt], "-fwhmy", 6) == 0)
00350 {
00351 nopt++;
00352 if (nopt >= argc) AlphaSim_error ("need argument after -fwhmy ");
00353 sscanf (argv[nopt], "%f", &fval);
00354 if (fval < 0.0)
00355 AlphaSim_error ("illegal argument after -fwhmy ");
00356 if (fval > 0.0) *filter = 1;
00357 *sigmay = FWHM_TO_SIGMA(fval);
00358 nopt++;
00359 continue;
00360 }
00361
00362
00363
00364 if (strncmp(argv[nopt], "-fwhmz", 6) == 0)
00365 {
00366 nopt++;
00367 if (nopt >= argc) AlphaSim_error ("need argument after -fwhmz ");
00368 sscanf (argv[nopt], "%f", &fval);
00369 if (fval < 0.0)
00370 AlphaSim_error ("illegal argument after -fwhmz ");
00371 if (fval > 0.0) *filter = 1;
00372 *sigmaz = FWHM_TO_SIGMA(fval);
00373 nopt++;
00374 continue;
00375 }
00376
00377
00378
00379 if (strncmp(argv[nopt], "-sigma", 7) == 0)
00380 {
00381 nopt++;
00382 if (nopt >= argc) AlphaSim_error ("need argument after -sigma ");
00383 sscanf (argv[nopt], "%f", &fval);
00384 if (fval < 0.0)
00385 AlphaSim_error ("illegal argument after -sigma ");
00386 if (fval > 0.0) *filter = 1;
00387 *sigmax = fval;
00388 *sigmay = fval;
00389 *sigmaz = fval;
00390 nopt++;
00391 continue;
00392 }
00393
00394
00395
00396 if (strncmp(argv[nopt], "-sigmax", 7) == 0)
00397 {
00398 nopt++;
00399 if (nopt >= argc) AlphaSim_error ("need argument after -sigmax ");
00400 sscanf (argv[nopt], "%f", &fval);
00401 if (fval < 0.0)
00402 AlphaSim_error ("illegal argument after -sigmax ");
00403 if (fval > 0.0) *filter = 1;
00404 *sigmax = fval;
00405 nopt++;
00406 continue;
00407 }
00408
00409
00410
00411 if (strncmp(argv[nopt], "-sigmay", 7) == 0)
00412 {
00413 nopt++;
00414 if (nopt >= argc) AlphaSim_error ("need argument after -sigmay ");
00415 sscanf (argv[nopt], "%f", &fval);
00416 if (fval < 0.0)
00417 AlphaSim_error ("illegal argument after -sigmay ");
00418 if (fval > 0.0) *filter = 1;
00419 *sigmay = fval;
00420 nopt++;
00421 continue;
00422 }
00423
00424
00425
00426 if (strncmp(argv[nopt], "-sigmaz", 7) == 0)
00427 {
00428 nopt++;
00429 if (nopt >= argc) AlphaSim_error ("need argument after -sigmaz ");
00430 sscanf (argv[nopt], "%f", &fval);
00431 if (fval < 0.0)
00432 AlphaSim_error ("illegal argument after -sigmaz ");
00433 if (fval > 0.0) *filter = 1;
00434 *sigmaz = fval;
00435 nopt++;
00436 continue;
00437 }
00438
00439
00440
00441 if (strncmp(argv[nopt], "-egfw", 5) == 0)
00442 {
00443 *egfw = 1;
00444 nopt++;
00445 continue;
00446 }
00447
00448
00449
00450 if (strncmp(argv[nopt], "-power", 6) == 0)
00451 {
00452 *power = 1;
00453 nopt++;
00454 continue;
00455 }
00456
00457
00458
00459 if (strncmp(argv[nopt], "-ax", 3) == 0)
00460 {
00461 nopt++;
00462 if (nopt >= argc) AlphaSim_error ("need argument after -ax ");
00463 sscanf (argv[nopt], "%d", &ival);
00464 if (ival <= 0)
00465 AlphaSim_error ("illegal argument after -ax ");
00466 *ax = ival;
00467 nopt++;
00468 continue;
00469 }
00470
00471
00472
00473 if (strncmp(argv[nopt], "-ay", 3) == 0)
00474 {
00475 nopt++;
00476 if (nopt >= argc) AlphaSim_error ("need argument after -ay ");
00477 sscanf (argv[nopt], "%d", &ival);
00478 if (ival <= 0)
00479 AlphaSim_error ("illegal argument after -ay ");
00480 *ay = ival;
00481 nopt++;
00482 continue;
00483 }
00484
00485
00486
00487 if (strncmp(argv[nopt], "-az", 3) == 0)
00488 {
00489 nopt++;
00490 if (nopt >= argc) AlphaSim_error ("need argument after -az ");
00491 sscanf (argv[nopt], "%d", &ival);
00492 if (ival <= 0)
00493 AlphaSim_error ("illegal argument after -az ");
00494 *az = ival;
00495 nopt++;
00496 continue;
00497 }
00498
00499
00500
00501 if (strncmp(argv[nopt], "-zsep", 5) == 0)
00502 {
00503 nopt++;
00504 if (nopt >= argc) AlphaSim_error ("need argument after -zsep ");
00505 sscanf (argv[nopt], "%f", &fval);
00506 if (fval < 0.0)
00507 AlphaSim_error ("illegal argument after -zsep ");
00508 *zsep = fval;
00509 nopt++;
00510 continue;
00511 }
00512
00513
00514
00515 if (strncmp(argv[nopt], "-rmm", 4) == 0)
00516 {
00517 nopt++;
00518 if (nopt >= argc) AlphaSim_error ("need argument after -rmm ");
00519 sscanf (argv[nopt], "%f", &fval);
00520 if (fval <= 0.0)
00521 AlphaSim_error ("illegal argument after -rmm ");
00522 *rmm = fval;
00523 nopt++;
00524 continue;
00525 }
00526
00527
00528
00529 if (strncmp(argv[nopt], "-pthr", 5) == 0)
00530 {
00531 nopt++;
00532 if (nopt >= argc) AlphaSim_error ("need argument after -pthr ");
00533 sscanf (argv[nopt], "%f", &fval);
00534 if ((fval <= 0.0) || (fval > 1.0))
00535 AlphaSim_error ("illegal argument after -pthr ");
00536 *pthr = fval;
00537 nopt++;
00538 continue;
00539 }
00540
00541
00542
00543 if (strncmp(argv[nopt], "-iter", 5) == 0)
00544 {
00545 nopt++;
00546 if (nopt >= argc) AlphaSim_error ("need argument after -iter ");
00547 sscanf (argv[nopt], "%d", &ival);
00548 if (ival <= 0)
00549 AlphaSim_error ("illegal argument after -iter ");
00550 *niter = ival;
00551 nopt++;
00552 continue;
00553 }
00554
00555
00556
00557 if (strncmp(argv[nopt], "-quiet", 6) == 0)
00558 {
00559 *quiet = 1;
00560 nopt++;
00561 continue;
00562 }
00563
00564
00565
00566 if (strncmp(argv[nopt], "-out", 4) == 0)
00567 {
00568 nopt++;
00569 if (nopt >= argc) AlphaSim_error ("need argument after -out ");
00570 *outfilename = AFMALL( char, sizeof(char) * MAX_NAME_LENGTH);
00571 strcpy (*outfilename, argv[nopt]);
00572 nopt++;
00573 continue;
00574 }
00575
00576
00577
00578 AlphaSim_error ("unrecognized command line option ");
00579 }
00580
00581
00582
00583 if (mask_vol != NULL)
00584 {
00585 *nx = mask_nx; *ny = mask_ny; *nz = mask_nz;
00586 *dx = fabs(mask_dx); *dy = fabs(mask_dy); *dz = fabs(mask_dz);
00587 }
00588
00589 }
00590
00591
00592
00593
00594
00595
00596
00597 void check_for_valid_inputs (int nx, int ny, int nz,
00598 float dx, float dy, float dz, int filter,
00599 float sigmax, float sigmay, float sigmaz,
00600 int power, int ax, int ay, int az, float zsep,
00601 float rmm, float pthr, int niter,
00602 char * outfilename)
00603 {
00604 FILE * fout = NULL;
00605 char message[MAX_NAME_LENGTH];
00606
00607 if (nx <= 0) AlphaSim_error ("Illegal value for nx ");
00608 if (ny <= 0) AlphaSim_error ("Illegal value for ny ");
00609 if (nz <= 0) AlphaSim_error ("Illegal value for nz ");
00610 if (dx <= 0.0) AlphaSim_error ("Illegal value for dx ");
00611 if (dy <= 0.0) AlphaSim_error ("Illegal value for dy ");
00612 if (dz <= 0.0) AlphaSim_error ("Illegal value for dz ");
00613 if (filter && ((sigmax <= 0.0) || (sigmay <= 0.0) || (sigmaz <= 0.0)))
00614 AlphaSim_error ("Illegal Gaussian filter width specification ");
00615 if (power && ((ax <= 0) || (ay <= 0) || (az <= 0)))
00616 AlphaSim_error ("Illegal dimensions for activation region ");
00617 if (power && (zsep <= 0.0)) AlphaSim_error ("Illegal value for zsep ");
00618 if ( (rmm < dx) && (rmm < dy) && (rmm < dz) )
00619 AlphaSim_error ("Cluster connection radius is too small ");
00620 if ((pthr <= 0.0) || (pthr > 1.0))
00621 AlphaSim_error ("Illegal value for pthr ");
00622 if (niter <= 0) AlphaSim_error ("Illegal value for niter ");
00623
00624 if (outfilename != NULL)
00625 {
00626
00627 fout = fopen (outfilename, "r");
00628 if (fout != NULL)
00629 {
00630 sprintf (message, "Output file %s already exists. ", outfilename);
00631 AlphaSim_error (message);
00632 }
00633 }
00634
00635 }
00636
00637
00638
00639
00640
00641
00642
00643 void initialize (int argc, char ** argv,
00644 int * nx, int * ny, int * nz,
00645 float * dx, float * dy, float * dz,
00646 int * filter, float * sigmax, float * sigmay, float * sigmaz,
00647 int * egfw, float * avgsx, float * avgsy, float * avgsz,
00648 int * power, int * ax, int * ay, int * az, float * zsep,
00649 float * rmm, float * pthr, int * niter, int * quiet,
00650 char ** outfilename, long * count,
00651 double * sum, double * sumsq, float * power_thr,
00652 float ** fim, float ** arfim,
00653 long ** freq_table, long ** max_table)
00654
00655 {
00656 int i;
00657 int nxyz;
00658
00659 int which;
00660 double p, q, z, mean, sd;
00661 int status;
00662 double bound;
00663
00664
00665
00666 get_options(argc, argv,
00667 nx, ny, nz, dx, dy, dz, filter, sigmax, sigmay, sigmaz,
00668 egfw, power, ax, ay, az, zsep, rmm, pthr, niter, quiet,
00669 outfilename);
00670
00671
00672
00673 check_for_valid_inputs (*nx, *ny, *nz, *dx, *dy, *dz, *filter,
00674 *sigmax, *sigmay, *sigmaz,
00675 *power, *ax, *ay, *az, *zsep,
00676 *rmm, *pthr, *niter, *outfilename);
00677
00678
00679
00680 nxyz = (*nx) * (*ny) * (*nz);
00681 *fim = (float *) malloc(nxyz * sizeof(float));
00682 if (*fim == NULL)
00683 AlphaSim_error ("memory allocation error");
00684
00685
00686
00687 if (*power)
00688 {
00689 nxyz = (*ax) * (*ay) * (*az);
00690 *arfim = (float *) malloc(nxyz * sizeof(float));
00691 if (*arfim == NULL)
00692 AlphaSim_error ("memory allocation error");
00693 }
00694
00695
00696
00697 *freq_table = (long *) malloc( MAX_CLUSTER_SIZE * sizeof(long) );
00698 if (*freq_table == NULL)
00699 AlphaSim_error ("memory allocation error");
00700 for (i = 0; i < MAX_CLUSTER_SIZE; i++)
00701 (*freq_table)[i] = 0;
00702
00703
00704
00705 *max_table = (long *) malloc( MAX_CLUSTER_SIZE * sizeof(long) );
00706 if (*max_table == NULL)
00707 AlphaSim_error ("memory allocation error");
00708 for (i = 0; i < MAX_CLUSTER_SIZE; i++)
00709 (*max_table)[i] = 0;
00710
00711
00712
00713 *count = 0;
00714 *sum = 0.0;
00715 *sumsq = 0.0;
00716
00717
00718
00719 if (*power)
00720 {
00721 which = 2;
00722 q = *pthr;
00723 p = 1.0 - q;
00724 mean = 0.0;
00725 sd = 1.0;
00726 cdfnor (&which, &p, &q, &z, &mean, &sd, &status, &bound);
00727
00728 z = *zsep - z;
00729 which = 1;
00730 cdfnor (&which, &p, &q, &z, &mean, &sd, &status, &bound);
00731 *power_thr = p;
00732 }
00733
00734
00735
00736 if (*egfw)
00737 {
00738 *avgsx = 0.0; *avgsy = 0.0; *avgsz = 0.0;
00739 }
00740
00741
00742 srand48 (1234567);
00743
00744 }
00745
00746
00747
00748
00749
00750
00751
00752 float uniform ()
00753 {
00754 return ( (float)drand48() );
00755 }
00756
00757
00758
00759
00760
00761
00762
00763 void normal (float * n1, float * n2)
00764 {
00765 float u1, u2;
00766 float r;
00767
00768
00769 u1 = 0.0;
00770 while (u1 <= 0.0)
00771 {
00772 u1 = uniform();
00773 }
00774 u2 = uniform();
00775
00776 r = sqrt(-2.0*log(u1));
00777 *n1 = r * cos(2.0*PI*u2);
00778 *n2 = r * sin(2.0*PI*u2);
00779 }
00780
00781
00782
00783
00784
00785
00786
00787 void activation_region (int nx, int ny, int nz, int ax, int ay, int az,
00788 int * xbot, int * xtop, int * ybot, int * ytop,
00789 int * zbot, int * ztop)
00790 {
00791 *xbot = nx/2 - ax/2;
00792 *xtop = *xbot + ax;
00793 *ybot = ny/2 - ay/2;
00794 *ytop = *ybot + ay;
00795 *zbot = nz/2 - az/2;
00796 *ztop = *zbot + az;
00797 }
00798
00799
00800
00801
00802
00803
00804
00805 void generate_image (int nx, int ny, int nz, int power,
00806 int ax, int ay, int az, float zsep, float * fim)
00807 {
00808 int nxy, nxyz;
00809 int nxyzdiv2;
00810 int ixyz;
00811 float n1, n2;
00812 int xbot, xtop, ybot, ytop, zbot, ztop;
00813 int ix, jy, kz;
00814
00815
00816
00817 nxy = nx * ny;
00818 nxyz = nxy * nz;
00819 nxyzdiv2 = nxyz / 2;
00820
00821
00822 for (ixyz = 0; ixyz < nxyzdiv2; ixyz++)
00823 {
00824 normal(&n1, &n2);
00825 fim[ixyz] = n1;
00826 fim[ixyz+nxyzdiv2] = n2;
00827 }
00828 normal(&n1, &n2);
00829 fim[nxyz-1] = n1;
00830
00831
00832 if (power)
00833 {
00834
00835 activation_region (nx, ny, nz, ax, ay, az,
00836 &xbot, &xtop, &ybot, &ytop, &zbot, &ztop);
00837
00838
00839 for (ixyz = 0; ixyz < nxyz; ixyz++)
00840 {
00841 IJK_TO_THREE (ixyz, ix, jy, kz, nx, nxy);
00842 if ( (ix >= xbot) && (ix < xtop)
00843 && (jy >= ybot) && (jy < ytop)
00844 && (kz >= zbot) && (kz < ztop) )
00845 fim[ixyz] += zsep;
00846 }
00847 }
00848 }
00849
00850
00851
00852
00853
00854
00855
00856 void gaussian_filter (int nx, int ny, int nz, float dx, float dy, float dz,
00857 float rmm, float sigmax, float sigmay, float sigmaz,
00858 float * fim)
00859 {
00860
00861
00862 EDIT_blur_volume_3d (nx, ny, nz, dx, dy, dz,
00863 MRI_float, fim, sigmax, sigmay, sigmaz);
00864
00865 }
00866
00867
00868
00869
00870
00871
00872
00873 void estimate_gfw (int nx, int ny, int nz, float dx, float dy, float dz,
00874 int niter, int quiet, float * fim,
00875 float * avgsx, float * avgsy, float * avgsz)
00876 {
00877 int nxy, nxyz;
00878 int ixyz;
00879 int ix, jy, kz, ixyz2;
00880 float fsum, fsq, var;
00881 float dfdx, dfdxsum, dfdxsq, varxx;
00882 float dfdy, dfdysum, dfdysq, varyy;
00883 float dfdz, dfdzsum, dfdzsq, varzz;
00884 int countx, county, countz;
00885 float sx, sy, sz;
00886 float arg;
00887
00888
00889
00890 nxy = nx * ny;
00891 nxyz = nxy * nz;
00892
00893
00894
00895 fsum = 0.0;
00896 fsq = 0.0;
00897 for (ixyz = 0; ixyz < nxyz; ixyz++)
00898 {
00899 fsum += fim[ixyz];
00900 fsq += fim[ixyz] * fim[ixyz];
00901 }
00902 var = (fsq - (fsum * fsum)/nxyz) / (nxyz-1);
00903
00904
00905
00906 dfdxsum = 0.0; dfdysum = 0.0; dfdzsum = 0.0;
00907 dfdxsq = 0.0; dfdysq = 0.0; dfdzsq = 0.0;
00908 countx = 0; county = 0; countz = 0;
00909 for (ixyz = 0; ixyz < nxyz; ixyz++)
00910 {
00911 IJK_TO_THREE (ixyz, ix, jy, kz, nx, nxy);
00912
00913 if (ix+1 < nx)
00914 {
00915 ixyz2 = THREE_TO_IJK (ix+1, jy, kz, nx, nxy);
00916 dfdx = (fim[ixyz2] - fim[ixyz]) / 1.0;
00917 dfdxsum += dfdx;
00918 dfdxsq += dfdx * dfdx;
00919 countx += 1;
00920 }
00921
00922 if (jy+1 < ny)
00923 {
00924 ixyz2 = THREE_TO_IJK (ix, jy+1, kz, nx, nxy);
00925 dfdy = (fim[ixyz2] - fim[ixyz]) / 1.0;
00926 dfdysum += dfdy;
00927 dfdysq += dfdy * dfdy;
00928 county += 1;
00929 }
00930
00931 if (kz+1 < nz)
00932 {
00933 ixyz2 = THREE_TO_IJK (ix, jy, kz+1, nx, nxy);
00934 dfdz = (fim[ixyz2] - fim[ixyz]) / 1.0;
00935 dfdzsum += dfdz;
00936 dfdzsq += dfdz * dfdz;
00937 countz += 1;
00938 }
00939
00940 }
00941
00942
00943 if (countx < 2)
00944 varxx = 0.0;
00945 else
00946 varxx = (dfdxsq - (dfdxsum * dfdxsum)/countx) / (countx-1);
00947
00948 if (county < 2)
00949 varyy = 0.0;
00950 else
00951 varyy = (dfdysq - (dfdysum * dfdysum)/county) / (county-1);
00952
00953 if (countz < 2)
00954 varzz = 0.0;
00955 else
00956 varzz = (dfdzsq - (dfdzsum * dfdzsum)/countz) / (countz-1);
00957
00958
00959
00960 arg = 1.0 - 0.5*(varxx/var);
00961 if ( (arg <= 0.0) || (varxx == 0.0) )
00962 sx = 0.0;
00963 else
00964 sx = sqrt( -1.0 / (4.0*log(arg)) ) * dx;
00965
00966 arg = 1.0 - 0.5*(varyy/var);
00967 if ( (arg <= 0.0) || (varyy == 0.0) )
00968 sy = 0.0;
00969 else
00970 sy = sqrt( -1.0 / (4.0*log(arg)) ) * dy;
00971
00972 arg = 1.0 - 0.5*(varzz/var);
00973 if ( (arg <= 0.0) || (varzz == 0.0) )
00974 sz = 0.0;
00975 else
00976 sz = sqrt( -1.0 / (4.0*log(arg)) ) * dz;
00977
00978
00979 *avgsx += sx / niter;
00980 *avgsy += sy / niter;
00981 *avgsz += sz / niter;
00982
00983
00984 if (!quiet)
00985 {
00986 printf ("var =%f \n", var);
00987 printf ("varxx=%f varyy=%f varzz=%f \n", varxx, varyy, varzz);
00988 printf (" sx=%f sy=%f sz=%f \n", sx, sy, sz);
00989 }
00990 }
00991
00992
00993
00994
00995
00996
00997
00998 void get_activation_region (int nx, int ny, int nz, int ax, int ay, int az,
00999 float pthr, float zsep, float * fim, float * arfim)
01000 {
01001 int nxy, nxyz;
01002 int axy, axyz;
01003 int ixyz, ixyz2;
01004 int xbot, xtop, ybot, ytop, zbot, ztop;
01005 int ix, jy, kz;
01006
01007
01008
01009 nxy = nx * ny;
01010 nxyz = nxy * nz;
01011 axy = ax * ay;
01012 axyz = axy * az;
01013
01014
01015
01016 activation_region (nx, ny, nz, ax, ay, az,
01017 &xbot, &xtop, &ybot, &ytop, &zbot, &ztop);
01018
01019
01020 for (ixyz = 0; ixyz < axyz; ixyz++)
01021 {
01022 IJK_TO_THREE (ixyz, ix, jy, kz, ax, axy);
01023 ix += xbot;
01024 jy += ybot;
01025 kz += zbot;
01026 ixyz2 = THREE_TO_IJK (ix, jy, kz, nx, nxy);
01027 arfim[ixyz] = fim[ixyz2];
01028 }
01029
01030
01031 }
01032
01033
01034
01035
01036
01037
01038
01039
01040 float pcalc (int nx, int ny, int nz, float * fim, float zthr)
01041 {
01042 int nxyz;
01043 int ixyz;
01044 int pcount;
01045 float p;
01046
01047
01048
01049 nxyz = nx * ny * nz;
01050
01051 pcount = 0;
01052 for (ixyz = 0; ixyz < nxyz; ixyz++)
01053 if (fim[ixyz] > zthr) pcount ++;
01054 p = (float)pcount / (float)nxyz;
01055
01056 return (p);
01057 }
01058
01059
01060
01061
01062
01063
01064
01065 void threshold_data (int nx, int ny, int nz, float * fim,
01066 float pthr, long * count, double * sum, double * sumsq,
01067 int quiet, int iter)
01068 {
01069 const float EPSILON = 1.0e-8;
01070 int ixyz;
01071 int nxyz;
01072 float zthr;
01073 float pact;
01074
01075 int which;
01076 double p, q, z, mean, sd;
01077 int status;
01078 double bound;
01079
01080
01081
01082 nxyz = nx * ny * nz;
01083
01084
01085
01086 if (*count < 1.0e+09)
01087 {
01088 *count += nxyz;
01089 for (ixyz = 0; ixyz < nxyz; ixyz++)
01090 {
01091 *sum += fim[ixyz];
01092 *sumsq += fim[ixyz] * fim[ixyz];
01093 }
01094 }
01095
01096
01097
01098 which = 2;
01099 p = 1.0 - pthr;
01100 q = pthr;
01101 mean = (*sum) / (*count);
01102 sd = sqrt(((*sumsq) - ((*sum) * (*sum))/(*count)) / ((*count)-1));
01103 cdfnor (&which, &p, &q, &z, &mean, &sd, &status, &bound);
01104 zthr = z;
01105
01106 if (!quiet)
01107 {
01108 pact = pcalc (nx, ny, nz, fim, zthr);
01109 printf ("pthr=%f zthr=%f pact=%f ", pthr, zthr, pact);
01110 }
01111
01112
01113
01114 for (ixyz = 0; ixyz < nxyz; ixyz++)
01115 if (fim[ixyz] > zthr)
01116 fim[ixyz] = 1.0;
01117 else
01118 fim[ixyz] = 0.0;
01119
01120
01121 }
01122
01123
01124
01125
01126
01127
01128
01129 void apply_mask (int nx, int ny, int nz, float * fim)
01130
01131 {
01132 int ixyz;
01133 int nxyz;
01134
01135
01136
01137 nxyz = nx * ny * nz;
01138
01139
01140
01141 for (ixyz = 0; ixyz < nxyz; ixyz++)
01142 if (! mask_vol[ixyz])
01143 fim[ixyz] = 0.0;
01144
01145
01146 }
01147
01148
01149
01150
01151
01152
01153
01154 void identify_clusters (int nx, int ny, int nz,
01155 float dx, float dy, float dz,
01156 float rmm, float * fim, int quiet,
01157 long * freq_table, long * max_table)
01158
01159
01160
01161
01162
01163
01164
01165
01166
01167
01168
01169 {
01170 MCW_cluster_array * clar;
01171 MCW_cluster * cl;
01172 int nxy;
01173 int nxyz;
01174 int iclu, ipt;
01175 int size, max_size;
01176 int count;
01177
01178
01179
01180 nxy = nx * ny;
01181
01182
01183 clar = MCW_find_clusters( nx,ny,nz , dx,dy,dz ,
01184 MRI_float , fim , rmm ) ;
01185
01186
01187 if ((clar == NULL) || (clar->num_clu == 0))
01188 {
01189 if (!quiet)
01190 printf ("NumCl=%4d MaxClSz=%4d\n", 0, 0);
01191 if (clar != NULL) DESTROY_CLARR(clar);
01192 }
01193 else
01194 {
01195 max_size = 0;
01196 for (iclu = 0; iclu < clar->num_clu; iclu++)
01197 {
01198 cl = clar->clar[iclu] ;
01199 if( cl == NULL ) continue ;
01200
01201 size = cl->num_pt;
01202
01203 if (size < MAX_CLUSTER_SIZE)
01204 freq_table[size]++;
01205 else
01206 freq_table[MAX_CLUSTER_SIZE-1]++;
01207
01208 if (size > max_size)
01209 max_size = size;
01210
01211 }
01212
01213 if (max_size < MAX_CLUSTER_SIZE)
01214 max_table[max_size]++;
01215 else
01216 max_table[MAX_CLUSTER_SIZE-1]++;
01217
01218 if (!quiet)
01219 printf ("NumCl=%4d MaxClSz=%4d\n", clar->num_clu, max_size);
01220
01221 DESTROY_CLARR(clar);
01222 }
01223
01224 }
01225
01226
01227
01228
01229
01230
01231
01232 void output_results (int nx, int ny, int nz, float dx, float dy, float dz,
01233 int filter, float sigmax, float sigmay, float sigmaz,
01234 int egfw, float avgsx, float avgsy, float avgsz,
01235 int power, int ax, int ay, int az, float zsep,
01236 float rmm, float pthr, int niter, char * outfilename,
01237 long * freq_table, long * max_table)
01238 {
01239 const float EPSILON = 1.0e-6;
01240 int i, j;
01241 float divisor;
01242 float * prob_table;
01243 float * alpha_table;
01244 float * cum_prop_table;
01245 long total_num_clusters;
01246 char message[MAX_NAME_LENGTH];
01247 FILE * fout;
01248
01249
01250
01251 prob_table = (float *) malloc( MAX_CLUSTER_SIZE * sizeof(float) );
01252 if (prob_table == NULL)
01253 AlphaSim_error ("memory allocation error");
01254 for (i = 1; i < MAX_CLUSTER_SIZE; i++)
01255 prob_table[i] = 0.0;
01256
01257
01258 alpha_table = (float *) malloc( MAX_CLUSTER_SIZE * sizeof(float) );
01259 if (alpha_table == NULL)
01260 AlphaSim_error ("memory allocation error");
01261 for (i = 1; i < MAX_CLUSTER_SIZE; i++)
01262 alpha_table[i] = 0.0;
01263
01264
01265 cum_prop_table = (float *) malloc( MAX_CLUSTER_SIZE * sizeof(float) );
01266 if (cum_prop_table == NULL)
01267 AlphaSim_error ("memory allocation error");
01268 for (i = 1; i < MAX_CLUSTER_SIZE; i++)
01269 cum_prop_table[i] = 0.0;
01270
01271 total_num_clusters = 0;
01272 for (i = 1; i < MAX_CLUSTER_SIZE; i++)
01273 total_num_clusters += freq_table[i];
01274
01275 if (power)
01276 divisor = (float)(niter) * ax * ay * az;
01277 else
01278 if (mask_vol)
01279 divisor = (float)(niter) * mask_ngood;
01280 else
01281 divisor = (float)(niter) * nx * ny * nz;
01282
01283 for (i = 1; i < MAX_CLUSTER_SIZE; i++)
01284 {
01285 prob_table[i] = i * freq_table[i] / divisor;
01286 alpha_table[i] = (float)max_table[i] / (float)niter;
01287 cum_prop_table[i] = (float)freq_table[i] / (float)total_num_clusters;
01288 }
01289
01290 for (i = 1; i < MAX_CLUSTER_SIZE-1; i++)
01291 {
01292 j = MAX_CLUSTER_SIZE - i;
01293 prob_table[j-1] += prob_table[j];
01294 alpha_table[j-1] += alpha_table[j];
01295 cum_prop_table[i+1] += cum_prop_table[i];
01296 }
01297
01298
01299
01300 if (outfilename == NULL)
01301 fout = stdout;
01302 else
01303 {
01304
01305 fout = fopen (outfilename, "r");
01306 if (fout != NULL)
01307 {
01308 sprintf (message, "file %s already exists. ", outfilename);
01309 AlphaSim_error (message);
01310 }
01311
01312
01313 fout = fopen (outfilename, "w");
01314 if (fout == NULL)
01315 {
01316 AlphaSim_error ("unable to write file ");
01317 }
01318 }
01319
01320
01321 fprintf (fout, "\n\n");
01322 fprintf (fout, "Program: %s \n", PROGRAM_NAME);
01323 fprintf (fout, "Author: %s \n", PROGRAM_AUTHOR);
01324 fprintf (fout, "Initial Release: %s \n", PROGRAM_INITIAL);
01325 fprintf (fout, "Latest Revision: %s \n", PROGRAM_LATEST);
01326 fprintf (fout, "\n");
01327
01328 fprintf (fout, "Data set dimensions: \n");
01329 fprintf (fout, "nx = %5d ny = %5d nz = %5d (voxels)\n", nx, ny, nz);
01330 fprintf (fout, "dx = %5.2f dy = %5.2f dz = %5.2f (mm)\n", dx, dy, dz);
01331
01332 if (mask_vol)
01333 fprintf (fout, "\nMask filename = %s \n", mask_filename);
01334 if (mask_vol && !power)
01335 fprintf (fout, "Voxels in mask = %5d \n", mask_ngood);
01336
01337 fprintf (fout, "\nGaussian filter widths: \n");
01338 fprintf (fout, "sigmax = %5.2f FWHMx = %5.2f \n",
01339 sigmax, sigmax * 2.0*sqrt(2.0*log(2.0)));
01340 fprintf (fout, "sigmay = %5.2f FWHMy = %5.2f \n",
01341 sigmay, sigmay * 2.0*sqrt(2.0*log(2.0)));
01342 fprintf (fout, "sigmaz = %5.2f FWHMz = %5.2f \n\n",
01343 sigmaz, sigmaz * 2.0*sqrt(2.0*log(2.0)));
01344
01345 if (egfw)
01346 {
01347 fprintf (fout, "Estimated Gaussian filter widths: \n");
01348 fprintf (fout, "Ave sx = %f Ave sy = %f Ave sz = %f \n\n",
01349 avgsx, avgsy, avgsz);
01350 }
01351
01352 if (power)
01353 {
01354 fprintf (fout, "Activation Region for Power Calculations: \n");
01355 fprintf (fout, "ax = %5d ay = %5d az = %5d (voxels) \n",
01356 ax, ay, az);
01357 fprintf (fout, "z separation = %f \n\n", zsep);
01358 }
01359
01360 fprintf (fout, "Cluster connection radius: rmm = %5.2f \n\n", rmm);
01361 fprintf (fout, "Threshold probability: pthr = %e \n\n", pthr);
01362 fprintf (fout, "Number of Monte Carlo iterations = %5d \n\n", niter);
01363 if (!power)
01364 fprintf (fout, "Cl Size Frequency Cum Prop p/Voxel"
01365 " Max Freq Alpha\n");
01366 else
01367 fprintf (fout, "Cl Size Frequency Cum Prop p/Voxel"
01368 " Max Freq Power\n");
01369 for (i = 1; i < MAX_CLUSTER_SIZE; i++)
01370 if (alpha_table[i] < EPSILON)
01371 break;
01372 else
01373 fprintf (fout, "%7d %12ld %10.6f %10.8f %7ld %10.6f\n",
01374 i, freq_table[i], cum_prop_table[i], prob_table[i],
01375 max_table[i], alpha_table[i]);
01376
01377 fclose(fout);
01378
01379 }
01380
01381
01382
01383
01384
01385
01386
01387 void terminate (float ** fim, float ** arfim,
01388 long ** freq_table, long ** max_table)
01389 {
01390 if (*fim != NULL)
01391 { free (*fim); *fim = NULL; }
01392
01393 if (*arfim != NULL)
01394 { free (*arfim); *arfim = NULL; }
01395
01396 if (*freq_table != NULL)
01397 { free (*freq_table); *freq_table = NULL; }
01398
01399 if (*max_table != NULL)
01400 { free (*max_table); *max_table = NULL; }
01401 }
01402
01403
01404
01405
01406
01407
01408
01409 int main (int argc, char ** argv)
01410 {
01411 int nx;
01412 int ny;
01413 int nz;
01414 float dx;
01415 float dy;
01416 float dz;
01417 int filter;
01418 float sigmax;
01419 float sigmay;
01420 float sigmaz;
01421 int egfw;
01422 float avgsx;
01423 float avgsy;
01424 float avgsz;
01425 int power;
01426 int ax;
01427 int ay;
01428 int az;
01429 float zsep;
01430 float rmm;
01431 float pthr;
01432 int niter;
01433 int quiet;
01434 char * outfilename;
01435
01436 long count;
01437 double sum, sumsq;
01438 int iter;
01439 float power_thr;
01440
01441 float * fim = NULL;
01442 float * arfim = NULL;
01443 long * freq_table = NULL;
01444 long * max_table = NULL;
01445
01446
01447
01448 printf ("\n\n");
01449 printf ("Program: %s \n", PROGRAM_NAME);
01450 printf ("Author: %s \n", PROGRAM_AUTHOR);
01451 printf ("Initial Release: %s \n", PROGRAM_INITIAL);
01452 printf ("Latest Revision: %s \n", PROGRAM_LATEST);
01453 printf ("\n");
01454
01455
01456
01457 initialize (argc, argv,
01458 &nx, &ny, &nz, &dx, &dy, &dz, &filter, &sigmax, &sigmay, &sigmaz,
01459 &egfw, &avgsx, &avgsy, &avgsz, &power, &ax, &ay, &az, &zsep,
01460 &rmm, &pthr, &niter, &quiet, &outfilename, &count, &sum, &sumsq,
01461 &power_thr, &fim, &arfim, &freq_table, &max_table);
01462
01463
01464
01465 for (iter = 1; iter <= niter; iter++)
01466 {
01467 if (!quiet) printf ("Iter =%5d \n", iter);
01468
01469
01470 generate_image (nx, ny, nz, power, ax, ay, az, zsep, fim);
01471
01472
01473
01474 if (filter) gaussian_filter (nx, ny, nz, dx, dy, dz, rmm,
01475 sigmax, sigmay, sigmaz, fim);
01476
01477
01478
01479 if (egfw) estimate_gfw (nx, ny, nz, dx, dy, dz,
01480 niter, quiet, fim, &avgsx, &avgsy, &avgsz);
01481
01482
01483
01484
01485 if (power) get_activation_region (nx, ny, nz, ax, ay, az, pthr, zsep,
01486 fim, arfim);
01487
01488
01489
01490 if (power) threshold_data (ax, ay, az,
01491 arfim, power_thr, &count, &sum, &sumsq,
01492 quiet, iter);
01493 else
01494 threshold_data (nx, ny, nz,
01495 fim, pthr, &count, &sum, &sumsq,
01496 quiet, iter);
01497
01498
01499
01500 if (mask_vol && (!power)) apply_mask (nx, ny, nz, fim);
01501
01502
01503
01504 if (power)
01505 identify_clusters (ax, ay, az, dx, dy, dz, rmm, arfim, quiet,
01506 freq_table, max_table);
01507 else
01508 identify_clusters (nx, ny, nz, dx, dy, dz, rmm, fim, quiet,
01509 freq_table, max_table);
01510
01511 }
01512
01513
01514
01515 output_results (nx, ny, nz, dx, dy, dz, filter, sigmax, sigmay, sigmaz,
01516 egfw, avgsx, avgsy, avgsz, power, ax, ay, az, zsep,
01517 rmm, pthr, niter, outfilename, freq_table, max_table);
01518
01519
01520
01521 terminate (&fim, &arfim, &freq_table, &max_table);
01522
01523 exit(0);
01524 }
01525
01526
01527
01528
01529
