if((ptr)==NULL) \
( UN_error ("Cannot allocate memory") )

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 }

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 }

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 }

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 }

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 }

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 }

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 }

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 }

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 }

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 }

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 }

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 }

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 }

00498 {
00499   int i;
00500 
00501 
00502   for (i = 0;  i < rpts;  i++)
00503     {
00504       c[i] = a[i] - b[i];
00505     }
00506 
00507 }

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 }

00078 {
00079   fprintf (stderr, "%s Error: %s \n", PROGRAM_NAME, message);
00080   exit(1);
00081 }

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 }

00319 {
00320 }

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 }

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 }