00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022 #define PROGRAM_NAME "3dUniformize"
00023 #define PROGRAM_AUTHOR "B. D. Ward"
00024 #define PROGRAM_INITIAL "28 January 2000"
00025 #define PROGRAM_LATEST "16 April 2003"
00026
00027
00028
00029
00030
00031
00032 #include "mrilib.h"
00033 #include "matrix.h"
00034
00035
00036
00037
00038
00039
00040
00041 #define MAX_STRING_LENGTH 80
00042
00043 static THD_3dim_dataset * anat_dset = NULL;
00044 char * commandline = NULL ;
00045
00046 int input_datum = MRI_short ;
00047 int quiet = 0 ;
00048 #define USE_QUIET
00049
00050 typedef struct UN_options
00051 {
00052 char * anat_filename;
00053 char * prefix_filename;
00054 Boolean quiet;
00055 int lower_limit;
00056 int rpts;
00057 int spts;
00058 int nbin;
00059 int npar;
00060 } UN_options;
00061
00062
00063
00064
00065
00066
00067
00068 #include "matrix.c"
00069 #include "estpdf3.c"
00070
00071
00072
00073
00074
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
00087
00088 #define MTEST(ptr) \
00089 if((ptr)==NULL) \
00090 ( UN_error ("Cannot allocate memory") )
00091
00092
00093
00094
00095
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
00120
00121
00122 void initialize_options
00123 (
00124 UN_options * option_data
00125 )
00126
00127 {
00128
00129
00130 option_data->anat_filename = NULL;
00131 option_data->prefix_filename = NULL;
00132 option_data->quiet = FALSE;
00133 option_data->lower_limit = 25;
00134 option_data->rpts = 200000;
00135 option_data->spts = 10000;
00136
00137 option_data->nbin = 250;
00138 option_data->npar = 35;
00139
00140 }
00141
00142
00143
00144
00145
00146
00147
00148 void get_options
00149 (
00150 int argc,
00151 char ** argv,
00152 UN_options * option_data
00153 )
00154
00155 {
00156 int nopt = 1;
00157 int ival, index;
00158 float fval;
00159 char message[MAX_STRING_LENGTH];
00160
00161
00162
00163 if (argc < 2 || strncmp(argv[1], "-help", 5) == 0) display_help_menu();
00164
00165
00166
00167 AFNI_logger (PROGRAM_NAME,argc,argv);
00168
00169
00170
00171 while (nopt < argc )
00172 {
00173
00174
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
00200
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
00230 if (strncmp(argv[nopt], "-quiet", 6) == 0)
00231 {
00232 option_data->quiet = TRUE;
00233 quiet = 1 ;
00234 nopt++;
00235 continue;
00236 }
00237
00238
00239
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
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
00266
00267
00268 void check_one_output_file
00269 (
00270 char * filename
00271 )
00272
00273 {
00274 char message[MAX_STRING_LENGTH];
00275 THD_3dim_dataset * new_dset=NULL;
00276 int ierror;
00277
00278
00279
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
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
00317 )
00318
00319 {
00320 }
00321
00322
00323
00324
00325
00326
00327
00328 void initialize_program
00329 (
00330 int argc,
00331 char ** argv,
00332 UN_options ** option_data,
00333 short ** sfim
00334 )
00335
00336 {
00337 int nxyz;
00338
00339
00340
00341 commandline = tross_commandline( PROGRAM_NAME , argc,argv ) ;
00342
00343
00344
00345 *option_data = (UN_options *) malloc (sizeof(UN_options));
00346 MTEST (*option_data);
00347
00348
00349
00350 initialize_options (*option_data);
00351
00352
00353
00354 get_options (argc, argv, *option_data);
00355
00356
00357
00358 verify_inputs (*option_data);
00359
00360
00361
00362 rand_initialize (1234567);
00363
00364
00365
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
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
00398
00399
00400
00401
00402
00403 void resample
00404 (
00405 UN_options * option_data,
00406 int * ir,
00407 float * vr
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
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
00443
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
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
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
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;
00585 int i, j, k;
00586 matrix x;
00587 matrix xtxinv;
00588 matrix xtxinvxt;
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
00601 matrix_initialize (&x);
00602 matrix_initialize (&xtxinv);
00603 matrix_initialize (&xtxinvxt);
00604 vector_initialize (&y);
00605 vector_initialize (&coef);
00606
00607
00608
00609 matrix_create (spts, p, &x);
00610 vector_create (spts, &y);
00611 xrow = (float *) malloc (sizeof(float) * p);
00612
00613
00614
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
00628
00629
00630
00631
00632 {
00633
00634 matrix xt, xtx;
00635 int ok;
00636
00637
00638
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
00657 matrix_destroy (&xtx);
00658 matrix_destroy (&xt);
00659
00660 }
00661
00662
00663
00664
00665
00666
00667
00668
00669 vector_multiply (xtxinvxt, y, &coef);
00670
00671
00672
00673
00674
00675 for (ip = 0; ip < p; ip++)
00676 {
00677 fpar[ip] = coef.elts[ip];
00678 }
00679
00680
00681
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
00695
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
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];
00750 Boolean ok = TRUE;
00751 char filename[MAX_STRING_LENGTH];
00752
00753
00754
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
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
00775 for (ip = 0; ip < npar; ip++)
00776 {
00777 fpar[ip] = 0.0;
00778 gpar[ip] = 0.0;
00779 }
00780
00781
00782
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
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
00803 for (iter = 1; iter <= 5; iter++)
00804 {
00805
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
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
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
00827 subtract (rpts, vr, fs, ur);
00828 }
00829
00830
00831
00832 for (ip = 0; ip < npar; ip++)
00833 fpar[ip] += gpar[ip];
00834
00835
00836
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
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
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
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
00915 rpts = option_data->rpts;
00916 npar = option_data->npar;
00917
00918
00919
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
00926 resample (option_data, ir, vr);
00927
00928
00929
00930 estimate_field (option_data, ir, vr, fpar);
00931
00932
00933
00934 remove_field (option_data, fpar, sfim);
00935
00936
00937
00938 free (ir); ir = NULL;
00939 free (vr); vr = NULL;
00940 free (fpar); fpar = NULL;
00941
00942 }
00943
00944
00945
00946
00947
00948
00949
00950
00951
00952 void write_afni_data
00953 (
00954 UN_options * option_data,
00955 short * sfim
00956 )
00957
00958 {
00959 int nxyz;
00960 int ii;
00961 THD_3dim_dataset * new_dset=NULL;
00962 int ierror;
00963 int ibuf[32];
00964 float fbuf[MAX_STAT_AUX];
00965 float fimfac;
00966 int output_datum;
00967 char * filename;
00968 byte *bfim = NULL ;
00969
00970
00971
00972 filename = option_data->prefix_filename;
00973 nxyz = DSET_NX(anat_dset) * DSET_NY(anat_dset) * DSET_NZ(anat_dset);
00974
00975
00976
00977 new_dset = EDIT_empty_copy( anat_dset ) ;
00978
00979
00980
00981 tross_Copy_History( anat_dset , new_dset ) ;
00982 if( commandline != NULL )
00983 tross_Append_History( new_dset , commandline ) ;
00984
00985
00986
00987 THD_delete_3dim_dataset (anat_dset, False); anat_dset = NULL ;
00988
00989
00990
00991
00992 output_datum = MRI_short ;
00993
00994 if( input_datum == MRI_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 ){
01001 for( ii=0 ; ii < nxyz ; ii++ ) bfim[ii] = (byte) sfim[ii] ;
01002 } else {
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
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
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));
01044
01045 fimfac = 1.0;
01046
01047
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
01065 THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
01066
01067 }
01068
01069
01070
01071
01072
01073
01074
01075 int main
01076 (
01077 int argc,
01078 char ** argv
01079 )
01080
01081 {
01082 UN_options * option_data = NULL;
01083 short * sfim = NULL;
01084
01085
01086 { int ii ;
01087 for( ii=1 ; ii < argc ; ii++ ){
01088 if( strcmp(argv[ii],"-quiet") == 0 ){ quiet = 1; break; }
01089 }
01090 }
01091
01092
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
01104 initialize_program (argc, argv, &option_data, &sfim);
01105
01106
01107
01108 uniformize (option_data, sfim);
01109
01110
01111
01112 write_afni_data (option_data, sfim);
01113
01114
01115 exit(0);
01116
01117 }
01118
01119
01120
01121
01122
01123
01124
01125