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 }

00113 {
00114    fprintf (stderr, "%s Error: %s \n", PROGRAM_NAME, message);
00115    exit(1);
00116 }

01131 {
01132   int ixyz;
01133   int nxyz;
01134 
01135 
01136   /*----- initialize local variables -----*/
01137   nxyz = nx * ny * nz;
01138 
01139 
01140   /*----- apply mask to volume data -----*/
01141   for (ixyz = 0;  ixyz < nxyz;  ixyz++)
01142     if (! mask_vol[ixyz])
01143       fim[ixyz] = 0.0;
01144 
01145 
01146 }

00603 {
00604   FILE * fout = NULL;
00605   char message[MAX_NAME_LENGTH];     /* error message */
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       /*----- see if output file already exists -----*/
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 }

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 }

00876 {
00877   int nxy, nxyz;                /* total number of voxels */
00878   int ixyz;                     /* voxel index */
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   /*----- initialize local variables -----*/
00890   nxy = nx * ny;
00891   nxyz = nxy * nz;
00892 
00893 
00894   /*----- estimate the variance of the data -----*/
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   /*----- estimate the partial derivatives -----*/
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   /*----- estimate the variance of the partial derivatives -----*/
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   /*----- now estimate the equivalent Gaussian filter width -----*/
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   /*-----  save results  -----*/
00979   *avgsx += sx / niter;
00980   *avgsy += sy / niter;
00981   *avgsz += sz / niter;
00982 
00983   /*-----  output results  -----*/
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 }

00859 {
00860 
00861   /*----- use Gaussian blur routine -----*/ 
00862   EDIT_blur_volume_3d (nx, ny, nz, dx, dy, dz, 
00863                        MRI_float, fim, sigmax, sigmay, sigmaz);
00864 
00865 }

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   /*----- initialize local variables -----*/
00817   nxy = nx * ny;
00818   nxyz = nxy * nz;
00819   nxyzdiv2 = nxyz / 2;
00820 
00821   /*----- generate random image -----*/
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   /*----- if power calculation, generate "island" of activation -----*/
00832   if (power)
00833     {
00834       /*--- calculate dimensions of activation region ---*/
00835       activation_region (nx, ny, nz, ax, ay, az, 
00836                          &xbot, &xtop, &ybot, &ytop, &zbot, &ztop);
00837 
00838       /*--- add z-score offset to voxels within activation region ---*/
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 }

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   /*----- initialize local variables -----*/
01009   nxy = nx * ny;
01010   nxyz = nxy * nz;
01011   axy = ax * ay;
01012   axyz = axy * az;
01013 
01014 
01015   /*--- calculate dimensions of activation region ---*/
01016   activation_region (nx, ny, nz, ax, ay, az, 
01017                      &xbot, &xtop, &ybot, &ytop, &zbot, &ztop);
01018 
01019   /*--- copy activation region ---*/
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 }

00171 {
00172   int nopt = 1;                  /* input option argument counter */
00173   int ival;                      /* integer input */
00174   float fval;                    /* float input */
00175   char message[MAX_NAME_LENGTH];         /* error message */
00176   int mask_nx, mask_ny, mask_nz, mask_nvox;   /* mask dimensions */
00177   float  mask_dx, mask_dy, mask_dz;            
00178 
00179   
00180   /*----- does user request help menu? -----*/
00181   if (argc < 2 || strncmp(argv[1], "-help", 5) == 0)  display_help_menu();  
00182   
00183   
00184   /*----- add to program log -----*/
00185   AFNI_logger (PROGRAM_NAME,argc,argv); 
00186 
00187 
00188   /*----- initialize the input options -----*/
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   /*----- main loop over input options -----*/
00194   while (nopt < argc )
00195     {
00196       
00197       /*-----   -nx n  -----*/
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       /*-----   -ny n  -----*/
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       /*-----   -nz n  -----*/
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       /*-----   -dx d   -----*/
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       /*-----   -dy d   -----*/
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       /*-----   -dz d   -----*/
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       /**** -mask mset [14 June 2000] ****/
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       /*-----   -fwhm s   -----*/
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       /*-----   -fwhmx sx   -----*/
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       /*-----   -fwhmy sy   -----*/
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       /*-----   -fwhmz sz   -----*/
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       /*-----   -sigma s   -----*/
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       /*-----   -sigmax sx   -----*/
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       /*-----   -sigmay sy   -----*/
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       /*-----   -sigmaz sz   -----*/
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       /*-----   -egfw   -----*/
00441       if (strncmp(argv[nopt], "-egfw", 5) == 0)
00442         {
00443           *egfw = 1;
00444           nopt++;
00445           continue;
00446         }
00447 
00448 
00449       /*-----   -power   -----*/
00450       if (strncmp(argv[nopt], "-power", 6) == 0)
00451         {
00452           *power = 1;
00453           nopt++;
00454           continue;
00455         }
00456 
00457 
00458       /*-----   -ax n  -----*/
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       /*-----   -ay n  -----*/
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       /*-----   -az n  -----*/
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       /*-----   -zsep z   -----*/
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       /*-----   -rmm r   -----*/
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       /*-----   -pthr p   -----*/
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       /*-----   -iter n  -----*/
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       /*-----   -quiet q  -----*/
00557       if (strncmp(argv[nopt], "-quiet", 6) == 0)
00558         {
00559           *quiet = 1;
00560           nopt++;
00561           continue;
00562         }
00563 
00564 
00565       /*-----   -out filename   -----*/
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       /*----- unknown command -----*/
00578       AlphaSim_error ("unrecognized command line option ");
00579     }
00580 
00581 
00582   /*----- If mask dataset is used, set dimensions accordingly -----*/
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 }

01169 {
01170   MCW_cluster_array * clar;
01171   MCW_cluster * cl;
01172   int nxy;
01173   int nxyz;                     /* total number of voxels */
01174   int iclu, ipt;
01175   int size, max_size;
01176   int count;
01177 
01178 
01179   /*----- initialize local variables -----*/
01180   nxy = nx * ny;
01181 
01182   /*----- create array of clusters -----*/
01183   clar = MCW_find_clusters( nx,ny,nz , dx,dy,dz ,
01184                             MRI_float , fim , rmm ) ;
01185 
01186   /*----- record cluster sizes -----*/
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 }

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   /*----- get command line inputs -----*/
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   /*----- check for valid inputs -----*/
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   /*----- allocate memory space for image data -----*/
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   /*-- if power calculation, allocate memory space for activation region --*/
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   /*----- allocate memory space and initialize frequency table -----*/   
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   /*----- allocate memory space and initialize max cluster size table -----*/  
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   /*----- initialize voxel intensity sums -----*/
00713   *count = 0;
00714   *sum = 0.0;  
00715   *sumsq = 0.0;
00716        
00717 
00718   /*----- calculate power threshold -----*/
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   /*----- initialize ave. est. gaussian filter widths -----*/
00736   if (*egfw) 
00737     {
00738       *avgsx = 0.0;   *avgsy = 0.0;   *avgsz = 0.0;
00739     }
00740 
00741   /*----- initialize random number generator -----*/
00742   srand48 (1234567);
00743 
00744 }

00133 {
00134   *nx = 0;                   /* number of voxels along x-axis */
00135   *ny = 0;                   /* number of voxels along y-axis */
00136   *nz = 0;                   /* number of voxels along z-axis */
00137   *dx = 0.0;                 /* voxel size along x-axis */
00138   *dy = 0.0;                 /* voxel size along y-axis */
00139   *dz = 0.0;                 /* voxel size along z-axis */
00140   *filter = 0;               /* flag for Gaussian filtering */
00141   *sigmax = 0.0;             /* Gaussian filter width, x-axis (1 sigma) */
00142   *sigmay = 0.0;             /* Gaussian filter width, y-axis (1 sigma) */
00143   *sigmaz = 0.0;             /* Gaussian filter width, z-axis (1 sigma) */
00144   *egfw = 0;                 /* flag for estimation of filter width */
00145   *power = 0;                /* flag for power calculations */
00146   *ax = 0;                   /* number of activation voxels along x-axis */
00147   *ay = 0;                   /* number of activation voxels along y-axis */
00148   *az = 0;                   /* number of activation voxels along z-axis */
00149   *zsep = 0.0;               /* z-score signal and noise separation */
00150   *rmm = 0.0;                /* cluster connection radius (mm) */
00151   *pthr = 0.0;               /* individual voxel threshold prob. */
00152   *niter = 0;                /* number of Monte Carlo simulations  */
00153   *quiet = 0;                /* generate screen output (default)  */
00154   *outfilename = NULL;       /* name of output file */
00155 }

01410 {
01411   int nx;                  /* number of voxels along x-axis */
01412   int ny;                  /* number of voxels along y-axis */
01413   int nz;                  /* number of voxels along z-axis */
01414   float dx;                /* voxel size along x-axis */
01415   float dy;                /* voxel size along y-axis */
01416   float dz;                /* voxel size along z-axis */
01417   int filter;              /* flag for Gaussian filtering */
01418   float sigmax;            /* Gaussian filter width, x-axis (1 sigma) */
01419   float sigmay;            /* Gaussian filter width, y-axis (1 sigma) */
01420   float sigmaz;            /* Gaussian filter width, z-axis (1 sigma) */
01421   int egfw;                /* flag for estimation of filter width */
01422   float avgsx;             /* est. Gaussian filter width, x-axis (1 sigma) */
01423   float avgsy;             /* est. Gaussian filter width, x-axis (1 sigma) */
01424   float avgsz;             /* est. Gaussian filter width, x-axis (1 sigma) */
01425   int power;               /* flag for perform power calculations */
01426   int ax;                  /* number of activation voxels along x-axis */
01427   int ay;                  /* number of activation voxels along y-axis */
01428   int az;                  /* number of activation voxels along z-axis */
01429   float zsep;              /* z-score separation between signal and noise */
01430   float rmm;               /* cluster connection radius (mm) */
01431   float pthr;              /* individual voxel threshold probability */
01432   int niter;               /* number of Monte Carlo simulations */
01433   int quiet;               /* set to 1 to suppress screen output */
01434   char * outfilename;      /* name of output file */
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   /*----- Identify software -----*/
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   /*----- program initialization -----*/
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   /*----- Monte Carlo iteration -----*/
01465   for (iter = 1;  iter <= niter;  iter++)
01466     {
01467       if (!quiet)  printf ("Iter =%5d  \n", iter);
01468 
01469       /*----- generate volume of random voxel intensities -----*/
01470       generate_image (nx, ny, nz, power, ax, ay, az, zsep, fim);
01471 
01472       
01473       /*----- apply gaussian filter to volume data -----*/
01474       if (filter)  gaussian_filter (nx, ny, nz, dx, dy, dz, rmm,
01475                                     sigmax, sigmay, sigmaz, fim);
01476 
01477 
01478       /*----- estimate equivalent gaussian filter width -----*/
01479       if (egfw)  estimate_gfw (nx, ny, nz, dx, dy, dz, 
01480                                niter, quiet, fim, &avgsx, &avgsy, &avgsz);
01481 
01482 
01483       /*----- if power calculation, get volume corresponding to   -----*/
01484       /*----- activation region and corresponding power threshold -----*/
01485       if (power)  get_activation_region (nx, ny, nz, ax, ay, az, pthr, zsep, 
01486                                          fim, arfim);
01487 
01488 
01489       /*----- apply threshold to volume data -----*/
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       /*----- apply mask to volume data -----*/
01500       if (mask_vol && (!power))  apply_mask (nx, ny, nz, fim);
01501 
01502 
01503       /*----- identify clusters -----*/
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   /*----- generate requested output -----*/
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   /*----- terminate program -----*/
01521   terminate (&fim, &arfim, &freq_table, &max_table);
01522 
01523   exit(0);
01524 }

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 }

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];     /* error message */
01247   FILE * fout;
01248 
01249   
01250   /*----- allocate memory space for probability table -----*/   
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   /*----- allocate memory space for alpha table -----*/   
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   /*----- allocate memory space for cum. prop. of cluster size table  -----*/ 
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   /*----- if output file has not been specified, use stdout -----*/
01300   if (outfilename == NULL)
01301     fout = stdout;
01302   else
01303     {
01304       /*----- see if output file already exists -----*/
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       /*----- open file for output -----*/
01313       fout = fopen (outfilename, "w");
01314       if (fout == NULL)
01315         { 
01316           AlphaSim_error ("unable to write file ");
01317         }
01318     }
01319 
01320   /*----- print out the results -----*/
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 }

01041 {
01042   int nxyz;
01043   int ixyz;
01044   int pcount;
01045   float p;
01046   
01047   
01048   /*----- initialize local variables -----*/
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 }

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 }

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   /*----- initialize local variables -----*/
01082   nxyz = nx * ny * nz;
01083 
01084 
01085   /*----- update sums -----*/
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   /*----- calculate z-threshold -----*/
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   /*----- apply threshold to image data -----*/
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 }

00753 {
00754   return ( (float)drand48() );
00755 }