00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006 
00007 /*---------------------------------------------------------------------------*/
00008 /*
00009   This program performs 2d image registration of slices contained in an AFNI
00010   3d+time dataset.  This program was adapted from plug_imreg.c and imreg.c.
00011 
00012 
00013   File:    2dImReg.c
00014   Author:  B. Douglas Ward
00015   Date:    04 February 1998
00016 
00017 
00018   Mod:     Added routines to write the registration parameters, and the RMS 
00019            error, to user specified ASCII files.
00020   Date:    20 March 1998
00021 
00022   Mod:     Added option to change dx and dy output format from pixels to mm.
00023   Date:    24 March 1998
00024 
00025   Mod:     No longer assume that input and base datasets have same length.
00026            This problem was reported by Paul Reber.
00027   Date:    3 June 1998
00028 
00029   Mod:     Routine eval_registration extended to include byte and float datum
00030            types for base and input datasets.
00031   Date:    3 June 1998
00032 
00033   Mod:     Corrected problem with base image memory deallocation.
00034   Date:    20 July 1998
00035 
00036   Mod:     Added changes for incorporating History notes.
00037   Date:    10 September 1999
00038 
00039   Mod:     Set MAX_NAME_LENGTH equal to THD_MAX_NAME.
00040   Date:    02 December 2002
00041 
00042   Mod:     If nsl == 0, set num_slices equal to nz.   [4 occurrences]
00043   Date:    06 October 2003  [rickr]
00044 */
00045 
00046 /*---------------------------------------------------------------------------*/
00047 
00048 #define PROGRAM_NAME    "2dImReg"                   /* name of this program */
00049 #define PROGRAM_INITIAL "04 Feb 1998"     /* date of initial program release */
00050 #define PROGRAM_LATEST  "02 Dec 2002"     /* date of latest program revision */
00051 
00052 /*---------------------------------------------------------------------------*/
00053 
00054 #include "mrilib.h"
00055 #include "matrix.h"
00056 
00057 #define MAX_NAME_LENGTH THD_MAX_NAME  /* max. string length for file names */ 
00058 #define STATE_DIM 4                   /* number of registration parameters */ 
00059 
00060 
00061 /*----- Global variables -----*/ 
00062 int * t_to_z = NULL;           /* convert time-order to z-order of slices */  
00063 int * z_to_t = NULL;           /* convert z-order to time-order of slices */
00064 
00065 
00066 static char * commandline = NULL ;       /* command line for history notes */
00067  
00068 
00069 typedef struct IR_options      /* user input options */
00070 {
00071   char * input_filename;       /* file name for input 3d+time dataset */
00072   char * base_filename;        /* file name for reference (base) volume */
00073   int base_vol_index;          /* image number for base volume */
00074   int nofine;                  /* boolean for no fine fit */
00075   float blur;                  /* FWHM of blurring prior to registration */
00076   float dxy;                   /* convergence tolerance for translations */
00077   float dphi;                  /* convergence tolerance for rotations */
00078   char * new_prefix;           /* prefix name for registered dataset */
00079   char * dprefix;              /* prefix name for registration parameters */
00080   int dmm;                     /* change dx and dy output from pixels to mm */
00081   char * rprefix;              /* prefix name for volume RMS error */
00082   int debug;                   /* write additional output to screen */
00083 } IR_options;
00084 
00085 
00086 #include "matrix.c"
00087 
00088 
00089 /*---------------------------------------------------------------------------*/
00090 /*
00091   Routine to display 2dImReg help menu.
00092 */
00093 
00094 void display_help_menu()
00095 {
00096   printf 
00097     (
00098      "This program performs 2d image registration.  Image alignment is      \n"
00099      "performed on a slice-by-slice basis for the input 3d+time dataset,    \n"
00100      "relative to a user specified base image.                              \n"
00101      "                                                                      \n"
00102      "Usage:                                                                \n"
00103      "2dImReg                                                               \n"
00104      "-input fname           Filename of input 3d+time dataset to process   \n"
00105      "-basefile fname        Filename of 3d+time dataset for base image     \n"
00106      "                         (default = current input dataset)            \n"
00107      "-base num              Time index for base image  (0 <= num)          \n"
00108      "                         (default:  num = 3)                          \n"
00109      "-nofine                Deactivate fine fit phase of image registration\n"
00110      "                         (default:  fine fit is active)               \n"
00111      "-fine blur dxy dphi    Set fine fit parameters                        \n"
00112      "   where:                                                             \n"
00113      "     blur = FWHM of blurring prior to registration (in pixels)        \n"
00114      "               (default:  blur = 1.0)                                 \n"
00115      "     dxy  = Convergence tolerance for translations (in pixels)        \n"
00116      "               (default:  dxy  = 0.07)                                \n"
00117      "     dphi = Convergence tolerance for rotations (in degrees)          \n"
00118      "               (default:  dphi = 0.21)                                \n"
00119      "                                                                      \n"
00120      "-prefix pname     Prefix name for output 3d+time dataset              \n"
00121      "                                                                      \n"
00122      "-dprefix dname    Write files 'dname'.dx, 'dname'.dy, 'dname'.psi     \n"
00123      "                    containing the registration parameters for each   \n"
00124      "                    slice in chronological order.                     \n"
00125      "                    File formats:                                     \n"
00126      "                      'dname'.dx:    time(sec)   dx(pixels)           \n"
00127      "                      'dname'.dy:    time(sec)   dy(pixels)           \n"
00128      "                      'dname'.psi:   time(sec)   psi(degrees)         \n"
00129      "-dmm              Change dx and dy output format from pixels to mm    \n"
00130      "                                                                      \n"
00131      "-rprefix rname    Write files 'rname'.oldrms and 'rname'.newrms       \n"
00132      "                    containing the volume RMS error for the original  \n"
00133      "                    and the registered datasets, respectively.        \n"
00134      "                    File formats:                                     \n"
00135      "                      'rname'.oldrms:   volume(number)   rms_error    \n"
00136      "                      'rname'.newrms:   volume(number)   rms_error    \n"
00137      "                                                                      \n"
00138      "-debug            Lots of additional output to screen                 \n"
00139     );
00140   
00141   exit(0);
00142 }
00143 
00144 
00145 /*---------------------------------------------------------------------------*/
00146 /*
00147    Print error message and stop.
00148 */
00149 
00150 void IR_error 
00151 (
00152   char * message               /* error message to be displayed */
00153 )
00154 
00155 {
00156   fprintf (stderr, "%s Error: %s \n", PROGRAM_NAME, message);
00157   exit(1);
00158 }
00159 
00160 
00161 /*---------------------------------------------------------------------------*/
00162 /*
00163   Routine to initialize the input options.
00164 */
00165 
00166 void initialize_options
00167 (
00168   IR_options ** opt            /* user input options */
00169 )
00170 
00171 {
00172   (*opt) = (IR_options *) malloc (sizeof(IR_options));
00173 
00174   (*opt)->input_filename = NULL;
00175   (*opt)->base_filename = NULL;
00176   (*opt)->base_vol_index = 3;       
00177   (*opt)->nofine = 0;
00178   (*opt)->blur = 1.0;     
00179   (*opt)->dxy  = 0.07;    
00180   (*opt)->dphi = 0.21;   
00181   (*opt)->new_prefix = NULL;
00182   (*opt)->dprefix = NULL;
00183   (*opt)->rprefix = NULL;
00184   (*opt)->dmm = 0;
00185   (*opt)->debug = 0;
00186 }
00187 
00188 
00189 /*---------------------------------------------------------------------------*/
00190 /*
00191    Routine to get user specified input options.
00192 */
00193 
00194 void get_user_inputs 
00195 (
00196   int argc,                       /* number of input arguments */
00197   char ** argv,                   /* array of input arguments */ 
00198   IR_options ** option_data       /* user input options */
00199 )
00200 
00201 {
00202   int nopt = 1;                   /* input option argument counter */
00203   int ival;
00204   float fval;
00205   char message[80];               /* error message */
00206 
00207   
00208   /*----- does user request help menu? -----*/
00209   if (argc < 2 || strncmp(argv[1], "-help", 5) == 0)  display_help_menu();  
00210   
00211 
00212    /*----- main loop over input options -----*/
00213   while (nopt < argc )
00214     {
00215 
00216       /*-----   -input filename   -----*/
00217       if (strncmp(argv[nopt], "-input", 6) == 0)
00218         {
00219           nopt++;
00220           if (nopt >= argc)  IR_error ("Need argument after -input ");
00221           (*option_data)->input_filename 
00222             = (char *) malloc (sizeof(char) * MAX_NAME_LENGTH);
00223           strcpy ((*option_data)->input_filename, argv[nopt]);
00224           nopt++;
00225           continue;
00226         }
00227      
00228 
00229       /*-----   -basefile filename   -----*/
00230       if (strncmp(argv[nopt], "-basefile", 9) == 0)
00231         {
00232           nopt++;
00233           if (nopt >= argc)  IR_error ("Need argument after -basefile ");
00234           (*option_data)->base_filename 
00235             = (char *) malloc (sizeof(char) * MAX_NAME_LENGTH);
00236           strcpy ((*option_data)->base_filename, argv[nopt]);
00237           nopt++;
00238           continue;
00239         }
00240      
00241 
00242       /*-----   -base num  -----*/
00243       if (strncmp(argv[nopt], "-base", 5) == 0)
00244         {
00245           nopt++;
00246           if (nopt >= argc)  IR_error ("Need argument after -base ");
00247           sscanf (argv[nopt], "%d", &ival);
00248           if (ival < 0) 
00249             IR_error ("Illegal argument after -base  ( must be >= 0 ) ");
00250           (*option_data)->base_vol_index = ival;
00251           nopt++;
00252           continue;
00253         }
00254 
00255 
00256       /*-----   -prefix pname   -----*/
00257       if (strncmp(argv[nopt], "-prefix", 7) == 0)
00258         {
00259           nopt++;
00260           if (nopt >= argc)  IR_error ("Need argument after -prefix ");
00261           (*option_data)->new_prefix 
00262             = (char *) malloc (sizeof(char) * MAX_NAME_LENGTH);
00263           strcpy ((*option_data)->new_prefix, argv[nopt]);
00264 
00265           nopt++;
00266           continue;
00267         }
00268      
00269 
00270       /*-----   -dprefix dname   -----*/
00271       if (strncmp(argv[nopt], "-dprefix", 8) == 0)
00272         {
00273           nopt++;
00274           if (nopt >= argc)  IR_error ("Need argument after -dprefix ");
00275           (*option_data)->dprefix 
00276             = (char *) malloc (sizeof(char) * MAX_NAME_LENGTH);
00277           strcpy ((*option_data)->dprefix, argv[nopt]);
00278 
00279           nopt++;
00280           continue;
00281         }
00282      
00283 
00284       /*-----   -rprefix rname   -----*/
00285       if (strncmp(argv[nopt], "-rprefix", 8) == 0)
00286         {
00287           nopt++;
00288           if (nopt >= argc)  IR_error ("Need argument after -rprefix ");
00289           (*option_data)->rprefix 
00290             = (char *) malloc (sizeof(char) * MAX_NAME_LENGTH);
00291           strcpy ((*option_data)->rprefix, argv[nopt]);
00292 
00293           nopt++;
00294           continue;
00295         }
00296      
00297 
00298       /*-----   -nofine -----*/
00299       if (strncmp(argv[nopt], "-nofine", 7) == 0)
00300         {
00301           (*option_data)->nofine = 1;
00302           nopt++;
00303           continue;
00304         }
00305 
00306            
00307       /*-----   -fine blur dxy dphi  -----*/
00308       if (strncmp(argv[nopt], "-fine", 5) == 0)
00309         {
00310           nopt++;
00311           if (nopt+2 >= argc)  IR_error ("Need 3 arguments after -fine ");
00312           (*option_data)->nofine = 0;
00313 
00314           sscanf (argv[nopt], "%f", &fval);
00315           if (fval <= 0.0) 
00316             IR_error ("Illegal argument for blur  ( must be > 0 ) ");
00317           (*option_data)->blur = fval;
00318           nopt++;
00319 
00320           sscanf (argv[nopt], "%f", &fval);
00321           if (fval <= 0.0)
00322             IR_error ("Illegal argument for dxy  ( must be > 0 ) ");
00323           (*option_data)->dxy = fval;
00324           nopt++;
00325 
00326           sscanf (argv[nopt], "%f", &fval);
00327           if (fval <= 0.0)
00328             IR_error ("Illegal argument for dphi  ( must be > 0 ) ");
00329           (*option_data)->dphi = fval;
00330           nopt++;
00331 
00332           continue;
00333         }
00334       
00335 
00336       /*-----   -dmm -----*/
00337       if (strncmp(argv[nopt], "-dmm", 4) == 0)
00338         {
00339           (*option_data)->dmm = 1;
00340           nopt++;
00341           continue;
00342         }
00343 
00344            
00345       /*-----   -debug -----*/
00346       if (strncmp(argv[nopt], "-debug", 6) == 0)
00347         {
00348           (*option_data)->debug = 1;
00349           nopt++;
00350           continue;
00351         }
00352 
00353            
00354       /*----- unknown command -----*/
00355       sprintf(message,"Unrecognized command line option: %s\n", argv[nopt]);
00356       IR_error (message);
00357 
00358 
00359     }
00360       
00361  
00362 }
00363 
00364 
00365 /*---------------------------------------------------------------------------*/
00366 /*
00367    Routine to read a 3d+time dataset.
00368 */
00369 
00370 void read_dataset 
00371 (
00372   char * filename,                /* file name of 3d+time dataset */
00373   THD_3dim_dataset ** dset        /* pointer to 3d+time dataset */
00374 )
00375 
00376 {
00377   char message[80];               /* error message */
00378 
00379 
00380   /*----- Open the 3d+time dataset -----*/
00381   *dset = THD_open_one_dataset (filename);
00382   if (*dset == NULL)  
00383     { 
00384       sprintf (message, 
00385                "Unable to open data file: %s", filename);
00386       IR_error (message);
00387     }
00388 
00389 }
00390 
00391 
00392 /*---------------------------------------------------------------------------*/
00393 /*
00394    Initialize the slice sequence arrays.
00395 */
00396 
00397 void initialize_slice_sequence 
00398 (
00399   IR_options * option_data,        /* user input options */
00400   THD_3dim_dataset * dset          /* pointer to 3d+time dataset */
00401 )
00402 
00403 {
00404   int num_slices;                  /* number of slices per volume */
00405   int ivolume;                     /* volume index number */
00406   int itemp;                       /* temporary variable */
00407   float ttemp;                     /* temporary variable */
00408   float * time_array = NULL;       /* array of slice acquisition times */
00409   int iz, i, j;                    /* index numbers */
00410   float z;                         /* slice z location */
00411 
00412 
00413   /*----- Initialize local variables -----*/
00414   num_slices = dset->taxis->nsl;
00415   ivolume = 0;
00416 
00417   if ( num_slices <= 0 )            /* 06 Oct 2003 [rickr] */
00418       num_slices = dset->daxes->nzz;
00419 
00420   /*----- Allocate memory for arrays -----*/
00421   t_to_z = (int *) malloc (sizeof(int) * num_slices);
00422   z_to_t = (int *) malloc (sizeof(int) * num_slices);
00423   time_array = (float *) malloc (sizeof(float) * num_slices);
00424 
00425 
00426   /*----- Initialize array of slice acquisition times -----*/
00427   for (iz = 0;  iz < num_slices;  iz++)
00428     {
00429       z = iz * dset->taxis->dz_sl + dset->taxis->zorg_sl;
00430       time_array[iz] = THD_timeof (ivolume, z, dset->taxis);
00431       t_to_z[iz] = iz;
00432     }
00433 
00434 
00435   /*----- Sort slice z-indices by increasing time -----*/
00436   for (i = 0;  i < num_slices-1;  i++)
00437     for (j = i+1;  j < num_slices;  j++)
00438       if (time_array[j] < time_array[i])
00439         {
00440           itemp = t_to_z[i];
00441           t_to_z[i] = t_to_z[j];
00442           t_to_z[j] = itemp;
00443 
00444           ttemp = time_array[i];
00445           time_array[i] = time_array[j];
00446           time_array[j] = ttemp;
00447         } 
00448 
00449 
00450   /*----- Sort slice time-indices by increasing z index -----*/
00451   for (i = 0;  i < num_slices;  i++)
00452     {
00453       j = t_to_z[i];
00454       z_to_t[j] = i;
00455     }
00456 
00457 
00458   /*----- Write out the slice ordering arrays -----*/
00459   if (option_data->debug)
00460     for (i = 0;  i < num_slices;  i++)
00461       printf ("time[%2d] = %12.3f   t_to_z[%2d] = %2d   z_to_t[%2d] = %2d\n",
00462               i, time_array[i], i, t_to_z[i], i, z_to_t[i]);
00463 
00464   
00465   /*----- Release memory -----*/
00466   free (time_array);   time_array = NULL;
00467 } 
00468 
00469 
00470 /*---------------------------------------------------------------------------*/
00471 /*
00472    Routine to initialize the array of state vectors.
00473 */
00474 
00475 void initialize_state_history 
00476 (
00477   THD_3dim_dataset * dset,            /* pointer to input 3d+time dataset */
00478   vector ** state_history             /* time series of state vectors */
00479 )
00480 
00481 {
00482   int num_slices;                     /* number of slices per volume */
00483   int ts_length;                      /* number of volumes */
00484   int num_vectors;                    /* total number of state vectors */
00485   int i;                              /* state vector index */
00486 
00487 
00488   /*----- Initialize local variables -----*/
00489   num_slices = dset->taxis->nsl;
00490 
00491   if ( num_slices <= 0 )              /* 06 Oct 2003 [rickr] */
00492       num_slices = dset->daxes->nzz;
00493 
00494   ts_length = DSET_NUM_TIMES(dset);
00495   num_vectors = ts_length * num_slices;
00496 
00497 
00498   /*----- Allocate memory for array of state vectors -----*/
00499   *state_history = (vector *) malloc (sizeof(vector) * num_vectors);
00500 
00501 
00502   /*----- Initialize array of state vectors -----*/
00503   for (i = 0;  i < num_vectors;  i++)
00504     {
00505       vector_initialize (&((*state_history)[i]));
00506       vector_create (STATE_DIM, &((*state_history)[i]));
00507     }
00508 }
00509 
00510  
00511 /*---------------------------------------------------------------------------*/
00512 /*
00513    Routine to initialize the RMS error arrays.
00514 */
00515 
00516 void initialize_rms_arrays 
00517 (
00518   THD_3dim_dataset * dset,       /* pointer to input 3d+time dataset */
00519   float ** old_rms_array,        /* volume RMS error for input dataset */
00520   float ** new_rms_array         /* volume RMS error for registered dataset */
00521 )
00522 
00523 {
00524   int ts_length;                 /* number of volumes */
00525 
00526 
00527   ts_length = DSET_NUM_TIMES(dset);
00528 
00529   /*----- Allocate space for RMS error arrays -----*/
00530   *old_rms_array = (float *) malloc (sizeof(float) * ts_length);
00531   *new_rms_array = (float *) malloc (sizeof(float) * ts_length);
00532   
00533 
00534 }
00535 
00536  
00537 /*---------------------------------------------------------------------------*/
00538 /*
00539   Routine to perform all program initialization.
00540 */
00541 
00542 void initialize_program
00543 (
00544   int argc,                       /* number of input arguments */
00545   char ** argv,                   /* array of input arguments */ 
00546   IR_options ** option_data,      /* user input options */
00547   vector ** state_history,        /* time series of state vectors */
00548   float ** old_rms_array,         /* volume RMS error for input dataset */
00549   float ** new_rms_array          /* volume RMS error for registered dataset */
00550 )
00551 
00552 {
00553   THD_3dim_dataset * dset = NULL;
00554 
00555 
00556   /*----- save command line for history notes -----*/
00557   commandline = tross_commandline( PROGRAM_NAME , argc,argv ) ;
00558 
00559 
00560   /*----- Initialize input options -----*/
00561   initialize_options (option_data);
00562 
00563 
00564   /*----- Get user inputs -----*/
00565   get_user_inputs (argc, argv, option_data);
00566 
00567 
00568   /*----- Read the input 3d+time dataset to be registered -----*/
00569   read_dataset ((*option_data)->input_filename, &dset);
00570 
00571   
00572   /*----- Initialize the z-slice time order arrays -----*/
00573   initialize_slice_sequence (*option_data, dset);
00574 
00575 
00576   /*----- Initialize the array of state vectors -----*/
00577   if ((*option_data)->dprefix != NULL)
00578     initialize_state_history (dset, state_history);
00579 
00580 
00581   /*----- Allocate space for RMS error arrays -----*/
00582   if ((*option_data)->rprefix != NULL)
00583     initialize_rms_arrays (dset, old_rms_array, new_rms_array);
00584 
00585 
00586   /*----- Release memory -----*/
00587   THD_delete_3dim_dataset (dset, False);   dset = NULL;
00588 
00589 }
00590 
00591 
00592 /*---------------------------------------------------------------------------*/
00593 /*
00594    Routine to make a copy of a dataset, with data attached.
00595    This routine is copied directly from afni_plugin.c
00596 */
00597 
00598 THD_3dim_dataset * copy_dset( THD_3dim_dataset * dset , char * new_prefix )
00599 {
00600    THD_3dim_dataset * new_dset ;
00601    int ival , ityp , nbytes , nvals ;
00602    void * new_brick , * old_brick ;
00603 
00604    /*-- sanity check --*/
00605 
00606    if( ! ISVALID_3DIM_DATASET(dset) ) return NULL ;
00607 
00608    /*-- make the empty copy --*/
00609 
00610    new_dset = EDIT_empty_copy( dset ) ;
00611 
00612    /*-- change its name? --*/
00613 
00614    if( new_prefix != NULL )
00615       EDIT_dset_items( new_dset ,
00616                           ADN_prefix , new_prefix ,
00617                           ADN_label1 , new_prefix ,
00618                        ADN_none ) ;
00619 
00620    /*-- make brick(s) for this dataset --*/
00621 
00622    THD_load_datablock( dset->dblk ) ;  /* make sure old one is in memory */
00623 
00624    nvals = DSET_NVALS(dset) ;
00625 
00626    for( ival=0 ; ival < nvals ; ival++ ){
00627       ityp      = DSET_BRICK_TYPE(new_dset,ival) ;   /* type of data */
00628       nbytes    = DSET_BRICK_BYTES(new_dset,ival) ;  /* how much data */
00629       new_brick = malloc( nbytes ) ;                 /* make room */
00630 
00631       if( new_brick == NULL ){
00632         THD_delete_3dim_dataset( new_dset , False ) ;
00633         return NULL ;
00634       }
00635 
00636       EDIT_substitute_brick( new_dset , ival , ityp , new_brick ) ;
00637 
00638       /*-- copy data from old brick to new brick --*/
00639 
00640       old_brick = DSET_BRICK_ARRAY(dset,ival) ;
00641 
00642       if( old_brick == NULL ){
00643          THD_delete_3dim_dataset( new_dset , False ) ;
00644          return NULL ;
00645       }
00646 
00647       memcpy( new_brick , old_brick , nbytes ) ;
00648    }
00649 
00650    return new_dset ;
00651 }
00652 
00653 
00654 /*---------------------------------------------------------------------------*/
00655 /*
00656   Check whether output file already exists.
00657 */
00658 
00659 void check_output_file 
00660 (
00661   THD_3dim_dataset * dset,      /* input afni data set pointer */
00662   char * filename               /* output file name */
00663 )
00664 
00665 {
00666   THD_3dim_dataset * new_dset=NULL;   /* output afni data set pointer */
00667   int ierror;                         /* number of errors in editing data */
00668   
00669   
00670   /*-- make an empty copy of the input dataset --*/
00671   new_dset = EDIT_empty_copy( dset ) ;
00672   
00673   
00674   ierror = EDIT_dset_items( new_dset ,
00675                             ADN_prefix , filename ,
00676                             ADN_label1 , filename ,
00677                             ADN_self_name , filename ,
00678                             ADN_type , ISHEAD(dset) ? HEAD_FUNC_TYPE : 
00679                                                       GEN_FUNC_TYPE ,
00680                             ADN_none ) ;
00681   
00682   if( ierror > 0 ){
00683     fprintf(stderr,
00684             "*** %d errors in attempting to create output dataset!\n", ierror ) ;
00685     exit(1) ;
00686   }
00687   
00688   if( THD_is_file(new_dset->dblk->diskptr->header_name) ){
00689     fprintf(stderr,
00690             "*** Output dataset file %s already exists--cannot continue!\a\n",
00691             new_dset->dblk->diskptr->header_name ) ;
00692     exit(1) ;
00693   }
00694   
00695   /*----- deallocate memory -----*/   
00696   THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
00697   
00698 }
00699 
00700 
00701 /*---------------------------------------------------------------------------*/
00702 /*
00703   Evaluate accuracy of registration.
00704 */
00705 
00706 void eval_registration
00707 (
00708   IR_options * option_data,       /* user input options */
00709   THD_3dim_dataset * old_dset,    /* input dataset */
00710   THD_3dim_dataset * new_dset,    /* output datasets */
00711   THD_3dim_dataset * base_dset,   /* base image dataset */
00712   int base,                       /* base volume index */
00713   float * old_rms_array,          /* volume RMS error for input dataset */
00714   float * new_rms_array           /* volume RMS error for registered dataset */
00715   )
00716      
00717 {
00718   int nx, ny, nz, nxyz;
00719   int ix, jy, kz;
00720   int ixyz;
00721   int datum, base_datum;
00722   float old_sse, new_sse;
00723   float diff;
00724   float old_rmse, new_rmse;
00725   int ivolume, num_volumes;
00726   byte  * bold = NULL, * bnew = NULL, * bbase = NULL;
00727   short * sold = NULL, * snew = NULL, * sbase = NULL;
00728   float * fold = NULL, * fnew = NULL, * fbase = NULL;
00729   float float_base, float_old, float_new;
00730 
00731 
00732   /*----- Initialize local variables -----*/
00733   nx = old_dset->daxes->nxx;
00734   ny = old_dset->daxes->nyy;
00735   nz = old_dset->daxes->nzz;
00736   nxyz = nx * ny * nz;
00737   num_volumes = DSET_NUM_TIMES (old_dset);
00738   datum       = DSET_BRICK_TYPE (new_dset,0) ;
00739   base_datum  = DSET_BRICK_TYPE (base_dset,0);
00740 
00741 
00742   /*----- Set base dataset pointer depending on base datum type -----*/
00743   switch ( base_datum )
00744     { 
00745     case MRI_byte:
00746       bbase = (byte *) DSET_ARRAY(base_dset,base);  break;
00747       
00748     case MRI_short:
00749       sbase = (short *) DSET_ARRAY(base_dset,base);  break;
00750 
00751     case MRI_float:
00752       fbase = (float *) DSET_ARRAY(base_dset,base);  break;
00753     }
00754 
00755 
00756   for (ivolume = 0;  ivolume < num_volumes; ivolume++)
00757     {
00758       old_sse = 0.0;
00759       new_sse = 0.0;
00760 
00761       /*----- Set old and new dataset pointers depending on datum type -----*/
00762       switch ( datum )
00763         {  
00764         case MRI_byte:
00765           bold = (byte *) DSET_ARRAY(old_dset,ivolume);  
00766           bnew = (byte *) DSET_ARRAY(new_dset,ivolume);  break;
00767           
00768         case MRI_short:
00769           sold = (short *) DSET_ARRAY(old_dset,ivolume);  
00770           snew = (short *) DSET_ARRAY(new_dset,ivolume);  break;
00771           
00772         case MRI_float:
00773           fold = (float *) DSET_ARRAY(old_dset,ivolume);  
00774           fnew = (float *) DSET_ARRAY(new_dset,ivolume);  break;
00775         }
00776 
00777       
00778       for (kz = 0;  kz < nz;  kz++)
00779         for (jy = 0;  jy < ny;  jy++)
00780           for (ix = 0;  ix < nx;  ix++)
00781             {
00782               ixyz = ix + jy*nx + kz*nx*ny;
00783 
00784               /*----- Get base voxel floating point value -----*/
00785               switch ( base_datum )
00786                 { 
00787                 case MRI_byte:   float_base = (float) bbase[ixyz];  break;
00788                   
00789                 case MRI_short:  float_base = (float) sbase[ixyz];  break;
00790                   
00791                 case MRI_float:  float_base = fbase[ixyz];  break;
00792                 }
00793               
00794               /*----- Get old and new voxel floating point value -----*/
00795               switch ( datum )
00796                 {  
00797                 case MRI_byte:
00798                   float_old = (float) bold[ixyz];  
00799                   float_new = (float) bnew[ixyz];  break;
00800                   
00801                 case MRI_short:
00802                   float_old = (float) sold[ixyz];  
00803                   float_new = (float) snew[ixyz];  break;
00804                   
00805                 case MRI_float:
00806                   float_old = fold[ixyz];  
00807                   float_new = fnew[ixyz];  break;
00808                 }
00809 
00810               diff = float_old - float_base;
00811               old_sse += diff*diff;
00812               diff = float_new - float_base;
00813               new_sse += diff*diff;
00814             }
00815       
00816       old_rmse = sqrt (old_sse / nxyz);
00817       new_rmse = sqrt (new_sse / nxyz);
00818 
00819       if (option_data->debug)
00820         printf ("Volume = %d   OLD RMSE = %f   NEW RMSE = %f \n", 
00821                 ivolume, old_rmse, new_rmse);
00822       
00823       old_rms_array[ivolume] = old_rmse;
00824       new_rms_array[ivolume] = new_rmse;
00825 
00826     }
00827 
00828 }
00829 
00830 
00831 /*---------------------------------------------------------------------------*/
00832 /*
00833   Main routine for this program.
00834   If the return string is not NULL, some error transpired, and
00835   the program will display the error message.
00836 */
00837 
00838 #define FREEUP(x) do{if((x) != NULL){free((x)); (x)=NULL;}}while(0)
00839 
00840 #define FREE_WORKSPACE                             \
00841   do{ FREEUP(bptr) ; FREEUP(sptr) ; FREEUP(fptr) ; \
00842       FREEUP(bout) ; FREEUP(sout) ; FREEUP(fout) ; \
00843       FREEUP(dxar) ; FREEUP(dyar) ; FREEUP(phiar); \
00844     } while(0) ;
00845 
00846 
00847 char * IMREG_main 
00848 (
00849   IR_options * opt,
00850   vector * state_history,
00851   float * old_rms_array,
00852   float * new_rms_array
00853 )
00854 {
00855    THD_3dim_dataset * old_dset , * new_dset ;  /* input and output datasets */
00856    THD_3dim_dataset * base_dset;               /* base image dataset */
00857    char * new_prefix ;                         /* string from user */
00858    int base , ntime , datum , nx,ny,nz , ii,kk , npix ;
00859    float                      dx,dy,dz ;
00860    int base_datum;
00861    MRI_IMARR * ims_in , * ims_out ;
00862    MRI_IMAGE * im , * imbase ;
00863 
00864    byte   ** bptr = NULL , * bbase = NULL, ** bout = NULL ;
00865    short  ** sptr = NULL , * sbase = NULL, ** sout = NULL ; 
00866    float  ** fptr = NULL , * fbase = NULL, ** fout = NULL ;
00867 
00868    float * dxar = NULL , * dyar = NULL , * phiar = NULL ;
00869 
00870    int it;
00871 
00872    /*--------------------------------------------------------------------*/
00873    /*----- Check batch command inputs to see if they are reasonable -----*/
00874 
00875    old_dset = THD_open_one_dataset(opt->input_filename) ;   
00876                                                       /* get ptr to dataset */
00877    if( old_dset == NULL )
00878       return "*************************\n"
00879              "Cannot find Input Dataset\n"
00880              "*************************"  ;
00881 
00882    ntime = DSET_NUM_TIMES(old_dset) ;
00883    if( ntime < 2 )
00884       return "*****************************\n"
00885              "Dataset has only 1 time point\n"
00886              "*****************************"  ;
00887 
00888    ii = DSET_NVALS_PER_TIME(old_dset) ;
00889    if( ii > 1 )
00890       return "************************************\n"
00891              "Dataset has > 1 value per time point\n"
00892              "************************************"  ;
00893 
00894    nx = old_dset->daxes->nxx ; dx = old_dset->daxes->xxdel ;
00895    ny = old_dset->daxes->nyy ; dy = old_dset->daxes->yydel ; npix = nx*ny ;
00896    nz = old_dset->daxes->nzz ; dz = old_dset->daxes->zzdel ;
00897 
00898    if( nx != ny || fabs(dx) != fabs(dy) ){
00899 
00900      if (opt->debug)
00901        fprintf(stderr,"\nIMREG: nx=%d ny=%d nz=%d  dx=%f dy=%f dz=%f\n",
00902                nx,ny,nz,dx,dy,dz ) ;
00903 
00904       return "***********************************\n"
00905              "Dataset does not have square slices\n"
00906              "***********************************"  ;
00907    }
00908 
00909    new_prefix = opt->new_prefix;     /* get string item (the output prefix) */
00910    if (new_prefix == NULL)           /* check if it is OK */
00911       return "************************\n"
00912              "Output Prefix is illegal\n"
00913              "************************"  ;
00914 
00915    /*----- Check whether output file already exists -----*/
00916    check_output_file (old_dset, new_prefix);
00917 
00918 
00919    /*--------- go to "base" input option ---------*/
00920 
00921    if (opt->base_filename == NULL)
00922      base_dset = old_dset;
00923    else
00924      {
00925        base_dset = THD_open_one_dataset(opt->base_filename) ;   
00926                                                   /* get ptr to base dataset */
00927        if( base_dset == NULL )
00928          return "************************\n"
00929                 "Cannot find Base Dataset\n"
00930                 "************************"  ;
00931 
00932        if ( (nx != base_dset->daxes->nxx) || (dx != base_dset->daxes->xxdel)
00933          || (ny != base_dset->daxes->nyy) || (dy != base_dset->daxes->yydel)
00934          || (nz != base_dset->daxes->nzz) || (dz != base_dset->daxes->zzdel) )
00935          {
00936            if (opt->debug)
00937                {
00938                  fprintf(stderr,"\nIMREG: Input Dataset:\n");
00939                  fprintf(stderr,"nx=%d ny=%d nz=%d  dx=%f dy=%f dz=%f\n",
00940                      nx,ny,nz,dx,dy,dz ) ;
00941 
00942                  fprintf(stderr,"\nIMREG: Base Dataset:\n");
00943                  fprintf(stderr,"nx=%d ny=%d nz=%d  dx=%f dy=%f dz=%f\n",
00944                          base_dset->daxes->nxx,   base_dset->daxes->nyy,
00945                          base_dset->daxes->nzz,   base_dset->daxes->xxdel,
00946                          base_dset->daxes->yydel, base_dset->daxes->zzdel) ;
00947                }
00948            return "*************************************************\n"
00949                   "Base Dataset is not compatible with Input Dataset\n"
00950                   "*************************************************"  ;
00951          }
00952      }
00953 
00954    base_datum = DSET_BRICK_TYPE(base_dset,0);
00955 
00956    base = opt->base_vol_index;
00957 
00958    if( base >= DSET_NUM_TIMES(base_dset))
00959       return "******************************\n"
00960              "Base image number is too large\n"
00961              "******************************"  ;
00962 
00963 
00964    /*--------- see if the "fine" input option is present --------*/
00965    if (opt->nofine)
00966      mri_align_params( 0 , 0.0,0.0,0.0 , 0.0,0.0,0.0 ) ;
00967    else{
00968       float fsig , fdxy , fdph ;
00969       fsig = opt->blur * 0.42466090;
00970       fdxy = opt->dxy;
00971       fdph = opt->dphi;
00972       mri_align_params( 0 , 0.0,0.0,0.0 , fsig,fdxy,fdph ) ;
00973 
00974       if (opt->debug)
00975         fprintf(stderr,"Set fine params = %f %f %f\n",fsig,fdxy,fdph) ; 
00976    }
00977 
00978    /*------------- ready to compute new dataset -----------*/
00979 
00980    if (opt->debug)
00981      fprintf(stderr,"IMREG: loading dataset\n") ;
00982 
00983 
00984    DSET_load( old_dset ) ;
00985    
00986    if (opt->base_filename != NULL)
00987      DSET_load (base_dset);
00988 
00989    /*** 1) Copy the dataset in toto ***/
00990 
00991    if (opt->debug)
00992      fprintf(stderr,"IMREG: Copying dataset\n") ;
00993 
00994 
00995    new_dset = copy_dset( old_dset , new_prefix ) ;
00996    if( new_dset == NULL )
00997       return "****************************\n"
00998              "Failed to copy input dataset\n"
00999              "****************************"  ;
01000 
01001    /*** 2) Make an array of empty images ***/
01002 
01003    if (opt->debug)
01004      fprintf(stderr,"IMREG: making empty images\n") ;
01005 
01006 
01007    datum = DSET_BRICK_TYPE(new_dset,0) ;
01008 
01009    INIT_IMARR(ims_in) ;
01010    for( ii=0 ; ii < ntime ; ii++ ){
01011       im = mri_new_vol_empty( nx , ny , 1 , datum ) ;
01012       ADDTO_IMARR(ims_in,im) ;
01013    }
01014 
01015    imbase = mri_new_vol_empty( nx , ny , 1 , base_datum ) ;
01016 
01017    dxar  = (float *) malloc( sizeof(float) * ntime ) ;
01018    dyar  = (float *) malloc( sizeof(float) * ntime ) ;
01019    phiar = (float *) malloc( sizeof(float) * ntime ) ;
01020 
01021    /*** 3) Get pointers to sub-bricks in old, base, and new datasets ***/
01022 
01023    if (opt->debug)
01024      fprintf(stderr,"IMREG: getting input brick pointers\n") ;
01025 
01026 
01027    switch( datum ){  /* pointer type depends on input datum type */
01028       case MRI_byte:
01029          bptr  = (byte **) malloc( sizeof(byte *) * ntime ) ;
01030          bout  = (byte **) malloc( sizeof(byte *) * ntime ) ;
01031          for( ii=0 ; ii < ntime ; ii++ ){
01032             bptr[ii]  = (byte *) DSET_ARRAY(old_dset,ii) ;
01033             bout[ii]  = (byte *) DSET_ARRAY(new_dset,ii) ;
01034          }
01035       break ;
01036 
01037       case MRI_short:
01038          sptr  = (short **) malloc( sizeof(short *) * ntime ) ;
01039          sout  = (short **) malloc( sizeof(short *) * ntime ) ;
01040          for( ii=0 ; ii < ntime ; ii++ ){
01041             sptr[ii]  = (short *) DSET_ARRAY(old_dset,ii) ;
01042             sout[ii]  = (short *) DSET_ARRAY(new_dset,ii) ;
01043          }
01044 
01045          if (opt->debug)
01046            fprintf(stderr,"IMREG: sptr[0] = %p  sout[0] = %p\n",
01047                    sptr[0],sout[0]) ;
01048 
01049       break ;
01050 
01051       case MRI_float:
01052          fptr  = (float **) malloc( sizeof(float *) * ntime ) ;
01053          fout  = (float **) malloc( sizeof(float *) * ntime ) ;
01054          for( ii=0 ; ii < ntime ; ii++ ){
01055             fptr[ii]  = (float *) DSET_ARRAY(old_dset,ii) ;
01056             fout[ii]  = (float *) DSET_ARRAY(new_dset,ii) ;
01057          }
01058       break ;
01059    }
01060 
01061    switch( base_datum ){  /* pointer type depends on base datum type */
01062       case MRI_byte:
01063         bbase = (byte *) DSET_ARRAY(base_dset,base) ; 
01064       break ;
01065 
01066       case MRI_short:
01067         sbase = (short *) DSET_ARRAY(base_dset,base) ;
01068         if (opt->debug)
01069           fprintf(stderr,"IMREG: sbase[%d] = %p \n", base, sbase) ;
01070       break ;
01071 
01072       case MRI_float:
01073         fbase = (float *) DSET_ARRAY(base_dset,base) ;
01074       break ;
01075    }
01076 
01077    /*** 4) Loop over slices ***/
01078 
01079    for( kk=0 ; kk < nz ; kk++ ){
01080 
01081       /*** 4a) Setup ims_in images to point to input slices ***/
01082 
01083      if (opt->debug)
01084        fprintf(stderr,"IMREG: slice %d -- setup input images\n",kk) ;
01085 
01086 
01087       for( ii=0 ; ii < ntime ; ii++ ){
01088          im = IMARR_SUBIMAGE(ims_in,ii) ;
01089          switch( datum ){
01090             case MRI_byte:  
01091               mri_fix_data_pointer( bptr[ii] + kk*npix, im ) ; break ;
01092             case MRI_short: 
01093               mri_fix_data_pointer( sptr[ii] + kk*npix, im ) ; break ;
01094             case MRI_float: 
01095               mri_fix_data_pointer( fptr[ii] + kk*npix, im ) ; break ;
01096          }
01097       }
01098 
01099       /*** 4b) Setup imbase to point to base image ***/
01100 
01101       if (opt->debug)
01102         fprintf(stderr,"IMREG: slice %d -- setup base image\n",kk) ;
01103 
01104 
01105       switch( base_datum ){
01106          case MRI_byte:  
01107            mri_fix_data_pointer( bbase + kk*npix, imbase ) ; break ;
01108          case MRI_short: 
01109            mri_fix_data_pointer( sbase + kk*npix, imbase ) ; break ;
01110          case MRI_float: 
01111            mri_fix_data_pointer( fbase + kk*npix, imbase ) ; break ;
01112       }
01113 
01114       /*** 4c) Register this slice at all times ***/
01115 
01116       if (opt->debug)
01117         fprintf(stderr,"IMREG: slice %d -- register\n",kk) ;
01118 
01119 
01120       ims_out = mri_align_dfspace( imbase , NULL , ims_in ,
01121                                    ALIGN_REGISTER_CODE , dxar,dyar,phiar ) ;
01122 
01123       if( ims_out == NULL )
01124         fprintf(stderr,"IMREG: mri_align_dfspace return NULL\n") ;
01125 
01126       
01127       /*----- Store the registration parameters -----*/
01128       if (opt->dprefix != NULL)
01129         for (ii = 0;  ii < ntime;  ii++)
01130           {
01131             it = ii*nz + z_to_t[kk];
01132             if (opt->dmm)
01133               {
01134                 state_history[it].elts[1] = dxar[ii] * dx;
01135                 state_history[it].elts[2] = dyar[ii] * dy;
01136               }
01137             else
01138               {
01139                 state_history[it].elts[1] = dxar[ii];
01140                 state_history[it].elts[2] = dyar[ii];
01141               }
01142 
01143             state_history[it].elts[3] = (180.0/PI)*phiar[ii];
01144           }
01145       
01146 
01147       /*** 4d) Put the output back in on top of the input;
01148                note that the output is always in MRI_float format ***/
01149 
01150       if (opt->debug)
01151         fprintf(stderr,"IMREG: slice %d -- put output back into dataset\n",kk);
01152 
01153 
01154       for( ii=0 ; ii < ntime ; ii++ ){
01155          switch( datum ){
01156            case MRI_byte:
01157               im = mri_to_mri( MRI_byte , IMARR_SUBIMAGE(ims_out,ii) ) ;
01158               memcpy( bout[ii] + kk*npix , MRI_BYTE_PTR(im) , 
01159                       sizeof(byte)*npix ) ;
01160               mri_free(im) ;
01161            break ;
01162 
01163            case MRI_short:
01164 
01165              if (opt->debug)
01166                if( ii==0 )
01167                  fprintf(stderr,"IMREG: conversion to short at ii=%d\n",ii) ;
01168 
01169               im = mri_to_mri( MRI_short , IMARR_SUBIMAGE(ims_out,ii) ) ;
01170 
01171               if (opt->debug)
01172                 if( ii==0 )
01173                   fprintf(stderr,"IMREG: copying to %p from %p\n",
01174                           sout[ii] + kk*npix,MRI_SHORT_PTR(im)) ;
01175 
01176 
01177               memcpy( sout[ii] + kk*npix , MRI_SHORT_PTR(im) , 
01178                       sizeof(short)*npix ) ;
01179 
01180               if (opt->debug)
01181                 if( ii==0 )
01182                   fprintf(stderr,"IMREG: freeing\n") ;
01183 
01184 
01185               mri_free(im) ;
01186            break ;
01187 
01188            case MRI_float:
01189               im = IMARR_SUBIMAGE(ims_out,ii) ;
01190               memcpy( fout[ii] + kk*npix , MRI_FLOAT_PTR(im) , 
01191                       sizeof(float)*npix ) ;
01192            break ;
01193          }
01194       }
01195 
01196       /*** 4e) Destroy the output images ***/
01197 
01198       if (opt->debug)
01199         fprintf(stderr,"IMREG: destroying aligned output\n") ;
01200 
01201 
01202       DESTROY_IMARR( ims_out ) ;
01203    }
01204 
01205    /*** 5) Destroy the empty images and other workspaces ***/
01206 
01207    if (opt->debug)
01208      fprintf(stderr,"IMREG: destroy workspaces\n") ;
01209 
01210 
01211    mri_clear_data_pointer(imbase) ; mri_free(imbase) ;
01212    for( ii=0 ; ii < ntime ; ii++ ){
01213       im = IMARR_SUBIMAGE(ims_in,ii) ;
01214       mri_clear_data_pointer(im) ;
01215    }
01216    DESTROY_IMARR(ims_in) ;
01217    FREE_WORKSPACE ;
01218 
01219    /*------------- write out the new dataset ------------*/
01220 
01221    if (opt->debug)
01222      fprintf(stderr,"IMREG: write new dataset to disk\n") ;
01223 
01224 
01225   /*----- Record history of dataset -----*/
01226    tross_Copy_History( old_dset , new_dset ) ;
01227    if( commandline != NULL )
01228      tross_Append_History( new_dset , commandline ) ;
01229 
01230 
01231   THD_load_statistics( new_dset ) ;
01232   THD_write_3dim_dataset( NULL,NULL , new_dset , True ) ;
01233   
01234 
01235   /*----- evaluate results -----*/
01236   if (opt->rprefix != NULL)
01237     eval_registration (opt, old_dset, new_dset, base_dset, base, 
01238                        old_rms_array, new_rms_array);
01239 
01240 
01241   /*----- deallocate memory -----*/   
01242   THD_delete_3dim_dataset( old_dset , False ) ; old_dset = NULL ;
01243   THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
01244   if (opt->base_filename != NULL)
01245     THD_delete_3dim_dataset( base_dset , False ) ; base_dset = NULL ;
01246     
01247 
01248 
01249    return NULL ;  /* null string returned means all was OK */
01250 }
01251 
01252 
01253 /*---------------------------------------------------------------------------*/
01254 /*
01255   Routine to write RMS error arrays. 
01256 */
01257 
01258 
01259 void output_rms_arrays
01260 (
01261   IR_options * option_data,     /* user input options */
01262   THD_3dim_dataset * dset,      /* pointer to input 3d+time dataset */
01263   float * old_rms_array,        /* volume RMS error for input dataset */
01264   float * new_rms_array         /* volume RMS error for registered dataset */
01265 )
01266 
01267 {
01268   int it;                                  /* time index */
01269   int ts_length;                           /* length of time series data */  
01270   char filename[MAX_NAME_LENGTH];          /* name of output file */
01271   FILE * old_rms_file , * new_rms_file;    /* file for volume RMS error */
01272 
01273 
01274   /*----- Initialize local variables -----*/
01275   ts_length = DSET_NUM_TIMES(dset);
01276 
01277 
01278   /*----- Open output files -----*/
01279   strcpy (filename, option_data->rprefix);
01280   strcat (filename, ".oldrms");
01281   old_rms_file = fopen (filename, "w");
01282   strcpy (filename, option_data->rprefix);
01283   strcat (filename, ".newrms");
01284   new_rms_file = fopen (filename, "w");
01285 
01286 
01287   /*----- Write volume RMS error data -----*/
01288   for (it = 0;  it < ts_length;  it++)
01289     {
01290       fprintf (old_rms_file, "%d  %f\n" , it, old_rms_array[it]);
01291       fprintf (new_rms_file, "%d  %f\n" , it, new_rms_array[it]);
01292     }
01293 
01294 
01295   /*----- Close output files -----*/
01296   fclose (old_rms_file);
01297   fclose (new_rms_file);
01298 
01299 }
01300 
01301 
01302 /*---------------------------------------------------------------------------*/
01303 /*
01304   Get the time corresponding to this particular slice.
01305 */
01306  
01307 float get_time  (int ivolume, int iz,  THD_3dim_dataset * dset)
01308 
01309 {
01310   float time;
01311   float z;
01312 
01313 
01314   z = iz * dset->taxis->dz_sl + dset->taxis->zorg_sl;
01315   time = THD_timeof (ivolume, z, dset->taxis);
01316   
01317   if (dset->taxis->units_type == UNITS_MSEC_TYPE)
01318     time /= 1000.0;
01319 
01320   return (time);
01321 }
01322 
01323 
01324 /*---------------------------------------------------------------------------*/
01325 /*
01326   Routine to write the registration parameter time series. 
01327 */
01328 
01329 
01330 void output_state_history
01331 (
01332   IR_options * option_data,     /* user input options */
01333   THD_3dim_dataset * dset,
01334   vector * state_history
01335 )
01336 
01337 {
01338   int iv;                           /* vector index */
01339   int ts_length;                    /* length of time series */
01340   int num_slices;                   /* number of slices in each volume */
01341   int num_vectors;                  /* total number of state vectors */
01342   int ivolume;                      /* volume index */
01343   int iz;                           /* z slice index */
01344   float t;                          /* time for current vector */
01345   char filename[MAX_NAME_LENGTH];   /* string for output file name */
01346 
01347   FILE * dx_file, * dy_file, * psi_file;   /* output files */
01348 
01349 
01350   /*----- Initialize local variables -----*/
01351   num_slices = dset->taxis->nsl;
01352   ts_length = DSET_NUM_TIMES(dset);
01353 
01354   if ( num_slices <= 0 )            /* 06 Oct 2003 [rickr] */
01355       num_slices = dset->daxes->nzz;
01356 
01357 
01358   /*----- Calculate total number of state vectors -----*/
01359   num_vectors = ts_length * num_slices;
01360 
01361 
01362   /*----- Open output files -----*/
01363   strcpy (filename, option_data->dprefix);
01364   strcat (filename, ".dx");
01365   dx_file = fopen (filename, "w");
01366 
01367   strcpy (filename, option_data->dprefix);
01368   strcat (filename, ".dy");
01369   dy_file = fopen (filename, "w");
01370 
01371   strcpy (filename, option_data->dprefix);
01372   strcat (filename, ".psi");
01373   psi_file = fopen (filename, "w");
01374 
01375   
01376   /*----- Write the registration parameters -----*/
01377   for (iv = 0;  iv < num_vectors;  iv++)
01378     {
01379       ivolume = iv / num_slices;
01380       iz = t_to_z[iv % num_slices];
01381       t = get_time (ivolume, iz, dset);
01382       fprintf (dx_file,  "%f   %f\n", t, state_history[iv].elts[1]);
01383       fprintf (dy_file,  "%f   %f\n", t, state_history[iv].elts[2]);
01384       fprintf (psi_file, "%f   %f\n", t, state_history[iv].elts[3]);
01385     }
01386 
01387 
01388   /*----- Close output files -----*/
01389   fclose (dx_file);
01390   fclose (dy_file);
01391   fclose (psi_file);
01392 }
01393 
01394 
01395 /*---------------------------------------------------------------------------*/
01396 /*
01397   Output results.
01398 */
01399 
01400 void output_results  
01401 (
01402   IR_options * option_data,     /* user input options */
01403   vector * state_history,       /* time series of state vectors */
01404   float * old_rms_array,        /* volume RMS error for input dataset */
01405   float * new_rms_array         /* volume RMS error for registered dataset */
01406 )
01407 
01408 {
01409   THD_3dim_dataset * dset;
01410 
01411 
01412   read_dataset (option_data->input_filename, &dset);
01413 
01414 
01415   /*----- Write the time series of state parameters -----*/
01416   if (option_data->dprefix != NULL)
01417     output_state_history (option_data, dset, state_history);
01418 
01419 
01420   /*----- Write user specified auxiliary time series data -----*/
01421   if (option_data->rprefix != NULL)
01422     output_rms_arrays (option_data, dset, old_rms_array, new_rms_array);
01423 
01424 
01425   THD_delete_3dim_dataset (dset, False);   dset = NULL;
01426 
01427 }
01428 
01429 
01430 /*---------------------------------------------------------------------------*/
01431 /*
01432   Terminate the program.
01433 */
01434 
01435 void terminate_program  
01436 (
01437   IR_options ** option_data,
01438   vector ** state_history,       /* time series of state vectors */
01439   float ** old_rms_array,        /* volume RMS error for input dataset */
01440   float ** new_rms_array         /* volume RMS error for registered dataset */
01441 )
01442 
01443 {
01444   THD_3dim_dataset * dset = NULL;   /* pointer to input 3d+time dataset */
01445   int num_slices;                   /* number of slices in each volume */
01446   int ts_length;                    /* length of time series */
01447   int num_vectors;                  /* total number of state vectors */
01448   int i;                            /* index */
01449 
01450 
01451   /*----- Initialize local variables -----*/
01452   read_dataset ((*option_data)->input_filename, &dset);
01453   num_slices = dset->taxis->nsl;
01454 
01455   if ( num_slices <= 0 )            /* 06 Oct 2003 [rickr] */
01456       num_slices = dset->daxes->nzz;
01457 
01458   ts_length = DSET_NUM_TIMES(dset);
01459   num_vectors = ts_length * num_slices;
01460   THD_delete_3dim_dataset (dset, False);   dset = NULL;
01461 
01462 
01463   /*----- Release memory -----*/
01464   free (*option_data);     *option_data = NULL;
01465   free (t_to_z);           t_to_z = NULL;
01466   free (z_to_t);           z_to_t = NULL;
01467 
01468   
01469   if (*old_rms_array != NULL)
01470     { free (*old_rms_array);   *old_rms_array = NULL; }
01471   if (*new_rms_array != NULL)
01472     { free (*new_rms_array);   *new_rms_array = NULL; }
01473 
01474 
01475   /*----- Deallocate memory for array of state vectors -----*/
01476   if (*state_history != NULL)
01477     {
01478       for (i = 0;  i < num_vectors;  i++)
01479         vector_destroy (&((*state_history)[i]));
01480       free (*state_history);   *state_history = NULL;
01481     }
01482 
01483 }
01484 
01485 
01486 /*---------------------------------------------------------------------------*/
01487 
01488 int main
01489 (
01490   int argc,                    /* number of input arguments */
01491   char ** argv                 /* array of input arguments */ 
01492 )
01493 
01494 {
01495   IR_options * option_data = NULL;     /* user input options */
01496   char * chptr;                        /* error message from processing */
01497   vector * state_history = NULL;       /* time series of state vectors */
01498   float * old_rms_array = NULL;        /* original data volume RMS error */
01499   float * new_rms_array = NULL;        /* registered data volume RMS error */
01500 
01501 
01502    
01503   /*----- Identify software -----*/
01504   printf ("\n\n");
01505   printf ("Program:          %s \n", PROGRAM_NAME);
01506   printf ("Initial Release:  %s \n", PROGRAM_INITIAL);
01507   printf ("Latest Revision:  %s \n", PROGRAM_LATEST);
01508   printf ("\n");
01509 
01510   
01511   /*----- Program initialization -----*/
01512   initialize_program (argc, argv, &option_data, &state_history, 
01513                       &old_rms_array, &new_rms_array);
01514 
01515 
01516   /*----- Register all slices in the dataset -----*/
01517   chptr = IMREG_main (option_data, state_history,
01518                       old_rms_array, new_rms_array);
01519 
01520 
01521   /*----- Check for processing errors -----*/
01522   if (chptr != NULL)   
01523     {
01524       printf ("%s \n\n", chptr);
01525       exit(1);
01526     }
01527 
01528 
01529   /*----- Output the results -----*/
01530   output_results (option_data, state_history, 
01531                   old_rms_array, new_rms_array);
01532 
01533 
01534   /*----- Terminate the program -----*/
01535   terminate_program (&option_data, &state_history,
01536                      &old_rms_array, &new_rms_array);
01537  
01538 }
01539 
01540 
01541