00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2001, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006 
00007 /*---------------------------------------------------------------------------*/
00008 /*
00009   Program to correct for image intensity non-uniformity.
00010 
00011   File:    3dUniformize.c
00012   Author:  B. Douglas Ward
00013   Date:    28 January 2000
00014 
00015   Mod:     Added call to AFNI_logger.
00016   Date:    15 August 2001
00017 
00018 */
00019 
00020 /*---------------------------------------------------------------------------*/
00021 
00022 #define PROGRAM_NAME "3dUniformize"                  /* name of this program */
00023 #define PROGRAM_AUTHOR "B. D. Ward"                        /* program author */
00024 #define PROGRAM_INITIAL "28 January 2000" /* date of initial program release */
00025 #define PROGRAM_LATEST  "16 April 2003"   /* date of latest program revision */
00026 
00027 /*---------------------------------------------------------------------------*/
00028 /*
00029   Include header files.
00030 */
00031 
00032 #include "mrilib.h"
00033 #include "matrix.h"
00034 
00035 
00036 /*---------------------------------------------------------------------------*/
00037 /*
00038   Global variables, constants, and data structures.
00039 */
00040 
00041 #define MAX_STRING_LENGTH 80
00042 
00043 static THD_3dim_dataset * anat_dset = NULL;     /* input anatomical dataset  */
00044 char * commandline = NULL ;                /* command line for history notes */
00045 
00046 int input_datum = MRI_short ;              /* 16 Apr 2003 - RWCox */
00047 int quiet       = 0 ;                      /* ditto */
00048 #define USE_QUIET
00049  
00050 typedef struct UN_options
00051 { 
00052   char * anat_filename;       /* file name for input anat dataset */
00053   char * prefix_filename;     /* prefix name for output dataset */
00054   Boolean quiet;              /* flag for suppress screen output */
00055   int lower_limit;    /* lower limit for voxel intensity */
00056   int rpts;           /* #voxels in sub-sampled image (for pdf) */
00057   int spts;           /* #voxels in subsub-sampled image (for field poly.) */
00058   int nbin;           /* #bins for pdf estimation */
00059   int npar;           /* #parameters for field polynomial */
00060 } UN_options;
00061 
00062 
00063 /*---------------------------------------------------------------------------*/
00064 /*
00065   Include source code files.
00066 */
00067 
00068 #include "matrix.c"
00069 #include "estpdf3.c"
00070 
00071 
00072 /*---------------------------------------------------------------------------*/
00073 /*
00074    Print error message and stop.
00075 */
00076 
00077 void UN_error (char * message)
00078 {
00079   fprintf (stderr, "%s Error: %s \n", PROGRAM_NAME, message);
00080   exit(1);
00081 }
00082 
00083 
00084 /*---------------------------------------------------------------------------*/
00085 
00086 /** macro to test a malloc-ed pointer for validity **/
00087 
00088 #define MTEST(ptr) \
00089 if((ptr)==NULL) \
00090 ( UN_error ("Cannot allocate memory") )
00091      
00092 
00093 /*---------------------------------------------------------------------------*/
00094 /*
00095    Routine to display 3dUniformize help menu.
00096 */
00097 
00098 void display_help_menu()
00099 {
00100   printf 
00101     (
00102      "This program corrects for image intensity non-uniformity.\n\n"
00103      "Usage: \n"
00104      "3dUniformize  \n"
00105      "-anat filename    Filename of anat dataset to be corrected            \n"
00106      "                                                                      \n"
00107      "[-quiet]          Suppress output to screen                           \n"
00108      "                                                                      \n"
00109      "-prefix pname     Prefix name for file to contain corrected image     \n"
00110       );
00111 
00112   printf("\n" MASTER_SHORTHELP_STRING ) ;
00113   exit(0);
00114 }
00115 
00116 
00117 /*---------------------------------------------------------------------------*/
00118 /*
00119   Routine to initialize the input options.
00120 */
00121  
00122 void initialize_options 
00123 (
00124   UN_options * option_data    /* uniformization program options */
00125 )
00126  
00127 {
00128 
00129   /*----- initialize default values -----*/
00130   option_data->anat_filename = NULL;    /* file name for input anat dataset */
00131   option_data->prefix_filename = NULL;  /* prefix name for output dataset */
00132   option_data->quiet = FALSE;           /* flag for suppress screen output */
00133   option_data->lower_limit = 25;        /* voxel intensity lower limit */
00134   option_data->rpts = 200000;   /* #voxels in sub-sampled image (for pdf) */
00135   option_data->spts = 10000;    /* #voxels in subsub-sampled image 
00136                                    (for field polynomial estimation) */
00137   option_data->nbin = 250;      /* #bins for pdf estimation */
00138   option_data->npar = 35;       /* #parameters for field polynomial */
00139 
00140 }
00141 
00142 
00143 /*---------------------------------------------------------------------------*/
00144 /*
00145   Routine to get user specified input options.
00146 */
00147 
00148 void get_options
00149 (
00150   int argc,                        /* number of input arguments */
00151   char ** argv,                    /* array of input arguments */ 
00152   UN_options * option_data         /* uniformization program options */
00153 )
00154 
00155 {
00156   int nopt = 1;                     /* input option argument counter */
00157   int ival, index;                  /* integer input */
00158   float fval;                       /* float input */
00159   char message[MAX_STRING_LENGTH];  /* error message */
00160 
00161 
00162   /*----- does user request help menu? -----*/
00163   if (argc < 2 || strncmp(argv[1], "-help", 5) == 0)  display_help_menu();  
00164    
00165   
00166   /*----- add to program log -----*/
00167   AFNI_logger (PROGRAM_NAME,argc,argv); 
00168 
00169 
00170   /*----- main loop over input options -----*/
00171   while (nopt < argc )
00172     {
00173 
00174       /*-----   -anat filename   -----*/
00175       if (strncmp(argv[nopt], "-anat", 5) == 0)
00176         {
00177           nopt++;
00178           if (nopt >= argc)  UN_error ("need argument after -anat ");
00179           option_data->anat_filename = 
00180             malloc (sizeof(char) * MAX_STRING_LENGTH);
00181           MTEST (option_data->anat_filename);
00182           strcpy (option_data->anat_filename, argv[nopt]);
00183 
00184           anat_dset = THD_open_dataset (option_data->anat_filename);
00185           if (!ISVALID_3DIM_DATASET (anat_dset))
00186             {
00187               sprintf (message, "Can't open dataset: %s\n", 
00188                        option_data->anat_filename); 
00189               UN_error (message); 
00190             } 
00191           THD_load_datablock (anat_dset->dblk); 
00192           if (DSET_ARRAY(anat_dset,0) == NULL)
00193             {
00194               sprintf (message, "Can't access data: %s\n", 
00195                        option_data->anat_filename); 
00196               UN_error (message); 
00197             }
00198 
00199           /** RWCox [16 Apr 2003]
00200               If input is a byte dataset, make a short copy of it. **/
00201 
00202           if( DSET_BRICK_TYPE(anat_dset,0) == MRI_byte ){
00203 
00204             THD_3dim_dataset *qset ;
00205             register byte *bar ; register short *sar ;
00206             register int ii,nvox ;
00207 
00208             fprintf(stderr,"++ WARNING: converting input dataset from byte to short\n") ;
00209             qset = EDIT_empty_copy(anat_dset) ;
00210             nvox = DSET_NVOX(anat_dset) ;
00211             bar  = (byte *) DSET_ARRAY(anat_dset,0) ;
00212             sar  = (short *)malloc(sizeof(short)*nvox) ;
00213             for( ii=0 ; ii < nvox ; ii++ ) sar[ii] = (short) bar[ii] ;
00214             EDIT_substitute_brick( qset , 0 , MRI_short , sar ) ;
00215             DSET_delete(anat_dset) ; anat_dset = qset ; input_datum = MRI_byte ;
00216 
00217           } else if ( DSET_BRICK_TYPE(anat_dset,0) != MRI_short ){
00218 
00219             fprintf(stderr,"** ERROR: input dataset not short or byte type!\n") ;
00220             exit(1) ;
00221 
00222           }
00223 
00224           nopt++;
00225           continue;
00226         }
00227       
00228 
00229       /*-----   -quiet  -----*/
00230       if (strncmp(argv[nopt], "-quiet", 6) == 0)
00231         {
00232           option_data->quiet = TRUE;
00233           quiet = 1 ;                /* 16 Apr 2003 */
00234           nopt++;
00235           continue;
00236         }
00237 
00238 
00239       /*-----   -prefix prefixname   -----*/
00240       if (strncmp(argv[nopt], "-prefix", 7) == 0)
00241         {
00242           nopt++;
00243           if (nopt >= argc)  UN_error ("need argument after -prefix ");
00244           option_data->prefix_filename = 
00245             malloc (sizeof(char) * MAX_STRING_LENGTH);
00246           MTEST (option_data->prefix_filename);
00247           strcpy (option_data->prefix_filename, argv[nopt]);
00248           nopt++;
00249           continue;
00250         }
00251       
00252 
00253       /*----- unknown command -----*/
00254       sprintf(message,"Unrecognized command line option: %s\n", argv[nopt]);
00255       UN_error (message);
00256       
00257     }
00258 
00259   
00260 }
00261 
00262 
00263 /*---------------------------------------------------------------------------*/
00264 /*
00265   Routine to check whether one output file already exists.
00266 */
00267 
00268 void check_one_output_file 
00269 (
00270   char * filename                   /* name of output file */
00271 )
00272 
00273 {
00274   char message[MAX_STRING_LENGTH];    /* error message */
00275   THD_3dim_dataset * new_dset=NULL;   /* output afni data set pointer */
00276   int ierror;                         /* number of errors in editing data */
00277 
00278   
00279   /*----- make an empty copy of input dataset -----*/
00280   new_dset = EDIT_empty_copy ( anat_dset );
00281   
00282   
00283   ierror = EDIT_dset_items( new_dset ,
00284                             ADN_prefix , filename ,
00285                             ADN_label1 , filename ,
00286                             ADN_self_name , filename ,
00287                             ADN_none ) ;
00288   
00289   if( ierror > 0 )
00290     {
00291       sprintf (message,
00292                "*** %d errors in attempting to create output dataset!\n", 
00293                ierror);
00294       UN_error (message);
00295     }
00296   
00297   if( THD_is_file(new_dset->dblk->diskptr->header_name) )
00298     {
00299       sprintf (message,
00300                "Output dataset file %s already exists"
00301                "--cannot continue!\a\n",
00302                new_dset->dblk->diskptr->header_name);
00303       UN_error (message);
00304     }
00305   
00306   /*----- deallocate memory -----*/   
00307   THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
00308   
00309 }
00310 
00311 
00312 /*---------------------------------------------------------------------------*/
00313 
00314 void verify_inputs 
00315 (  
00316   UN_options * option_data         /* uniformization program options */
00317 )
00318 
00319 {
00320 }
00321 
00322 
00323 /*---------------------------------------------------------------------------*/
00324 /*
00325   Program initialization.
00326 */
00327 
00328 void initialize_program 
00329 (
00330   int argc,                        /* number of input arguments */
00331   char ** argv,                    /* array of input arguments */ 
00332   UN_options ** option_data,       /* uniformization program options */
00333   short ** sfim                    /* output image volume */
00334 )
00335 
00336 {
00337   int nxyz;                        /* #voxels in input dataset */
00338 
00339 
00340   /*----- Save command line for history notes -----*/
00341   commandline = tross_commandline( PROGRAM_NAME , argc,argv ) ;
00342 
00343 
00344   /*----- Allocate memory for input options -----*/
00345   *option_data = (UN_options *) malloc (sizeof(UN_options));
00346   MTEST (*option_data);
00347 
00348   
00349   /*----- Initialize the input options -----*/
00350   initialize_options (*option_data); 
00351 
00352 
00353   /*----- Get operator inputs -----*/
00354   get_options (argc, argv, *option_data);
00355 
00356 
00357   /*----- Verify that inputs are acceptable -----*/
00358   verify_inputs (*option_data);
00359 
00360 
00361   /*----- Initialize random number generator -----*/
00362   rand_initialize (1234567);
00363 
00364 
00365   /*----- Allocate memory for output volume -----*/
00366   nxyz = DSET_NX(anat_dset) * DSET_NY(anat_dset) * DSET_NZ(anat_dset);
00367   *sfim = (short *) malloc (sizeof(short) * nxyz);
00368   MTEST (*sfim);
00369 }
00370 
00371 
00372 /*---------------------------------------------------------------------------*/
00373 /*
00374   Write time series data to specified file.
00375 */
00376 
00377 void ts_write (char * filename, int ts_length, float * data)
00378 {
00379   int it;
00380   FILE * outfile = NULL;
00381 
00382   outfile = fopen (filename, "w");
00383 
00384 
00385   for (it = 0;  it < ts_length;  it++)
00386     {
00387       fprintf (outfile, "%f ", data[it]);
00388       fprintf (outfile, " \n");
00389     }
00390 
00391   fclose (outfile);
00392 }
00393 
00394 
00395 /*---------------------------------------------------------------------------*/
00396 /*
00397   Resample the original image at randomly selected voxels (whose intensity
00398   value is greater than the specified lower limit, to exclude voxels outside
00399   the brain).  Take the logarithm of the intensity values for the selected 
00400   voxels.
00401 */
00402 
00403 void resample 
00404 (
00405   UN_options * option_data,
00406   int * ir,                         /* voxel indices for resampled image */
00407   float * vr                        /* resampled image data (logarithms) */
00408 )
00409 
00410 {
00411   short * anat_data = NULL;
00412   int nxyz;
00413   int rpts;
00414   int lower_limit;
00415   int it, k;
00416 
00417 
00418   /*----- Initialize local variables -----*/
00419   nxyz = DSET_NX(anat_dset) * DSET_NY(anat_dset) * DSET_NZ(anat_dset);
00420   anat_data = (short *) DSET_BRICK_ARRAY(anat_dset,0);
00421   lower_limit = option_data->lower_limit;
00422   rpts = option_data->rpts;
00423 
00424   it = 0;
00425   while (it < rpts)
00426     {
00427       k = rand_uniform (0, nxyz);
00428       if ( (k >= 0) && (k < nxyz) && (anat_data[k] > lower_limit) )
00429         {
00430           ir[it] = k;
00431           vr[it] = log (anat_data[k] + rand_uniform (0.0,1.0));
00432           it++;
00433         }
00434     }
00435 
00436   return;
00437 }
00438 
00439 
00440 /*---------------------------------------------------------------------------*/
00441 /*
00442   Create intensity map that will tend to concentrate values around the
00443   means of the gray and white matter distributions.
00444 */
00445 
00446 void create_map (pdf vpdf, float * pars, float * vtou)
00447 
00448 {
00449   int ibin;
00450   float v;
00451 
00452   for (ibin = 0;  ibin < vpdf.nbin;  ibin++)
00453     {
00454       v = PDF_ibin_to_xvalue (vpdf, ibin);
00455           
00456       if ((v > pars[4]-2.0*pars[5]) && (v < 0.5*(pars[4]+pars[7])))
00457         vtou[ibin] = pars[4];
00458       else if ((v > 0.5*(pars[4]+pars[7])) && (v < pars[7]+2.0*pars[8]))
00459         vtou[ibin] = pars[7];
00460       else
00461         vtou[ibin] = v;
00462 
00463     }
00464   
00465 }
00466 
00467 
00468 /*---------------------------------------------------------------------------*/
00469 /*
00470   Use the intensity map to transform values of voxel intensities.
00471 */
00472 
00473 void map_vtou (pdf vpdf, int rpts, float * vr, float * vtou, float * ur)
00474 
00475 {
00476   int i, ibin;
00477   float v;
00478 
00479 
00480   for (i = 0;  i < rpts;  i++)
00481     {
00482       v = vr[i];
00483       ibin = PDF_xvalue_to_ibin (vpdf, v);
00484       
00485       if ((ibin >= 0) && (ibin < vpdf.nbin))
00486         ur[i] = vtou[ibin];
00487       else
00488         ur[i] = v;
00489     }
00490 
00491 }
00492 
00493 
00494 /*---------------------------------------------------------------------------*/
00495 
00496 void subtract (int rpts, float * a, float * b, float * c)
00497 
00498 {
00499   int i;
00500 
00501 
00502   for (i = 0;  i < rpts;  i++)
00503     {
00504       c[i] = a[i] - b[i];
00505     }
00506 
00507 }
00508 
00509 
00510 /*---------------------------------------------------------------------------*/
00511 /*
00512   Create one row of the X matrix.
00513 */
00514 
00515 void create_row (int ixyz, int nx, int ny, int nz, float * xrow)
00516 
00517 {
00518   int ix, jy, kz;
00519   float x, y, z, x2, y2, z2, x3, y3, z3, x4, y4, z4;
00520 
00521 
00522   IJK_TO_THREE (ixyz, ix, jy, kz, nx, nx*ny); 
00523 
00524 
00525   x = (float) ix / (float) nx - 0.5;
00526   y = (float) jy / (float) ny - 0.5;
00527   z = (float) kz / (float) nz - 0.5;
00528 
00529   x2 = x*x;   x3 = x*x2;   x4 = x2*x2;
00530   y2 = y*y;   y3 = y*y2;   y4 = y2*y2;
00531   z2 = z*z;   z3 = z*z2;   z4 = z2*z2;
00532 
00533 
00534   xrow[0] = 1.0;
00535   xrow[1] = x;
00536   xrow[2] = y;
00537   xrow[3] = z;
00538   xrow[4] = x*y;
00539   xrow[5] = x*z;
00540   xrow[6] = y*z;
00541   xrow[7] = x2;
00542   xrow[8] = y2;
00543   xrow[9] = z2;
00544   xrow[10] = x*y*z;
00545   xrow[11] = x2*y;
00546   xrow[12] = x2*z;
00547   xrow[13] = y2*x;
00548   xrow[14] = y2*z;
00549   xrow[15] = z2*x;
00550   xrow[16] = z2*y;
00551   xrow[17] = x3;
00552   xrow[18] = y3;
00553   xrow[19] = z3;
00554   xrow[20] = x2*y*z;
00555   xrow[21] = x*y2*z;
00556   xrow[22] = x*y*z2;
00557   xrow[23] = x2*y2;
00558   xrow[24] = x2*z2;
00559   xrow[25] = y2*z2;
00560   xrow[26] = x3*y;
00561   xrow[27] = x3*z;
00562   xrow[28] = x*y3;
00563   xrow[29] = y3*z;
00564   xrow[30] = x*z3;
00565   xrow[31] = y*z3;
00566   xrow[32] = x4;
00567   xrow[33] = y4;
00568   xrow[34] = z4;
00569 
00570 
00571   return;
00572 }
00573 
00574 
00575 /*---------------------------------------------------------------------------*/
00576 /*
00577   Approximate the distortion field with a polynomial function in 3 dimensions.
00578 */
00579 
00580 void poly_field (int nx, int ny, int nz, int rpts, int * ir, float * fr, 
00581                  int spts, int npar, float * fpar)
00582 
00583 {
00584   int p;                   /* number of parameters in the full model */ 
00585   int i, j, k;
00586   matrix x;                      /* independent variable matrix */
00587   matrix xtxinv;                 /* matrix:  1/(X'X)       */
00588   matrix xtxinvxt;               /* matrix:  (1/(X'X))X'   */
00589   vector y;
00590   vector coef;
00591   float * xrow = NULL;
00592   int ip;
00593   int iter;
00594   float f;
00595 
00596 
00597   p = npar;
00598 
00599 
00600   /*----- Initialize matrices and vectors -----*/
00601   matrix_initialize (&x);
00602   matrix_initialize (&xtxinv);
00603   matrix_initialize (&xtxinvxt);
00604   vector_initialize (&y);
00605   vector_initialize (&coef);
00606 
00607 
00608   /*----- Allocate memory -----*/
00609   matrix_create (spts, p, &x);
00610   vector_create (spts, &y);
00611   xrow = (float *) malloc (sizeof(float) * p); 
00612  
00613 
00614   /*----- Set up the X matrix and Y vector -----*/
00615   for (i = 0;  i < spts;  i++)
00616     {
00617       k = rand_uniform (0, rpts);
00618       create_row (ir[k], nx, ny, nz, xrow);
00619 
00620       for (j = 0;  j < p;  j++)
00621         x.elts[i][j] = xrow[j];
00622       y.elts[i] = fr[k];
00623     }
00624 
00625 
00626   /*  
00627       matrix_sprint ("X matrix = ", x);
00628       vector_sprint ("Y vector = ", y);
00629   */
00630 
00631 
00632   {
00633     /*----- calculate various matrices which will be needed later -----*/
00634     matrix xt, xtx;            /* temporary matrix calculation results */
00635     int ok;                    /* flag for successful matrix inversion */
00636     
00637     
00638     /*----- initialize matrices -----*/
00639     matrix_initialize (&xt);
00640     matrix_initialize (&xtx);
00641     
00642         
00643     matrix_transpose (x, &xt);
00644     matrix_multiply (xt, x, &xtx);
00645     ok = matrix_inverse (xtx, &xtxinv);
00646     
00647     if (ok)
00648       matrix_multiply (xtxinv, xt, &xtxinvxt);
00649     else
00650       {
00651         matrix_sprint ("X matrix = ", x);
00652         matrix_sprint ("X'X matrix = ", xtx);
00653         UN_error ("Improper X matrix  (cannot invert X'X) ");
00654       }
00655     
00656     /*----- dispose of matrices -----*/
00657     matrix_destroy (&xtx);
00658     matrix_destroy (&xt);
00659     
00660   }
00661 
00662 
00663   /*
00664     matrix_sprint ("1/(X'X)     = ", xtxinv);
00665     matrix_sprint ("(1/(X'X))X' = ", xtxinvxt);
00666     vector_sprint ("Y data  = ", y);
00667   */
00668   
00669   vector_multiply (xtxinvxt, y, &coef);
00670   /*
00671     vector_sprint ("Coef    = ", coef);
00672   */
00673   
00674 
00675   for (ip = 0;  ip < p;  ip++)
00676     {
00677      fpar[ip] = coef.elts[ip];
00678     }
00679   
00680 
00681   /*----- Dispose of matrices and vectors -----*/
00682   matrix_destroy (&x);
00683   matrix_destroy (&xtxinv);
00684   matrix_destroy (&xtxinvxt);
00685   vector_destroy (&y);
00686   vector_destroy (&coef);
00687 
00688   
00689 }
00690 
00691 
00692 /*---------------------------------------------------------------------------*/
00693 /*
00694   Use the 3-dimensional polynomial function to estimate the distortion field
00695   at each point.
00696 */
00697 
00698 float warp_image (int npar, float * fpar, int nx, int ny, int nz,
00699                   int rpts, int * ir, float * fs)
00700 {
00701   int i, j;
00702   float x;
00703   float * xrow;
00704   float max_warp;
00705 
00706 
00707   xrow = (float *) malloc (sizeof(float) * npar); 
00708 
00709 
00710   max_warp = 0.0;
00711 
00712   for (i = 0;  i < rpts;  i++)
00713     {
00714       create_row (ir[i], nx, ny, nz, xrow);
00715 
00716       fs[i] = 0.0;
00717             
00718       for (j = 1;  j < npar;  j++)
00719         fs[i] += fpar[j] * xrow[j];
00720 
00721       if (fabs(fs[i]) > max_warp)
00722         max_warp = fabs(fs[i]);
00723     }
00724 
00725 
00726   free (xrow);   xrow = NULL;
00727 
00728 
00729   return (max_warp);
00730 }
00731 
00732 
00733 /*---------------------------------------------------------------------------*/
00734 /*
00735   Find polynomial approximation to the distortion field.
00736 */
00737 
00738 void estimate_field (UN_options * option_data, 
00739                      int * ir, float * vr, float * fpar)
00740 {
00741   float * ur = NULL, * us = NULL, * fr = NULL, * fs = NULL, * wr = NULL;
00742   float * vtou = NULL;
00743   float * gpar;
00744   int iter = 0;
00745   int ip;
00746   int it;
00747   int nx, ny, nz, nxy, nxyz;
00748   int rpts, spts, nbin, npar;
00749   float parameters [DIMENSION];    /* parameters for PDF estimation */
00750   Boolean ok = TRUE;               /* flag for successful PDF estimation */
00751   char filename[MAX_STRING_LENGTH];
00752 
00753 
00754   /*----- Initialize local variables -----*/
00755   nx = DSET_NX(anat_dset);  ny = DSET_NY(anat_dset);  nz = DSET_NZ(anat_dset);
00756   nxy = nx*ny;   nxyz = nxy*nz;
00757   rpts = option_data->rpts;
00758   spts = option_data->spts;
00759   nbin = option_data->nbin;
00760   npar = option_data->npar;
00761   
00762 
00763 
00764   /*----- Allocate memory -----*/
00765   ur   = (float *) malloc (sizeof(float) * rpts);   MTEST (ur);
00766   us   = (float *) malloc (sizeof(float) * rpts);   MTEST (us);
00767   fr   = (float *) malloc (sizeof(float) * rpts);   MTEST (fr);
00768   fs   = (float *) malloc (sizeof(float) * rpts);   MTEST (fs);
00769   wr   = (float *) malloc (sizeof(float) * rpts);   MTEST (wr);
00770   gpar = (float *) malloc (sizeof(float) * npar);   MTEST (gpar);
00771   vtou = (float *) malloc (sizeof(float) * nbin);   MTEST (vtou);
00772 
00773 
00774   /*----- Initialize polynomial coefficients -----*/
00775   for (ip = 0;  ip < npar;  ip++)
00776     {
00777       fpar[ip] = 0.0;
00778       gpar[ip] = 0.0;
00779     }
00780 
00781 
00782   /*----- Estimate pdf for resampled data -----*/
00783   PDF_initialize (&p);
00784   PDF_float_to_pdf (rpts, vr, nbin, &p);
00785 
00786   if( !quiet ){
00787    sprintf (filename, "p%d.1D", iter);
00788    PDF_write_file (filename, p);
00789   }
00790 
00791 
00792   /*----- Estimate gross field distortion -----*/
00793   poly_field (nx, ny, nz, rpts, ir, vr, spts, npar, fpar);
00794   warp_image (npar, fpar, nx, ny, nz, rpts, ir, fs);
00795   subtract (rpts, vr, fs, ur);
00796  
00797   
00798   for (ip = 0;  ip < rpts;  ip++)
00799     vr[ip] = ur[ip];
00800 
00801 
00802   /*----- Iterate over field distortion for concentrating the PDF -----*/
00803   for (iter = 1;  iter <= 5;  iter++)
00804     {
00805       /*----- Estimate pdf for perturbed image ur -----*/
00806       estpdf_float (rpts, ur, nbin, parameters);
00807       PDF_sprint ("p", p);
00808       if( !quiet ){
00809        sprintf (filename, "p%d.1D", iter);
00810        PDF_write_file (filename, p);
00811       }
00812 
00813       /*----- Sharpen the pdf and produce modified image wr -----*/
00814       create_map (p, parameters, vtou);
00815       if( !quiet ){
00816        sprintf (filename, "vtou%d.1D", iter);
00817        ts_write (filename, p.nbin, vtou);
00818       }
00819       map_vtou (p, rpts, ur, vtou, wr);
00820 
00821       /*----- Estimate smooth distortion field fs -----*/
00822       subtract (rpts, vr, wr, fr);
00823       poly_field (nx, ny, nz, rpts, ir, fr, spts, npar, gpar);
00824       warp_image (npar, gpar, nx, ny, nz, rpts, ir, fs);
00825 
00826       /*----- Create perturbed image ur -----*/
00827       subtract (rpts, vr, fs, ur);
00828     }
00829   
00830 
00831   /*----- Accumulate distortion field polynomial coefficients -----*/
00832   for (ip = 0;  ip < npar;  ip++)
00833     fpar[ip] += gpar[ip];
00834 
00835 
00836   /*----- Deallocate memory -----*/
00837   free (ur);     ur = NULL;  
00838   free (us);     us = NULL;
00839   free (fr);     fr = NULL;
00840   free (fs);     fs = NULL;
00841   free (wr);     wr = NULL;
00842   free (gpar);   gpar = NULL;
00843   free (vtou);   vtou = NULL;
00844 
00845 
00846   return;
00847 }
00848 
00849 
00850 /*---------------------------------------------------------------------------*/
00851 /*
00852   Remove the nonuniformity field.
00853 */
00854 
00855 void remove_field (UN_options * option_data, float * fpar, short * sfim)
00856 {
00857   short * anat_data = NULL;
00858   int rpts;
00859   int npar;
00860   int lower_limit;
00861   int nx, ny, nz, nxyz;
00862   int ixyz, jpar;
00863   float x;
00864   float * xrow;
00865   float f;
00866 
00867 
00868   /*----- Initialize local variables -----*/
00869   nx = DSET_NX(anat_dset);  ny = DSET_NY(anat_dset);  nz = DSET_NZ(anat_dset);
00870   nxyz = nx*ny*nz;
00871   anat_data = (short *) DSET_BRICK_ARRAY(anat_dset,0);
00872   rpts = option_data->rpts;
00873   npar = option_data->npar;
00874   lower_limit = option_data->lower_limit;
00875 
00876   xrow = (float *) malloc (sizeof(float) * npar); 
00877 
00878 
00879   for (ixyz = 0;  ixyz < nxyz;  ixyz++)
00880     {
00881       if (anat_data[ixyz] > lower_limit) 
00882         {
00883           create_row (ixyz, nx, ny, nz, xrow);
00884           
00885           f = 0.0;
00886           for (jpar = 1;  jpar < npar;  jpar++)
00887             f += fpar[jpar] * xrow[jpar];
00888           
00889           sfim[ixyz] = exp( log(anat_data[ixyz]) - f);
00890         }
00891       else
00892         sfim[ixyz] = anat_data[ixyz];
00893     }
00894 
00895   
00896   return;
00897 }
00898 
00899 
00900 /*---------------------------------------------------------------------------*/
00901 /*
00902   Correct for image intensity nonuniformity.
00903 */
00904 
00905 void uniformize (UN_options * option_data, short * sfim)
00906 
00907 {
00908   int * ir = NULL;
00909   float * vr = NULL;
00910   float * fpar = NULL;
00911   int rpts, npar;
00912 
00913 
00914   /*----- Initialize local variables -----*/
00915   rpts = option_data->rpts;
00916   npar = option_data->npar;
00917 
00918 
00919   /*----- Allocate memory -----*/
00920   ir = (int *) malloc (sizeof(int) * rpts);         MTEST(ir);
00921   vr = (float *) malloc (sizeof(float) * rpts);     MTEST(vr);
00922   fpar = (float *) malloc (sizeof(float) * npar);   MTEST(fpar);
00923 
00924 
00925   /*----- Resample the data -----*/
00926   resample (option_data, ir, vr);
00927 
00928 
00929   /*----- Estimate the nonuniformity field -----*/
00930   estimate_field (option_data, ir, vr, fpar);
00931 
00932 
00933   /*----- Remove the nonuniformity field -----*/
00934   remove_field (option_data, fpar, sfim);
00935 
00936  
00937   /*----- Deallocate memory -----*/
00938   free (ir);     ir = NULL;
00939   free (vr);     vr = NULL;
00940   free (fpar);   fpar = NULL;
00941 
00942 }
00943 
00944 
00945 /*---------------------------------------------------------------------------*/
00946 /*
00947   Routine to write one AFNI dataset.
00948   
00949   
00950 */
00951 
00952 void write_afni_data 
00953 (
00954   UN_options * option_data,
00955   short * sfim
00956 )
00957 
00958 {
00959   int nxyz;                           /* number of voxels */
00960   int ii;                             /* voxel index */
00961   THD_3dim_dataset * new_dset=NULL;   /* output afni data set pointer */
00962   int ierror;                         /* number of errors in editing data */
00963   int ibuf[32];                       /* integer buffer */
00964   float fbuf[MAX_STAT_AUX];           /* float buffer */
00965   float fimfac;                       /* scale factor for short data */
00966   int output_datum;                   /* data type for output data */
00967   char * filename;                    /* prefix filename for output */
00968   byte *bfim = NULL ;                 /* 16 Apr 2003 */
00969 
00970 
00971   /*----- initialize local variables -----*/
00972   filename = option_data->prefix_filename;
00973   nxyz = DSET_NX(anat_dset) * DSET_NY(anat_dset) * DSET_NZ(anat_dset);
00974 
00975   
00976   /*-- make an empty copy of this dataset, for eventual output --*/
00977   new_dset = EDIT_empty_copy( anat_dset ) ;
00978 
00979 
00980   /*----- Record history of dataset -----*/
00981   tross_Copy_History( anat_dset , new_dset ) ;
00982   if( commandline != NULL )
00983      tross_Append_History( new_dset , commandline ) ;
00984 
00985   
00986   /*----- deallocate memory -----*/   
00987   THD_delete_3dim_dataset (anat_dset, False);   anat_dset = NULL ;
00988 
00989   /*-- 16 Apr 2003 - RWCox:
00990        see if we can convert output back to bytes, if input was bytes --*/
00991 
00992   output_datum = MRI_short ;             /* default, in sfim */
00993 
00994   if( input_datum == MRI_byte ){         /* if input was byte */
00995     short stop = sfim[0] ;
00996     for( ii=1 ; ii < nxyz ; ii++ )
00997       if( sfim[ii] > stop ) stop = sfim[ii] ;
00998     output_datum = MRI_byte ;
00999     bfim = malloc(sizeof(byte)*nxyz) ;
01000     if( stop <= 255 ){                   /* output fits into byte range */
01001       for( ii=0 ; ii < nxyz ; ii++ ) bfim[ii] = (byte) sfim[ii] ;
01002     } else {                             /* must scale output down */
01003       float sfac = 255.9 / stop ;
01004       fprintf(stderr,"++ WARNING: scaling by %g back down to byte data\n",sfac);
01005       for( ii=0 ; ii < nxyz ; ii++ ) bfim[ii] = (byte)(sfim[ii]*sfac) ;
01006     }
01007     free(sfim) ;
01008   }
01009  
01010   /*-- we now return control to your regular programming --*/ 
01011   ibuf[0] = output_datum ;
01012   
01013   ierror = EDIT_dset_items( new_dset ,
01014                             ADN_prefix , filename ,
01015                             ADN_label1 , filename ,
01016                             ADN_self_name , filename ,
01017                             ADN_datum_array , ibuf ,
01018                             ADN_malloc_type, DATABLOCK_MEM_MALLOC ,  
01019                             ADN_none ) ;
01020 
01021      
01022   if( ierror > 0 ){
01023     fprintf(stderr,
01024             "*** %d errors in attempting to create output dataset!\n", 
01025             ierror ) ;
01026     exit(1) ;
01027   }
01028 
01029 
01030   if( THD_is_file(new_dset->dblk->diskptr->header_name) ){
01031     fprintf(stderr,
01032             "*** Output dataset file %s already exists--cannot continue!\a\n",
01033             new_dset->dblk->diskptr->header_name ) ;
01034     exit(1) ;
01035   }
01036   
01037   
01038   /*----- attach bricks to new data set -----*/
01039 
01040   if( output_datum == MRI_short )
01041     mri_fix_data_pointer (sfim, DSET_BRICK(new_dset,0)); 
01042   else if( output_datum == MRI_byte )
01043     mri_fix_data_pointer (bfim, DSET_BRICK(new_dset,0));    /* 16 Apr 2003 */
01044 
01045   fimfac = 1.0;
01046 
01047   /*----- write afni data set -----*/
01048   if (!quiet)
01049     {
01050       printf ("\nWriting anatomical dataset: ");
01051       printf("%s\n", DSET_BRIKNAME(new_dset) ) ;
01052       printf("data type = %s\n",MRI_TYPE_name[output_datum]) ;
01053     }
01054 
01055 
01056   for( ii=0 ; ii < MAX_STAT_AUX ; ii++ ) fbuf[ii] = 0.0 ;
01057   (void) EDIT_dset_items( new_dset , ADN_stat_aux , fbuf , ADN_none ) ;
01058   fbuf[0] = (output_datum == MRI_short && fimfac != 1.0 ) ? fimfac : 0.0 ;
01059   (void) EDIT_dset_items( new_dset , ADN_brick_fac , fbuf , ADN_none ) ;
01060   THD_load_statistics( new_dset ) ;
01061   THD_write_3dim_dataset( NULL,NULL , new_dset , True ) ;
01062 
01063   
01064   /*----- deallocate memory -----*/   
01065   THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
01066   
01067 }
01068 
01069 
01070 /*---------------------------------------------------------------------------*/
01071 /*
01072   This is the main routine for program 3dUniformize.
01073 */
01074 
01075 int main
01076 (
01077   int argc,                /* number of input arguments */
01078   char ** argv             /* array of input arguments */ 
01079 )
01080 
01081 {
01082   UN_options * option_data = NULL;     /* uniformization program options */
01083   short * sfim = NULL;                 /* output uniformized image */
01084 
01085 
01086   { int ii ;                           /* 16 Apr 2003 */
01087     for( ii=1 ; ii < argc ; ii++ ){
01088       if( strcmp(argv[ii],"-quiet") == 0 ){ quiet = 1; break; }
01089     }
01090   }
01091 
01092   /*----- Identify software -----*/
01093   if( !quiet ){
01094    printf ("\n\n");
01095    printf ("Program: %s \n", PROGRAM_NAME);
01096    printf ("Author:  %s \n", PROGRAM_AUTHOR);
01097    printf ("Initial Release:  %s \n", PROGRAM_INITIAL);
01098    printf ("Latest Revision:  %s \n", PROGRAM_LATEST);
01099    printf ("\n");
01100   }
01101 
01102   
01103   /*----- Program initialization -----*/
01104   initialize_program (argc, argv, &option_data, &sfim);
01105 
01106 
01107   /*----- Perform uniformization -----*/
01108   uniformize (option_data, sfim);
01109 
01110 
01111   /*----- Write out the results -----*/
01112   write_afni_data (option_data, sfim);
01113   
01114 
01115   exit(0);
01116 
01117 }
01118 
01119 /*---------------------------------------------------------------------------*/
01120 
01121 
01122 
01123 
01124 
01125