00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2001, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006 
00007 /*---------------------------------------------------------------------------*/
00008 /*
00009   Program to resample a 'data parent' dataset to the grid defined by an
00010   'anat parent' dataset.
00011 
00012   File:    adwarp.c
00013   Author:  B. Douglas Ward
00014   Date:    02 April 1999
00015 
00016   Mod:     Added changes for incorporating History notes.
00017   Date:    10 September 1999
00018 
00019   Mod:     Added changes for proper byte ordering on output.
00020   Date:    23 Nov 1999 - RW Cox
00021 
00022   Mod:     Added -force option, and some checks related to it.
00023   Date:    13 Dec 1999 - RW Cox
00024 
00025   Mod:     Added 'data parent' sub-brick selector.  Also, allow operator to
00026            specify separate resampling modes for the functional sub-bricks
00027            and the threshold sub-bricks.
00028   Date:    25 February 2000
00029 
00030   Mod:     Added call to AFNI_logger.
00031   Date:    15 August 2001
00032 
00033 */
00034 
00035 /*---------------------------------------------------------------------------*/
00036 
00037 #define PROGRAM_NAME "adwarp.c"                      /* name of this program */
00038 #define PROGRAM_AUTHOR "R. W. Cox and B. D. Ward"          /* program author */
00039 #define PROGRAM_INITIAL "02 April 1999"   /* date of initial program release */
00040 #define PROGRAM_LATEST "15 August 2001"     /* date of last program revision */
00041 
00042 /*---------------------------------------------------------------------------*/
00043 
00044 #include "afni_warp.h"
00045 
00046 #define MAIN
00047 
00048 /*---------------------------------------------------------------------------*/
00049 
00050 void AFNI_copy_statistics( THD_3dim_dataset *dsold , THD_3dim_dataset *dsnew )
00051 {
00052    int ibr , nvold , nvnew ;
00053    THD_statistics *stold , *stnew ;
00054 
00055 ENTRY("AFNI_copy_statistics") ;
00056 
00057    if( !ISVALID_3DIM_DATASET(dsold) || !ISVALID_3DIM_DATASET(dsnew) ) EXRETURN ;
00058 
00059    nvold = dsold->dblk->nvals ;
00060    nvnew = dsnew->dblk->nvals ;
00061    stold = dsold->stats ;
00062    stnew = dsnew->stats ;
00063    if( !ISVALID_STATISTIC(stold) ) EXRETURN ;
00064 
00065    if( stnew == NULL ){
00066       dsnew->stats  = stnew = myXtNew( THD_statistics ) ;
00067       stnew->type   = STATISTICS_TYPE ;
00068       stnew->nbstat = nvnew ;
00069       stnew->bstat  = (THD_brick_stats *)
00070                         XtMalloc( sizeof(THD_brick_stats) * nvnew ) ;
00071       ADDTO_KILL(dsnew->kl,stnew) ;
00072       stnew->parent = (XtPointer) dsnew ;
00073    } else {
00074       stnew->nbstat = nvnew ;
00075       stnew->bstat  = (THD_brick_stats *)
00076                         XtRealloc( (char *) stnew->bstat ,
00077                                    sizeof(THD_brick_stats) * nvnew ) ;
00078    }
00079 
00080    for( ibr=0 ; ibr < nvnew ; ibr++ ){
00081       if( ibr < nvold )
00082          stnew->bstat[ibr] = stold->bstat[ibr] ;
00083       else
00084          INVALIDATE_BSTAT(stnew->bstat[ibr]) ;
00085    }
00086 
00087    EXRETURN ;
00088 }
00089 
00090 /*---------------------------------------------------------------------------*/
00091 
00092 #define PARENTIZE(ds,par) \
00093   if( ISVALID_3DIM_DATASET((ds)) ) (ds)->parent = (XtPointer) (par)
00094 
00095 /** macro to test a malloc-ed pointer for validity **/
00096 #define MTEST(ptr) \
00097   if((ptr)==NULL) \
00098   ( AW_error ("Cannot allocate memory") )
00099 
00100 /*---------------------------------------------------------------------------*/
00101 
00102 typedef struct adwarp_options
00103 {
00104   THD_3dim_dataset *aset;       /* anat parent dataset */
00105   THD_3dim_dataset *dset;       /* data parent dataset */
00106 
00107   char *prefix;                 /* prefix for output file */
00108 
00109   float dxyz;                   /* grid spacing in output dataset (mm) */
00110 
00111   int anat_resam_mode;          /* anatomical dataset resampling mode */
00112   int thr_resam_mode;           /* threshold sub-brick resampling mode */
00113   int func_resam_mode;          /* functional sub-brick resampling mode */
00114 
00115   int verbose ;                 /* RWCox: 06 Apr 1999 */
00116   int force   ;                 /* RWCox: 13 Dec 1999 */
00117 
00118 } adwarp_options;
00119 
00120 
00121 /*---------------------------------------------------------------------------*/
00122 /*
00123    Print error message and stop.
00124 */
00125 
00126 void AW_error (char *message)
00127 {
00128   fprintf (stderr, "%s Error: %s \n", PROGRAM_NAME, message);
00129   exit(1);
00130 }
00131 
00132 
00133 /*---------------------------------------------------------------------------*/
00134 /*
00135   Routine to display adwarp help menu.
00136 */
00137 
00138 void display_help_menu()
00139 {
00140    int ii ;
00141 
00142 
00143    printf
00144      ("Usage: adwarp [options]\n"
00145       "Resamples a 'data parent' dataset to the grid defined by an\n"
00146       "'anat parent' dataset.  The anat parent dataset must contain\n"
00147       "in its .HEAD file the coordinate transformation (warp) needed\n"
00148       "to bring the data parent dataset to the output grid.  This\n"
00149       "program provides a batch implementation of the interactive\n"
00150       "AFNI 'Write' buttons, one dataset at a time.\n"
00151       "\n"
00152       "  Example: adwarp -apar anat+tlrc -dpar func+orig\n"
00153       "\n"
00154       "  This will create dataset func+tlrc (.HEAD and .BRIK).\n"
00155       "\n"
00156       "Options (so to speak):\n"
00157       "----------------------\n"
00158       "-apar aset  = Set the anat parent dataset to 'aset'.  This\n"
00159       "                is a nonoptional option (must be present).\n"
00160       "\n"
00161       "-dpar dset  = Set the data parent dataset to 'dset'.  This\n"
00162       "                is a nonoptional option (must be present).\n"
00163       "              Note: dset may contain a sub-brick selector,\n"
00164       "              e.g.,  -dpar 'dset+orig[2,5,7]'             \n"
00165       "\n"
00166       "-prefix ppp = Set the prefix for the output dataset to 'ppp'.\n"
00167       "                The default is the prefix of 'dset'.\n"
00168       "\n"
00169       "-dxyz ddd   = Set the grid spacing in the output datset to\n"
00170       "                'ddd' mm.  The default is 1 mm.\n"
00171       "\n"
00172       "-verbose    = Print out progress reports.\n"
00173       "-force      = Write out result even if it means deleting\n"
00174       "                an existing dataset.  The default is not\n"
00175       "                to overwrite.\n"
00176       "\n"
00177       "-resam rrr  = Set resampling mode to 'rrr' for all sub-bricks\n"
00178       "                     --- OR ---                              \n"
00179       "-thr   rrr  = Set resampling mode to 'rrr' for threshold sub-bricks\n"
00180       "-func  rrr  = Set resampling mode to 'rrr' for functional sub-bricks\n"
00181       "\n"
00182       "The resampling mode 'rrr' must be one of the following:\n"
00183       ) ;
00184 
00185    for( ii=0 ; ii <= LAST_RESAM_TYPE ; ii++ )
00186       printf("                 %s = %s\n", RESAM_shortstr[ii],RESAM_typestr[ii] ) ;
00187 
00188 
00189    printf
00190      (
00191       "\n"
00192       "NOTE:  The default resampling mode is %s for all sub-bricks. \n" ,
00193       RESAM_shortstr[1]
00194       ) ;
00195 
00196    exit(0) ;
00197 }
00198 
00199 
00200 /*---------------------------------------------------------------------------*/
00201 /*
00202   Routine to initialize the input options.
00203 */
00204 
00205 
00206 void initialize_options (adwarp_options *option_data)
00207 {
00208   option_data->aset = NULL;
00209   option_data->dset = NULL;
00210   option_data->prefix = NULL;
00211 
00212   option_data->dxyz = 1.0;
00213 
00214   option_data->anat_resam_mode = 1;
00215   option_data->thr_resam_mode  = 1;
00216   option_data->func_resam_mode = 1;
00217 
00218   option_data->verbose = 0 ;  /* 06 Apr 1999 */
00219   option_data->force   = 0 ;  /* 13 Dec 1999 */
00220 }
00221 
00222 
00223 /*---------------------------------------------------------------------------*/
00224 /*
00225   Routine to get user specified input options.
00226 */
00227 
00228 void get_options
00229 (
00230   int argc,                        /* number of input arguments */
00231   char ** argv,                    /* array of input arguments */
00232   adwarp_options *option_data     /* adwarp program options */
00233 )
00234 
00235 {
00236   int nopt = 1;                     /* input option argument counter */
00237   int ival;                         /* integer input */
00238   float fval;                       /* float input */
00239   int ii;                           /* index */
00240   char message[THD_MAX_NAME];       /* error message */
00241 
00242 
00243   /*----- does user request help menu? -----*/
00244   if (argc < 2 || strncmp(argv[1], "-help", 5) == 0)  display_help_menu();
00245 
00246   AFNI_logger(PROGRAM_NAME,argc,argv) ;
00247 
00248 
00249   /*----- initialize the input options -----*/
00250   initialize_options (option_data);
00251 
00252 
00253   /*----- main loop over input options -----*/
00254   while (nopt < argc )
00255     {
00256 
00257       /*-----   -apar aset   -----*/
00258       if (strncmp(argv[nopt], "-apar", 5) == 0)
00259         {
00260           nopt++;
00261           if (nopt >= argc)  AW_error ("need argument after -apar ");
00262           option_data->aset = THD_open_one_dataset (argv[nopt]);
00263           if (! ISVALID_3DIM_DATASET(option_data->aset))
00264             AW_error ("Cannot read anat parent dataset.\n") ;
00265           nopt++;
00266           continue;
00267         }
00268 
00269       /*-----   -dpar dset   -----*/
00270       if (strncmp(argv[nopt], "-dpar", 5) == 0)
00271         {
00272           nopt++;
00273           if (nopt >= argc)  AW_error ("need argument after -dpar ");
00274           option_data->dset = THD_open_dataset (argv[nopt]);
00275           if (! ISVALID_3DIM_DATASET(option_data->dset))
00276             AW_error ("Cannot read data parent dataset.\n") ;
00277           nopt++;
00278           continue;
00279         }
00280 
00281       /*-----   -prefix ppp   -----*/
00282       if (strncmp(argv[nopt], "-prefix", 7) == 0)
00283         {
00284           nopt++;
00285           if (nopt >= argc)  AW_error ("need argument after -prefix ");
00286           option_data->prefix = AFMALL(char,sizeof(char)*THD_MAX_PREFIX);
00287           MTEST (option_data->prefix);
00288           MCW_strncpy (option_data->prefix, argv[nopt], THD_MAX_PREFIX);
00289 
00290      if( strstr(option_data->prefix,".nii") != NULL ){  /* 06 Apr 2005 */
00291        fprintf(stderr,"** You can't use adwarp to create a .nii file!\n") ;
00292        exit(1) ;
00293      }
00294 
00295           nopt++;
00296           continue;
00297         }
00298 
00299       /*-----  -verbose  -----*/
00300       if( strncmp(argv[nopt],"-verbose",5) == 0 ){  /* 06 Apr 1999 */
00301          option_data->verbose = 1 ;
00302          nopt++ ; continue ;
00303       }
00304 
00305       /*-----  -force  -----*/
00306       if( strncmp(argv[nopt],"-force",5) == 0 ){  /* 06 Apr 1999 */
00307          option_data->force = 1 ;
00308          nopt++ ; continue ;
00309       }
00310 
00311       /*-----   -dxyz ddd   -----*/
00312       if (strncmp(argv[nopt], "-dxyz", 5) == 0)
00313         {
00314           nopt++;
00315           if (nopt >= argc)  AW_error ("need argument after -dxyz ");
00316           sscanf (argv[nopt], "%f", &fval);
00317           if (fval <= 0.0)  AW_error ("Illegal argument after -dxyz ");
00318           option_data->dxyz = fval;
00319           nopt++;
00320           continue;
00321         }
00322 
00323       /*-----   -resam rrr  -----*/
00324       if (strncmp(argv[nopt], "-resam", 6) == 0)
00325         {
00326           nopt++;
00327           if (nopt >= argc)  AW_error ("need argument after -resam ");
00328           ival = -1;
00329           for( ii=0 ; ii <= LAST_RESAM_TYPE ; ii++ )
00330             if (strncmp(argv[nopt], RESAM_shortstr[ii], 2) == 0)
00331               ival = ii;
00332           if ((ival < 0) || (ival > LAST_RESAM_TYPE))
00333             AW_error ("illegal argument after -resam ");
00334           option_data->anat_resam_mode = ival;
00335           option_data->thr_resam_mode  = ival;
00336           option_data->func_resam_mode = ival;
00337           nopt++;
00338           continue;
00339         }
00340 
00341 
00342       /*-----   -thr rrr  -----*/
00343       if (strncmp(argv[nopt], "-thr", 4) == 0)
00344         {
00345           nopt++;
00346           if (nopt >= argc)  AW_error ("need argument after -thr ");
00347           ival = -1;
00348           for( ii=0 ; ii <= LAST_RESAM_TYPE ; ii++ )
00349             if (strncmp(argv[nopt], RESAM_shortstr[ii], 2) == 0)
00350               ival = ii;
00351           if ((ival < 0) || (ival > LAST_RESAM_TYPE))
00352             AW_error ("illegal argument after -thr ");
00353           option_data->thr_resam_mode  = ival;
00354           nopt++;
00355           continue;
00356         }
00357 
00358 
00359       /*-----   -func rrr  -----*/
00360       if (strncmp(argv[nopt], "-func", 5) == 0)
00361         {
00362           nopt++;
00363           if (nopt >= argc)  AW_error ("need argument after -func ");
00364           ival = -1;
00365           for( ii=0 ; ii <= LAST_RESAM_TYPE ; ii++ )
00366             if (strncmp(argv[nopt], RESAM_shortstr[ii], 2) == 0)
00367               ival = ii;
00368           if ((ival < 0) || (ival > LAST_RESAM_TYPE))
00369             AW_error ("illegal argument after -func ");
00370           option_data->func_resam_mode  = ival;
00371           nopt++;
00372           continue;
00373         }
00374 
00375 
00376       /*----- unknown command -----*/
00377       sprintf(message,"Unrecognized command line option: %s\n", argv[nopt]);
00378       AW_error (message);
00379 
00380     }
00381 
00382 
00383   /*----- Check for required inputs -----*/
00384   if (option_data->aset == NULL)
00385     AW_error ("Must specify anat parent dataset");
00386   if (option_data->dset == NULL)
00387     AW_error ("Must specify data parent dataset");
00388 
00389   /*-- 13 Dec 1999: check if datasets have the same view --*/
00390 
00391   if( option_data->aset->view_type == option_data->dset->view_type ){
00392       if( option_data->force ){
00393          if( option_data->prefix == NULL ){
00394             fprintf(stderr,
00395                     "** Error: -apar & -dpar are in same +view!\n"
00396                     "          This is illegal without the -prefix option!\n") ;
00397             exit(1) ;
00398          } else {
00399             fprintf(stderr,
00400                     "++ Warning: -apar & -dpar are in same +view!\n") ;
00401          }
00402       } else {
00403          fprintf(stderr,
00404                  "** Error: -apar & -dpar are in same +view!\n"
00405                  "          If this is OK, use -force and -prefix options.\n" ) ;
00406          exit(1) ;
00407       }
00408   }
00409 
00410 
00411 }
00412 
00413 
00414 /*---------------------------------------------------------------------------*/
00415 /*
00416    Make a warped dataset whose grid corresponds to the anat_parent and
00417    whose data comes from the data_parent.
00418    Note that the assumption is made that data_parent and the warp parent
00419    of the anat_parent are both in the same coordinate system (up to the
00420    to_dicomm transformation of their dataxes structs).
00421 
00422    This routine is adapted from AFNI_follower_dataset.
00423 */
00424 
00425 THD_3dim_dataset * adwarp_follower_dataset
00426 (
00427   adwarp_options   *option_data,    /* adwarp program options */
00428   THD_3dim_dataset *anat_parent,    /* dataset containing warp information */
00429   THD_3dim_dataset *data_parent     /* dataset to be warped */
00430 )
00431 {
00432   THD_3dim_dataset *new_dset ;
00433   int ii ;
00434 
00435 ENTRY("adwarp_follower_dataset") ;
00436 
00437 /* sanity checks */
00438 
00439   if( ! ISVALID_3DIM_DATASET(anat_parent) ||
00440       ! ISVALID_3DIM_DATASET(data_parent)   ) RETURN(NULL) ;
00441 
00442   /* make new dataset, copying appropriate fields from its various parents */
00443 
00444   new_dset = myXtNew( THD_3dim_dataset ) ; INIT_KILL( new_dset->kl ) ;
00445 
00446   new_dset->type      = data_parent->type;        /* same data type */
00447   new_dset->func_type = data_parent->func_type;
00448   new_dset->view_type = anat_parent->view_type;   /* but different view type */
00449 
00450   new_dset->anat_parent = anat_parent;            /* what else makes sense? */
00451 
00452   new_dset->tagset = NULL ;  /* Oct 1998 */
00453 
00454   MCW_strncpy( new_dset->anat_parent_name ,
00455                anat_parent->self_name , THD_MAX_NAME ) ;
00456 
00457   new_dset->anat_parent_idcode = anat_parent->idcode ;
00458 
00459    /* 11/09/94 addition: the data_parent may itself be a warp;
00460        in this case, we want the true warp parent to be the original data */
00461 
00462   new_dset->warp_parent =  (data_parent->warp_parent != NULL)
00463                          ? (data_parent->warp_parent) : (data_parent) ;
00464 
00465   MCW_strncpy( new_dset->warp_parent_name ,
00466                new_dset->warp_parent->self_name , THD_MAX_NAME ) ;
00467 
00468   new_dset->warp_parent_idcode = new_dset->warp_parent->idcode ;
00469 
00470   new_dset->idcode = MCW_new_idcode() ;
00471 
00472   /* make the actual warp from the warp_parent to this dataset */
00473 
00474   new_dset->vox_warp       = NULL ;
00475   new_dset->warp           = myXtNew( THD_warp ) ;
00476   *(new_dset->warp)        = IDENTITY_WARP ;  /* start with (Dicom) identity */
00477 
00478   new_dset->self_warp      = NULL ;           /* 26 Aug 2002 */
00479 
00480   /* follow the links backward from desired view to original view */
00481 
00482   AFNI_concatenate_warp( new_dset->warp , anat_parent->warp ) ;
00483   AFNI_concatenate_warp( new_dset->warp , data_parent->warp ) ;
00484 
00485   /* reset the bounds in the new warp to be the same as in the anat_parent */
00486 
00487   if( ISVALID_WARP(anat_parent->warp) &&
00488       anat_parent->warp->type == new_dset->warp->type ){
00489 
00490     switch( anat_parent->warp->type ){
00491         
00492     case WARP_AFFINE_TYPE:
00493       COPY_LMAP_BOUNDS( new_dset->warp->rig_bod.warp ,
00494                         anat_parent->warp->rig_bod.warp ) ;
00495       break ;
00496 
00497     case WARP_TALAIRACH_12_TYPE:
00498       for( ii=0 ; ii < 12 ; ii++ )
00499         COPY_LMAP_BOUNDS( new_dset->warp->tal_12.warp[ii] ,
00500                           anat_parent->warp->tal_12.warp[ii] ) ;
00501       break ;
00502     }
00503   }
00504 
00505   /* make up some names for this new dataset */
00506 
00507   MCW_strncpy( new_dset->self_name  ,
00508                new_dset->warp_parent->self_name , THD_MAX_NAME ) ;
00509   ii = strlen( new_dset->self_name ) ;
00510   new_dset->self_name[ii++] = '@' ;
00511   MCW_strncpy( &(new_dset->self_name[ii]) ,
00512                VIEW_typestr[new_dset->view_type] , THD_MAX_NAME-ii ) ;
00513 
00514   MCW_strncpy( new_dset->label1 , data_parent->label1 , THD_MAX_LABEL ) ;
00515   MCW_strncpy( new_dset->label2 , data_parent->label2 , THD_MAX_LABEL ) ;
00516 
00517   /* set the axes for this new dataset
00518      (same as anatomy parent, since that's the meaning of this routine) */
00519 
00520   new_dset->daxes         = myXtNew( THD_dataxes ) ;  /* copy data axes of */
00521   *(new_dset->daxes)      = *(anat_parent->daxes) ; /* anatomy parent */
00522 
00523   new_dset->wod_daxes     = NULL ;
00524   new_dset->wod_flag      = True ;
00525 
00526   /* 06 Aug 1996: added ability to use 3D+t datasets here */
00527 
00528   if( DSET_NUM_TIMES(data_parent) < 2 ){
00529     new_dset->taxis = NULL ;
00530   } else {
00531     new_dset->taxis  = myXtNew( THD_timeaxis ) ;  /* new */
00532     *(new_dset->taxis) = *(data_parent->taxis) ;  /* copy insides */
00533 
00534     new_dset->taxis->nsl     = 0 ;                      /* no slice stuff */
00535     new_dset->taxis->toff_sl = NULL ;
00536     new_dset->taxis->zorg_sl = 0.0 ;
00537     new_dset->taxis->dz_sl   = 0.0 ;
00538   }
00539 
00540   /* create a datablock and diskptr, in case the data is ever
00541      filled into memory (instead of wod) and written to disk */
00542 
00543   new_dset->dblk = myXtNew( THD_datablock ) ; INIT_KILL( new_dset->dblk->kl ) ;
00544 
00545   new_dset->dblk->type        = DATABLOCK_TYPE ;
00546   new_dset->dblk->nvals       = data_parent->dblk->nvals ;
00547   new_dset->dblk->malloc_type = DATABLOCK_MEM_UNDEFINED ;
00548   new_dset->dblk->natr        = new_dset->dblk->natr_alloc  = 0 ;
00549   new_dset->dblk->atr         = NULL ;
00550   new_dset->dblk->parent      = (XtPointer) new_dset ;
00551 
00552   if( data_parent->dblk->brick_lab == NULL ){
00553     THD_init_datablock_labels( new_dset->dblk ) ; /* 30 Nov 1997 */
00554   } else {
00555     THD_copy_datablock_auxdata( data_parent->dblk , new_dset->dblk ) ;
00556   }
00557 
00558   DSET_unlock(new_dset) ;  /* Feb 1998 */
00559 
00560   new_dset->dblk->diskptr               = myXtNew( THD_diskptr ) ;
00561   new_dset->dblk->diskptr->type         = DISKPTR_TYPE ;
00562   new_dset->dblk->diskptr->nvals        = data_parent->dblk->nvals ;
00563   new_dset->dblk->diskptr->rank         = 3 ;
00564   new_dset->dblk->diskptr->storage_mode = STORAGE_UNDEFINED ;
00565   new_dset->dblk->diskptr->byte_order   = THD_get_write_order() ;
00566                                                             /* 25 April 1998 */
00567   new_dset->dblk->diskptr->dimsizes[0]  = new_dset->daxes->nxx ;
00568   new_dset->dblk->diskptr->dimsizes[1]  = new_dset->daxes->nyy ;
00569   new_dset->dblk->diskptr->dimsizes[2]  = new_dset->daxes->nzz ;
00570 
00571   new_dset->dblk->brick_fac   = NULL ;  /* initialized below */
00572   new_dset->dblk->brick_bytes = NULL ;
00573   new_dset->dblk->brick       = NULL ;
00574   THD_init_datablock_brick( new_dset->dblk , -1 , data_parent->dblk ) ;
00575 
00576   new_dset->dblk->master_nvals = 0 ;     /* 11 Jan 1999 */
00577   new_dset->dblk->master_ival  = NULL ;
00578   new_dset->dblk->master_bytes = NULL ;
00579 
00580   /* create the names for storage on disk (if ever)
00581      -- note we put it in the same directory as the data_parent */
00582 
00583   if (option_data->prefix == NULL)
00584     THD_init_diskptr_names (new_dset->dblk->diskptr,
00585                             data_parent->dblk->diskptr->directory_name, NULL,
00586                             data_parent->dblk->diskptr->prefix,
00587                             new_dset->view_type, True);
00588   else
00589     THD_init_diskptr_names (new_dset->dblk->diskptr,
00590                             data_parent->dblk->diskptr->directory_name, NULL,
00591                             option_data->prefix,
00592                             new_dset->view_type, True);
00593 
00594 
00595   ADDTO_KILL( new_dset->dblk->kl , new_dset->dblk->diskptr ) ;
00596 
00597   /* oh yeah, set the new_dset kill list,
00598      copy statistics if available, and NULL out any unused stuff */
00599 
00600   ADDTO_KILL( new_dset->kl , new_dset->warp ) ;
00601   ADDTO_KILL( new_dset->kl , new_dset->daxes ) ;
00602   ADDTO_KILL( new_dset->kl , new_dset->dblk ) ;
00603 
00604   new_dset->stats = NULL ;
00605   AFNI_copy_statistics( data_parent , new_dset ) ;
00606 
00607   INIT_STAT_AUX( new_dset , MAX_STAT_AUX , data_parent->stat_aux ) ;
00608 
00609   new_dset->markers     = NULL ;  /* no markers */
00610   new_dset->death_mark  = 0 ;     /* don't kill me! */
00611 #ifdef ALLOW_DATASET_VLIST
00612   new_dset->pts         = NULL ;
00613 #endif
00614 
00615   PARENTIZE(new_dset,data_parent->parent) ;
00616 
00617   new_dset->tcat_list   = NULL ;  /* 03 Aug 2004 */
00618   new_dset->tcat_num    = 0 ;
00619   new_dset->tcat_len    = NULL ;
00620 
00621   RETURN(new_dset) ;
00622 }
00623 
00624 
00625 /*---------------------------------------------------------------------------*/
00626 /*
00627    (re)compute and (re)write a dataset to disk, with the indicated
00628    geometric parameters.  Note that the filenames for the save should
00629    already be initialized in the dset->dblk->diskptr.  Also, note
00630    that the dataset's "permanent" dataxes will be remade to fit the
00631    new geometry.
00632 
00633    This routine is adapted from AFNI_refashion_dataset.
00634 */
00635 
00636 
00637 Boolean adwarp_refashion_dataset
00638 (
00639   adwarp_options   *option_data,  /* adwarp program options */
00640   THD_3dim_dataset *dset,         /* new (output) dataset */
00641   THD_dataxes      *daxes         /* new dataset axes */
00642 )
00643 {
00644   THD_datablock *dblk  = dset->dblk ;
00645   THD_diskptr   *dkptr = dset->dblk->diskptr ;
00646   Boolean good ;
00647   int npix , nx,ny,nz,nv , kk , ival , code , nzv , dsiz , isfunc , cmode ;
00648   MRI_IMAGE *im ;
00649   void *imar ;
00650   FILE *far ;
00651   float brfac_save ;
00652   int resam_mode;
00653 
00654   int native_order , save_order ;  /* 23 Nov 1999 */
00655 
00656 ENTRY("adwarp_refashion_dataset") ;
00657 
00658   /* set up for warp-on-demand */
00659 
00660   dset->wod_daxes         = myXtNew(THD_dataxes) ; /* 02 Nov 1996 */
00661   dset->wod_daxes->type   = DATAXES_TYPE ;       /* 02 Nov 1996 */
00662   dset->vox_warp          = myXtNew(THD_warp) ;    /* 02 Nov 1996 */
00663 
00664   dset->self_warp         = NULL ;                 /* 26 Aug 2002 */
00665 
00666   *(dset->wod_daxes)      = *daxes ;            /* copy insides of daxes */
00667   dset->wod_flag          = True ;              /* mark for warp-on-demand */
00668   dset->vox_warp->type    = ILLEGAL_TYPE ;      /* mark for recomputation */
00669 
00670   /* copy the new geometric information into various places */
00671 
00672   *(dset->daxes)     = *(daxes) ;               /* Make daxes permanent */
00673   dkptr->dimsizes[0] = dset->daxes->nxx ;       /* Will cause trouble */
00674   dkptr->dimsizes[1] = dset->daxes->nyy ;       /* if diskptr and     */
00675   dkptr->dimsizes[2] = dset->daxes->nzz ;       /* daxes don't match! */
00676 
00677   /* write the header out */
00678 
00679   good = THD_write_3dim_dataset( NULL,NULL , dset , False ) ;
00680   if( !good ){
00681     fprintf(stderr,"\a\n*** cannot write dataset header ***\n") ;
00682 
00683     RETURN(False) ;
00684   }
00685   STATUS("wrote output header file") ;
00686 
00687   /* purge the datablock that now exists,
00688      then delete the file on disk that now exists (if any) */
00689 
00690   DSET_unlock( dset ) ; /* Feb 1998 */
00691   PURGE_DSET( dset ) ;
00692   COMPRESS_unlink(dkptr->brick_name) ;
00693 
00694   /* refashion its brick data structure,
00695      which requires first saving in a temporary
00696      array the datum type for each sub-brick    */
00697 
00698   { int ibr ; int *typ ;
00699     typ = (int *) XtMalloc( sizeof(int) * dblk->nvals ) ;
00700     for( ibr=0 ; ibr < dblk->nvals ; ibr++ )
00701       typ[ibr] = DBLK_BRICK_TYPE(dblk,ibr) ;
00702     THD_init_datablock_brick( dblk , dblk->nvals , typ ) ;
00703     myXtFree( typ ) ;
00704   }
00705 
00706   if( dblk->total_bytes > 500*1024*1024 ){
00707     int mb = (int)(dblk->total_bytes/(1024*1024)) ;
00708     fprintf(stderr,"++ WARNING: output filesize will be %d Mbytes!\n"
00709                    "++ SUGGEST: increasing voxel size to save disk space!\n",mb) ;
00710   }
00711 
00712   dkptr->storage_mode = STORAGE_UNDEFINED ;       /* just for now */
00713   dblk->malloc_type   = DATABLOCK_MEM_UNDEFINED ;
00714 
00715   /*--- open the output file
00716     (N.B.: much of the following code is from THD_write_datablock) ---*/
00717 
00718   /*-- create directory if necessary --*/
00719 
00720   if( ! THD_is_directory(dkptr->directory_name) ){
00721     kk = mkdir( dkptr->directory_name , THD_MKDIR_MODE ) ;
00722     if( kk != 0 ){
00723       fprintf(stderr,
00724             "\a\n*** cannot mkdir new directory: %s\n",dkptr->directory_name) ;
00725 
00726       RETURN(False) ;
00727     }
00728     STATUS("created subdirectory") ;
00729   }
00730 
00731   /*-- open output file --*/
00732 
00733   if( option_data->verbose )
00734     fprintf(stderr,"++ Opening output file %s\n",dkptr->brick_name) ;
00735 
00736   cmode = THD_get_write_compression() ;
00737   far = COMPRESS_fopen_write( dkptr->brick_name , cmode ) ;
00738   if( far == NULL ){
00739     fprintf(stderr,
00740             "\a\n*** cannot open output file %s\n",dkptr->brick_name) ;
00741 
00742     RETURN(False) ;
00743   }
00744   STATUS("created output brick file") ;
00745 
00746   /*--------- now, create each slice and write it out ----------*/
00747 
00748   nx = dset->daxes->nxx ;
00749   ny = dset->daxes->nyy ;  npix = nx*ny ;
00750   nz = dset->daxes->nzz ;
00751   nv = dkptr->nvals     ;  nzv  = nz*nv ;
00752 
00753   isfunc = ISFUNC(dset) ;  /* 09 Dec 1997 */
00754 
00755   if( ! isfunc )
00756     resam_mode = option_data->anat_resam_mode ;
00757 
00758    native_order = mri_short_order() ;                           /* 23 Nov 1999 */
00759    save_order   = (dkptr->byte_order > 0) ? dkptr->byte_order
00760                                           : THD_get_write_order() ;
00761 
00762   for( ival=0 ; ival < nv ; ival++ ){  /* for each sub-brick */
00763 
00764     if( option_data->verbose )
00765        fprintf(stderr,"++ Start sub-brick %d",ival) ;
00766 
00767     dsiz = mri_datum_size( DSET_BRICK_TYPE(dset,ival) ) ;
00768 
00769     /** force return of unscaled slices for output **/
00770 
00771     brfac_save                   = DBLK_BRICK_FACTOR(dblk,ival) ;
00772     DBLK_BRICK_FACTOR(dblk,ival) = 0.0 ;
00773 
00774     if( isfunc )
00775       resam_mode = (DSET_BRICK_STATCODE(dset,ival) > 0)  /* 09 Dec 1997 */
00776       ? option_data->thr_resam_mode
00777       : option_data->func_resam_mode ;
00778 
00779     for( kk=0 ; kk < nz ; kk++ ){  /* for each slice */
00780 
00781       im = AFNI_dataset_slice( dset , 3 , kk , ival , resam_mode ) ;
00782 STATUS("have new image") ;
00783 
00784       if( option_data->verbose && kk%7==0 ) fprintf(stderr,".");
00785 
00786       if( im == NULL ){
00787         fprintf(stderr,"\a\n*** failure to compute dataset slice %d\n",kk) ;
00788         COMPRESS_fclose(far) ;
00789         COMPRESS_unlink( dkptr->brick_name ) ;
00790 
00791         RETURN(False) ;
00792       }
00793 
00794 #ifdef AFNI_DEBUG
00795 { char str[256] ;
00796   sprintf(str,"writing slice %d: type=%s nx=%d ny=%d\n",
00797           kk,MRI_TYPE_NAME(im) , im->nx,im->ny ) ;
00798   STATUS(str) ; }
00799 #endif
00800 
00801         imar = mri_data_pointer(im) ;
00802          if( save_order != native_order ){                   /* 23 Nov 1999 */
00803             switch( im->kind ){
00804                case MRI_short:   mri_swap2(  npix,imar) ; break ;
00805                case MRI_float:
00806                case MRI_int:     mri_swap4(  npix,imar) ; break ;
00807                case MRI_complex: mri_swap4(2*npix,imar) ; break ;
00808             }
00809          }
00810         code = fwrite( imar , dsiz , npix , far ) ;
00811         mri_free(im) ;
00812         
00813         if( code != npix ){
00814           fprintf(stderr,
00815                   "\a\n*** failure to write dataset slice %d (is disk full?)\n",kk) ;
00816           COMPRESS_fclose(far) ;
00817           COMPRESS_unlink( dkptr->brick_name ) ;
00818         
00819           RETURN(False) ;
00820         }
00821         
00822     } /* end of loop over kk (z-direction) */
00823 
00824     if( option_data->verbose ) fprintf(stderr,"\n");
00825 
00826     /* restore the correct scaling of this sub-brick */
00827 
00828     DBLK_BRICK_FACTOR(dblk,ival) = brfac_save ;
00829 
00830   } /* end of loop over iv (nvals direction) */
00831   STATUS("all slices written") ;
00832 
00833   /*--------------------- done!!! ---------------------*/
00834 
00835   COMPRESS_fclose(far) ;
00836   STATUS("output file closed") ;
00837 
00838   /*--- do a little surgery on the dataset's storage flags ---*/
00839 
00840   dkptr->storage_mode = STORAGE_BY_BRICK ;
00841 #if MMAP_THRESHOLD > 0
00842   dblk->malloc_type   = (dblk->total_bytes > MMAP_THRESHOLD)
00843                       ? DATABLOCK_MEM_MMAP : DATABLOCK_MEM_MALLOC ;
00844 
00845   if( cmode >= 0 ) dblk->malloc_type = DATABLOCK_MEM_MALLOC ;
00846   DBLK_mmapfix(dblk) ;  /* 28 Mar 2005 */
00847 #else
00848   dblk->malloc_type   = DATABLOCK_MEM_MALLOC ;
00849 #endif
00850 
00851   /*--- recompute the statistics and rewrite the header to hold them ---*/
00852 
00853   STATUS("recomputing statistics") ;
00854 
00855   if( option_data->verbose ) fprintf(stderr,"++ Computing statistics\n");
00856   THD_load_statistics( dset ) ;
00857 
00858   STATUS("rewriting header") ;
00859 
00860   (void) THD_write_3dim_dataset( NULL,NULL , dset , False ) ;
00861 
00862   STATUS("purging datablock") ;
00863 
00864   PURGE_DSET( dset ) ;
00865 
00866   myXtFree(dset->wod_daxes) ; myXtFree(dset->vox_warp) ;  /* 02 Nov 1996 */
00867 
00868   RETURN(True) ;
00869 }
00870 
00871 
00872 /*---------------------------------------------------------------------------*/
00873 
00874 int main( int argc , char *argv[] )
00875 {
00876   adwarp_options   *option_data;           /* adwarp program options */
00877   THD_3dim_dataset *new_dset = NULL;       /* new (output) dataset */
00878   THD_dataxes new_daxes;                   /* new dataset axes */
00879 
00880   mainENTRY("adwarp main") ; machdep() ;
00881 
00882   /*----- Identify software -----*/
00883   printf ("\n\n");
00884   printf ("Program:          %s \n", PROGRAM_NAME);
00885   printf ("Author:           %s \n", PROGRAM_AUTHOR);
00886   printf ("Initial Release:  %s \n", PROGRAM_INITIAL);
00887   printf ("Latest Revision:  %s \n", PROGRAM_LATEST);
00888   printf ("\n");
00889 
00890   /*----- Allocate memory -----*/
00891   option_data = (adwarp_options *) malloc (sizeof(adwarp_options));
00892 
00893   /*----- Get user inputs -----*/
00894   get_options (argc, argv, option_data);
00895 
00896   /*----- Create empty shell of dataset from anat parent and data parent ----*/
00897   new_dset = adwarp_follower_dataset (option_data, option_data->aset,
00898                                       option_data->dset);
00899 
00900   /*--- 13 Dec 1999: check if output dataset already exists on disk ---*/
00901 
00902   if( THD_is_file(DSET_HEADNAME(new_dset)) ||
00903       THD_is_file(DSET_BRIKNAME(new_dset))    ){
00904 
00905       if( option_data->force ){
00906          fprintf(stderr,
00907                  "++ Warning: overwriting dataset %s and %s\n",
00908                  DSET_HEADNAME(new_dset), DSET_BRIKNAME(new_dset) ) ;
00909       } else {
00910          fprintf(stderr,
00911                  "** Error: can't overwrite dataset %s and %s\n"
00912                  "          unless you use the -force option!\n" ,
00913                  DSET_HEADNAME(new_dset), DSET_BRIKNAME(new_dset) ) ;
00914          exit(1) ;
00915       }
00916    }
00917 
00918   /*----- Record history of dataset -----*/
00919   tross_Copy_History( option_data->dset , new_dset ) ;
00920   tross_Make_History( PROGRAM_NAME , argc , argv , new_dset ) ;
00921 
00922   /*----- Allow for resampling to a new voxel size -----*/
00923   new_daxes.type = DATAXES_TYPE;
00924   THD_edit_dataxes (option_data->dxyz, option_data->aset->daxes, &new_daxes);
00925 
00926 
00927   /*----- Fill in the dataset and write out to disk -----*/
00928   adwarp_refashion_dataset (option_data, new_dset, &new_daxes);
00929 
00930   exit(0) ;
00931 }