自学教程：C++ DSET_load函数代码示例

/*---------------------------------------------------------------------  23 Feb 2012: Return the absolute value of the difference between                two volumes, divided by the number of voxels               and the number of sub-bricks. Voxels that are zero               in both sets are not counted.               Comparisons are done after conversion of data to double    return = -1.0 ERROR           =  0.0 Exactly the same-----------------------------------------------------------------------*/double THD_diff_vol_vals(THD_3dim_dataset *d1, THD_3dim_dataset *d2, int scl) {   double dd=0.0, denom=0.0;   int i=0, k=0;   double *a1=NULL, *a2=NULL;   MRI_IMAGE *b1 = NULL , *b2 = NULL;      ENTRY("THD_diff_vol_vals");      if (!d1 && !d2) RETURN(dd);   if (!d1 || !d2) RETURN(-1.0);   if (!EQUIV_GRIDS(d1,d2)) RETURN(-1.0);   if (DSET_NVALS(d1) != DSET_NVALS(d2)) RETURN(-1.0);     DSET_mallocize(d1) ; DSET_load(d1) ;   DSET_mallocize(d2) ; DSET_load(d2) ;   dd = 0.0; denom = 0;   for (i=0; i<DSET_NVALS(d1); ++i) {      b1 = THD_extract_double_brick(i, d1);      b2 = THD_extract_double_brick(i, d2);      a1 = MRI_DOUBLE_PTR(b1);      a2 = MRI_DOUBLE_PTR(b2);      for (k=0; k<DSET_NVOX(d1); ++k) {         dd += ABS(a1[k]-a2[k]);         if (a1[k]!=0.0 || a2[k]!=0.0) ++denom;      }      mri_clear_data_pointer(b1); mri_free(b1) ;      mri_clear_data_pointer(b2); mri_free(b2) ;   }   if (scl && denom>0.0) dd /= denom;      RETURN(dd);   }  

/* * for each input dataset name *    open (check dims, etc.) *    dilate (zeropad, make binary, dilate, unpad, apply) *    fill list of bytemask datasets * * also, count total volumes */int process_input_dsets(param_t * params){   THD_3dim_dataset * dset, * dfirst=NULL;   int                iset, nxyz;   ENTRY("process_input_dsets");   if( !params ) ERROR_exit("NULL inputs to PID");   if( params->ndsets <= 0 ) {      ERROR_message("process_input_dsets: no input datasets");      RETURN(1);   }   /* allocate space for dsets array */   params->dsets = (THD_3dim_dataset **)malloc(params->ndsets*                                               sizeof(THD_3dim_dataset*));   if( !params->dsets ) ERROR_exit("failed to allocate dset pointers");   if( params->verb ) INFO_message("processing %d input datasets...",                                   params->ndsets);      /* warn user of dilations */   if(params->verb && params->ndsets) {      int pad = needed_padding(&params->IND);      INFO_message("padding all datasets by %d (for dilations)", pad);   }   /* process the datasets */   nxyz = 0;   for( iset=0; iset < params->ndsets; iset++ ) {      /* open and verify dataset */      dset = THD_open_dataset(params->inputs[iset]);      if( !dset ) ERROR_exit("failed to open mask dataset '%s'",                             params->inputs[iset]);      DSET_load(dset);  CHECK_LOAD_ERROR(dset);      if( params->verb>1 ) INFO_message("loaded dset %s, with %d volumes",                                        DSET_PREFIX(dset), DSET_NVALS(dset));      if( nxyz == 0 ) { /* make an empty copy of the first dataset */         nxyz = DSET_NVOX(dset);         dfirst = EDIT_empty_copy(dset);      }      /* check for consistency in voxels and grid */      if( DSET_NVOX(dset) != nxyz ) ERROR_exit("nvoxel mis-match");      if( ! EQUIV_GRIDS(dset, dfirst) )         WARNING_message("grid from dset %s does not match that of dset %s",                         DSET_PREFIX(dset), DSET_PREFIX(dfirst));      /* apply dilations to all volumes, returning bytemask datasets */      params->dsets[iset] = apply_dilations(dset, &params->IND,1,params->verb);      if( ! params->dsets[iset] ) RETURN(1);   }    DSET_delete(dfirst); /* and nuke */   RETURN(0);}

int main( int argc , char *argv[] ){   THD_3dim_dataset *dset ; int aa,ll ; char *cpt ; float val ;   if( argc < 2 ){     printf("Usage: 3dSatCheck dataset [...]/n"            "/n"            "Prints the 'raw' initial transient (saturation) check/n"            "value for each dataset on the command line.  Round this/n"            "number to the nearest integer to get an estimate of/n"            "how many non-saturated time points start a dataset./n"           ) ;     exit(0) ;   }   for( aa=1 ; aa < argc ; aa++ ){     dset = THD_open_dataset( argv[aa] ) ; if( !ISVALID_DSET(dset) ) continue ;     if( DSET_NVALS(dset) < 9 ) continue ;     DSET_load(dset) ; if( !DSET_LOADED(dset) ) continue ;     val = THD_saturation_check( dset , NULL , 0,0 ) ;     ll = strlen(argv[aa]) ;     cpt = (ll <= 50) ? argv[aa] : argv[aa]+(ll-50) ;     INFO_message("%-50.50s = %.3f",cpt,val) ;     DSET_delete(dset) ;   }   exit(0) ;}

/*----------------------------------------------------------------------****  Main routine for this plugin (will be called from AFNI).****----------------------------------------------------------------------*/char * HEMISUB_main( PLUGIN_interface * plint ){    THD_3dim_dataset * dset, * new_dset;    MCW_idcode       * idc;    hemi_s             hs = { 0 };    char             * new_prefix;    char             * ret_string = NULL;    char             * tag;    if ( plint == NULL )	return  "------------------------/n"		"HEMISUB_main: NULL input/n"		"------------------------/n";    PLUTO_next_option( plint );    idc  = PLUTO_get_idcode( plint );    dset = PLUTO_find_dset( idc );    if( dset == NULL )	return "-------------------------------/n"	       "HEMISUB_main: bad input dataset/n"	       "-------------------------------";    DSET_load( dset );    PLUTO_next_option( plint );    new_prefix = PLUTO_get_string( plint );    if ( ! PLUTO_prefix_ok( new_prefix ) )	return  "------------------------/n"		"HEMISUB_main: bad prefix/n"		"------------------------/n";    if ( ( new_dset = PLUTO_copy_dset( dset, new_prefix ) ) == NULL )	return  "------------------------------------------/n"		"HEMISUB_main: failed to copy input dataset/n"		"------------------------------------------/n";    tag = PLUTO_get_optiontag( plint );    if ( tag && ! strcmp( tag, "Thresh Type" ) )    {	tag = PLUTO_get_string( plint );	if ( tag != NULL )	    hs.thresh_type = PLUTO_string_index( tag, NUM_T_OPTS, thresh_opts );    }    if ( ret_string = process_data( new_dset, &hs ) )	return  ret_string;    if ( PLUTO_add_dset( plint, new_dset, DSET_ACTION_MAKE_CURRENT ) )    {	THD_delete_3dim_dataset( new_dset, False );	return  "---------------------------------------/n"		"HEMISUB_main: failed to add new dataset/n"		"---------------------------------------/n";    }    return NULL;}

/* - moved from 3dhistog.c               18 Dec 2012 [rickr] */int THD_slow_minmax_dset(THD_3dim_dataset *dset, float *dmin, float *dmax,                         int iv_bot, int iv_top){   int iv_fim;   float fimfac;   float vbot , vtop, temp_fbot=1.0, temp_ftop=0.0 ;   DSET_load(dset);   for( iv_fim=iv_bot ; iv_fim <= iv_top ; iv_fim++ ){     /* minimum and maximum for sub-brick */     vbot = mri_min( DSET_BRICK(dset,iv_fim) ) ;     vtop = mri_max( DSET_BRICK(dset,iv_fim) ) ;     fimfac = DSET_BRICK_FACTOR(dset,iv_fim) ;     if (fimfac == 0.0)  fimfac = 1.0;     vbot *= fimfac ; vtop *= fimfac ;     /* update global min and max */     if ( temp_fbot > temp_ftop ) { /* first time, just copy */        temp_fbot = vbot;        temp_ftop = vtop;     } else {        if( vbot < temp_fbot ) temp_fbot = vbot;        if( vtop > temp_ftop ) temp_ftop = vtop;     }   }   *dmin = temp_fbot;   *dmax = temp_ftop;   return(0);}

int main( int argc , char *argv[] ){   THD_3dim_dataset *dset ;   MRI_IMAGE *im ;   int iarg=1 ;   if( strcmp(argv[1],"-nper") == 0 ){      nper = strtol( argv[2] , NULL , 10 ) ;      iarg = 3 ;   }   for( ; iarg < argc ; iarg++ ){     printf("======= dataset %s  nper=%d ========/n",argv[iarg],nper) ;     dset = THD_open_dataset(argv[iarg]) ;     if( dset == NULL ) continue ;     DSET_load(dset) ;     im = DSET_BRICK(dset,0) ;     im->dx = DSET_DX(dset) ;     im->dy = DSET_DY(dset) ;     im->dz = DSET_DZ(dset) ;     (void) THD_orient_guess( im ) ;     printf( "Data Axes Orientation:/n"             "  first  (x) = %s/n"             "  second (y) = %s/n"             "  third  (z) = %s/n" ,        ORIENT_typestr[dset->daxes->xxorient] ,        ORIENT_typestr[dset->daxes->yyorient] ,        ORIENT_typestr[dset->daxes->zzorient]  ) ;     DSET_delete(dset) ;   }   exit(0) ;}

MRI_IMAGE * THD_median_brick( THD_3dim_dataset *dset ){   int nvox , nvals , ii ;   MRI_IMAGE *tsim , *medim ;   float *medar ;   float *tsar ;  /* 05 Nov 2001 */ENTRY("THD_median_brick") ;   if( !ISVALID_DSET(dset) ) RETURN(NULL) ;   DSET_load(dset) ;   if( !DSET_LOADED(dset) ) RETURN(NULL) ;   nvals = DSET_NVALS(dset) ;   tsim  = DSET_BRICK(dset,0) ;   if( nvals == 1 ){     medim = mri_scale_to_float( DSET_BRICK_FACTOR(dset,0), tsim ) ;     RETURN(medim) ;   }   medim = mri_new_conforming( tsim , MRI_float ) ;   medar = MRI_FLOAT_PTR(medim) ;   nvox  = DSET_NVOX(dset) ;   tsar = (float *) calloc( sizeof(float),nvals+1 ) ; /* 05 Nov 2001 */   for( ii=0 ; ii < nvox ; ii++ ){     THD_extract_array( ii , dset , 0 , tsar ) ;     /* 05 Nov 2001 */     medar[ii] = qmed_float( nvals , tsar ) ;   }   free(tsar) ; RETURN(medim) ;}

MRI_IMARR * THD_medmad_bricks( THD_3dim_dataset *dset ){   int nvox , nvals , ii ;   MRI_IMAGE *tsim , *madim, *medim ;   float             *madar, *medar ;   MRI_IMARR *imar ;   float *tsar ;ENTRY("THD_medmad_bricks") ;   if( !ISVALID_DSET(dset) ) RETURN(NULL) ;   nvals = DSET_NVALS(dset) ; if( nvals == 1 ) RETURN(NULL) ;   DSET_load(dset) ;  if( !DSET_LOADED(dset) ) RETURN(NULL) ;   tsim  = DSET_BRICK(dset,0) ;   madim = mri_new_conforming( tsim , MRI_float ) ;   madar = MRI_FLOAT_PTR(madim) ;   medim = mri_new_conforming( tsim , MRI_float ) ;   medar = MRI_FLOAT_PTR(medim) ;   nvox  = DSET_NVOX(dset) ;   tsar = (float *) calloc( sizeof(float),nvals+1 ) ;   for( ii=0 ; ii < nvox ; ii++ ){     THD_extract_array( ii , dset , 0 , tsar ) ;     qmedmad_float( nvals , tsar , medar+ii , madar+ii ) ;   }   free(tsar) ;   INIT_IMARR(imar) ; ADDTO_IMARR(imar,medim) ; ADDTO_IMARR(imar,madim) ;   RETURN(imar) ;}

MRI_IMAGE * THD_rms_brick( THD_3dim_dataset *dset ){   int nvox , nvals , ii , jj ;   MRI_IMAGE *tsim , *medim ;   float *medar , sum,fac ;   float *tsar ;ENTRY("THD_rms_brick") ;   if( !ISVALID_DSET(dset) ) RETURN(NULL) ;   DSET_load(dset) ;   if( !DSET_LOADED(dset) ) RETURN(NULL) ;   nvals = DSET_NVALS(dset)   ; fac = 1.0 / nvals ;   tsim  = DSET_BRICK(dset,0) ;   if( nvals == 1 ){     medim = mri_scale_to_float( DSET_BRICK_FACTOR(dset,0), tsim ) ;     RETURN(medim) ;   }   medim = mri_new_conforming( tsim , MRI_float ) ;   medar = MRI_FLOAT_PTR(medim) ;   nvox  = DSET_NVOX(dset) ;   tsar = (float *) calloc( sizeof(float),nvals+1 ) ;   for( ii=0 ; ii < nvox ; ii++ ){     THD_extract_array( ii , dset , 0 , tsar ) ;     for( sum=0.0,jj=0 ; jj < nvals ; jj++ ) sum += tsar[jj]*tsar[jj] ;     medar[ii] = sqrtf(fac * sum) ;   }   free(tsar) ; RETURN(medim) ;}

double ENTROPY_dataset( THD_3dim_dataset *dset ){   if( !ISVALID_DSET(dset) ) return(0.0) ;   DSET_load(dset) ;   if( !DSET_LOADED(dset) ) return(0.0) ;   return ENTROPY_datablock( dset->dblk ) ;}

/* Calculates the average value for each voxel in dset across all timepoints.   PRE:   dset is a valid 3D+T dataset with at least 1 timepoint;   ignore is the number of timepoints to ignore at the beginning of dset;   0 <= ignore < number of timepoints in dset;   POST:   return value is NULL on error,   else, return value is an array of floats with the same dimensions as the   subbricks of dset, containing the voxel value averages;   Note: The caller is responsible for free()ing the returned block of memory   when done.   Note2: The complex datatype is not supported, and any such bricks will   result in an error (NULL return value).*/double * RIC_CalcVoxelMeans(const THD_3dim_dataset * dset, int ignore) {    double * avg;       /* The voxel averages to be returned */    float scalefactor; /* Current dset brick scaling factor */    int ival, nvals;   /* Current, number of dset timepoints */    int ivox, nvoxs;   /* Current, number of dset brick voxels */    /* Quick check of arguments */    if (!ISVALID_3DIM_DATASET(dset) || DSET_NVALS(dset) < 1 ||	ignore < 0 || ignore >= DSET_NVALS(dset)) {	return NULL;    }    /* Initialize */    DSET_load(dset);    nvals = DSET_NVALS(dset);    nvoxs = dset->daxes->nxx * dset->daxes->nyy * dset->daxes->nzz;    avg = malloc(sizeof(double) * nvoxs);    if (avg == NULL) {	return NULL;    }    /* Calculate the voxel averages; treat matrices as 1-D arrays */    /* Zero the voxel sums */    for (ivox = 0; ivox < nvoxs; ivox += 1) {	avg[ivox] = 0.0;    }    /* Sum each voxel across time (and hope there are not too many points) */    for (ival = ignore; ival < nvals; ival += 1) {	scalefactor = DSET_BRICK_FACTOR(dset, ival);	switch (DSET_BRICK_TYPE(dset, ival)) {	case MRI_short:	    RIC_CALCVOXELMEANS__DO_VOXSUM(short);	    break;	case MRI_byte:	    RIC_CALCVOXELMEANS__DO_VOXSUM(byte);	    break;	case MRI_float:	    RIC_CALCVOXELMEANS__DO_VOXSUM(float);	    break;	default: /* Unsupported datatype */	    free(avg);	    return NULL;	}    }    /* Divide by number of timepoints to get average */    nvals -= ignore;  /* We do not average over the ignored timepoints */    for (ivox = 0; ivox < nvoxs; ivox += 1) {	avg[ivox] /= nvals;    }    return avg;}

THD_3dim_dataset *Seg_load_dset_eng( char *set_name, char *view ) {   static char FuncName[]={"Seg_load_dset_eng"};   THD_3dim_dataset *dset=NULL, *sdset=NULL;   int i=0;   byte make_cp=0;   int verb=0;   char sprefix[THD_MAX_PREFIX+10], *stmp=NULL;      SUMA_ENTRY;      dset = THD_open_dataset( set_name );   if( !ISVALID_DSET(dset) ){     fprintf(stderr,"**ERROR: can't open dataset %s/n",set_name) ;     SUMA_RETURN(NULL);   }      DSET_mallocize(dset)   ; DSET_load(dset);      for (i=0; i<DSET_NVALS(dset); ++i) {      if (DSET_BRICK_TYPE(dset,i) != MRI_short) {         if (verb) INFO_message("Sub-brick %d in %s not of type short./n"                       "Creating new short copy of dset ",                        i, DSET_PREFIX(dset));         make_cp=1; break;      }   }      if (make_cp) {      if (!SUMA_ShortizeDset(&dset, -1.0)) {         SUMA_S_Err("**ERROR: Failed to shortize");         SUMA_RETURN(NULL);      }   }      if (DSET_IS_MASTERED(dset)) {      if (verb) INFO_message("Dset is mastered, making copy...");      stmp = SUMA_ModifyName(set_name, "append", ".cp", NULL);      sdset = dset;      dset = EDIT_full_copy(sdset, stmp);      free(stmp); DSET_delete(sdset); sdset = NULL;     }            if (view) {      if (view) {               if (!strstr(view,"orig"))             EDIT_dset_items(dset,ADN_view_type, VIEW_ORIGINAL_TYPE, ADN_none);          else  if (!strstr(view,"acpc"))             EDIT_dset_items(dset,ADN_view_type, VIEW_ACPCALIGNED_TYPE, ADN_none);         else  if (!strstr(view,"tlrc"))             EDIT_dset_items(dset ,ADN_view_type, VIEW_TALAIRACH_TYPE, ADN_none);         else SUMA_S_Errv("View of %s is rubbish", view);      }   }      SUMA_RETURN(dset);}

void get_options(  int argc,                        /* number of input arguments */  char ** argv                     /* array of input arguments */ ){  int nopt = 1;                     /* input option argument counter */  int ival, index;                  /* integer input */  float fval;                       /* float input */  char message[MAX_STRING_LENGTH];  /* error message */  /*----- does user request help menu? -----*/  if (argc < 2 || strncmp(argv[1], "-help", 5) == 0)  display_help_menu();       /*----- main loop over input options -----*/  while (nopt < argc )    {      /*-----   -anat filename   -----*/      if (strncmp(argv[nopt], "-anat", 5) == 0)	{	  nopt++;	  if (nopt >= argc)  estPDF_error ("need argument after -anat ");	  anat_filename = malloc (sizeof(char) * MAX_STRING_LENGTH);	  MTEST (anat_filename);	  strcpy (anat_filename, argv[nopt]);	  anat = THD_open_one_dataset (anat_filename);	  if (!ISVALID_3DIM_DATASET (anat))	    {	      sprintf (message, "Can't open dataset: %s/n", anat_filename); 	      estPDF_error (message); 	    } 	  DSET_load(anat); CHECK_LOAD_ERROR(anat);	  nopt++;	  continue;	}            /*----- unknown command -----*/      sprintf(message,"Unrecognized command line option: %s/n", argv[nopt]);      estPDF_error (message);          }  }

int fill_mask(options_t * opts){   THD_3dim_dataset * mset;   int nvox;ENTRY("fill_mask");   if( opts->automask ) {      if( opts->verb ) INFO_message("creating automask...");      opts->mask = THD_automask(opts->inset);      if( ! opts->mask ) {         ERROR_message("failed to apply -automask");         RETURN(1);      }      RETURN(0);   }   if( opts->mask_name ) {      if( opts->verb )         INFO_message("reading mask dset from %s...", opts->mask_name);      mset = THD_open_dataset( opts->mask_name );      if( ! mset ) ERROR_exit("cannot open mask dset '%s'", opts->mask_name);      nvox = DSET_NVOX(opts->inset);      if( DSET_NVOX(mset) != nvox ) {         ERROR_message("mask does not have the same voxel count as input");         RETURN(1);      }      /* fill mask array and mask_nxyz, remove mask dset */      DSET_load(mset); CHECK_LOAD_ERROR(mset);      opts->mask = THD_makemask(mset, 0, 1, 0);      DSET_delete(mset);      if( ! opts->mask ) {         ERROR_message("cannot make mask from '%s'", opts->mask_name);         RETURN(1);      }      if( opts->verb > 1 )         INFO_message("have mask with %d voxels", nvox);   }   RETURN(0);}

MRI_vectim * THD_dset_to_vectim_byslice( THD_3dim_dataset *dset, byte *mask ,                                         int ignore , int kzbot , int kztop  ){   byte *mmm ;   MRI_vectim *mrv=NULL ;   int kk,iv , nvals , nvox , nmask , nxy , nz ;ENTRY("THD_dset_to_vectim_byslice") ;                     if( !ISVALID_DSET(dset) ) RETURN(NULL) ;   DSET_load(dset) ; if( !DSET_LOADED(dset)  ) RETURN(NULL) ;   nvals = DSET_NVALS(dset) ; if( nvals <= 0 ) RETURN(NULL) ;   nvox  = DSET_NVOX(dset) ;   nxy = DSET_NX(dset) * DSET_NY(dset) ; nz = DSET_NZ(dset) ;   if( kzbot <  0  ) kzbot = 0 ;   if( kztop >= nz ) kztop = nz-1 ;   if( kztop < kzbot ) RETURN(NULL) ;   if( kzbot == 0 && kztop == nz-1 ){     mrv = THD_dset_to_vectim( dset , mask, ignore ) ; RETURN(mrv) ;   }   /* make a mask that includes cutting out un-desirable slices */   { int ibot , itop , ii ;#pragma omp critical (MALLOC)     mmm = (byte *)malloc(sizeof(byte)*nvox) ;     if( mask == NULL ) AAmemset( mmm ,    1 , sizeof(byte)*nvox ) ;     else               AAmemcpy( mmm , mask , sizeof(byte)*nvox ) ;     if( kzbot > 0 )       AAmemset( mmm               , 0 , sizeof(byte)*kzbot       *nxy ) ;     if( kztop < nz-1 )       AAmemset( mmm+(kztop+1)*nxy , 0 , sizeof(byte)*(nz-1-kztop)*nxy ) ;   }   /* and make the vectim using the standard function */   mrv = THD_dset_to_vectim( dset , mmm , ignore ) ;   free(mmm) ;   RETURN(mrv) ;}

int is_integral_sub_brick ( THD_3dim_dataset *dset, int isb, int check_values){   float mfac = 0.0;   void *vv=NULL;   if(   !ISVALID_DSET(dset)    ||            isb < 0                     ||            isb >= DSET_NVALS(dset)  ) {      fprintf(stderr,"** Bad dset or sub-brick index./n");      return (0) ;   }   if( !DSET_LOADED(dset) ) DSET_load(dset);   switch( DSET_BRICK_TYPE(dset,isb) ){      case MRI_short:      case MRI_byte:         if (check_values) {            mfac = DSET_BRICK_FACTOR(dset,isb) ;            if (mfac != 0.0f && mfac != 1.0f) return(0);         }         break;      case MRI_double:      case MRI_complex:      case MRI_float:         vv = (void *)DSET_ARRAY(dset,isb);         mfac = DSET_BRICK_FACTOR(dset,isb) ;         if (mfac != 0.0f && mfac != 1.0f) return(0);         if (!vv) {            fprintf(stderr,"** NULL array!/n");            return(0);         }         return(is_integral_data(DSET_NVOX(dset),                                 DSET_BRICK_TYPE(dset,isb),                                 DSET_ARRAY(dset,isb) ) );         break;      default:         return(0);   }   return(1);}

void THD_load_tcat( THD_datablock *dblk ){   int ivout , dd , iv ;   THD_3dim_dataset *dset_in , *dset_out ;   NI_str_array *sar ;ENTRY("THD_load_tcat") ;   if( !ISVALID_DBLK(dblk) ) EXRETURN ;   dset_out = (THD_3dim_dataset *)dblk->parent ;   if( !ISVALID_DSET(dset_out) ) EXRETURN ;   sar = NI_decode_string_list( dset_out->tcat_list , "~" ) ;   if( sar == NULL ) EXRETURN ;   if( sar->num != dset_out->tcat_num ){ NI_delete_str_array(sar); EXRETURN; }   ivout = 0 ;   for( dd=0 ; dd < sar->num ; dd++ ){     dset_in = THD_open_dataset( sar->str[dd] ) ;     if( dset_in == NULL ){       NI_delete_str_array(sar) ; DSET_unload(dset_out) ;       EXRETURN ;     }     DSET_mallocize(dset_in) ; DSET_load(dset_in) ;     if( !DSET_LOADED(dset_in) ){       NI_delete_str_array(sar) ; DSET_unload(dset_out) ; DSET_delete(dset_in) ;       EXRETURN ;     }     for( iv=0 ; iv < DSET_NVALS(dset_in) ; iv++ ){       EDIT_substitute_brick( dset_out , ivout ,                              DSET_BRICK_TYPE(dset_in,iv), DSET_ARRAY(dset_in,iv) );       mri_fix_data_pointer( NULL , DSET_BRICK(dset_in,iv) ) ;       EDIT_BRICK_FACTOR( dset_out , ivout , DSET_BRICK_FACTOR(dset_in,iv) ) ;       EDIT_BRICK_LABEL(dset_out, ivout,                         DSET_BRICK_LABEL(dset_in, iv)); /* ZSS Aug. 27 2012 */       ivout++ ;     }     DSET_delete(dset_in) ;   }   NI_delete_str_array(sar) ; EXRETURN ;}

MRI_IMAGE * THD_dset_to_1Dmri( THD_3dim_dataset *dset ){   MRI_IMAGE *im ; float *far ;   int nx , ny , ii ;ENTRY("THD_dset_to_1D") ;   if( !ISVALID_DSET(dset) ) RETURN(NULL) ;   DSET_load(dset) ;   if( !DSET_LOADED(dset) ) RETURN(NULL) ;   nx = DSET_NVALS(dset) ;   ny = DSET_NVOX(dset) ;   im = mri_new( nx , ny , MRI_float ) ; far = MRI_FLOAT_PTR(im) ;   for( ii=0 ; ii < ny ; ii++ )     THD_extract_array( ii , dset , 0 , far + ii*nx ) ;   RETURN(im) ;}

//.........这里部分代码省略.........          ERROR_message("Unknown argument on command line: '%s'",argv[iarg]) ;     suggest_best_prog_option(argv[0], argv[iarg]);     exit (1);   }   /*------- check options for completeness and consistency -----*/      if (nup == -1) {      if( iarg+1 >= argc )        ERROR_exit("need 'n' and 'dataset' on command line!") ;      nup = (int)strtod(argv[iarg++],NULL) ;      if( nup < 2 || nup > 320 )        ERROR_exit("3dUpsample rate '%d' is outside range 2..320",nup) ;   }    if (!dsetname) {      if( iarg >= argc )        ERROR_exit("need 'dataset' on command line!") ;      dsetname = argv[iarg];   }      inset = THD_open_dataset(dsetname) ;   if( !ISVALID_DSET(inset) )     ERROR_exit("3dUpsample can't open dataset '%s'", dsetname) ;   ntin = DSET_NVALS(inset) ; trin = DSET_TR(inset) ;   if( ntin < 2 )     ERROR_exit("dataset '%s' has only 1 value per voxel?!",dsetname) ;   nvox = DSET_NVOX(inset) ;   if( verb ) INFO_message("loading input dataset into memory") ;   DSET_load(inset) ; CHECK_LOAD_ERROR(inset) ;   /*------ create output dataset ------*/   ntout = ntin * nup ; trout = trin / nup ;   /* scaling factor for output */   fac = NULL; maxtop = 0.0;   if (MRI_IS_INT_TYPE(datum)) {      fac = (float *)calloc(DSET_NVALS(inset), sizeof(float));      ofac = (float *)calloc(ntout, sizeof(float));      for (ii=0; ii<DSET_NVALS(inset); ++ii) {         top = MCW_vol_amax( DSET_NVOX(inset),1,1 ,                              DSET_BRICK_TYPE(inset,ii),                              DSET_BRICK_ARRAY(inset,ii) ) ;         if (DSET_BRICK_FACTOR(inset, ii))             top = top * DSET_BRICK_FACTOR(inset,ii);         fac[ii] = (top > MRI_TYPE_maxval[datum]) ?                         top/MRI_TYPE_maxval[datum] : 0.0 ;         if (top > maxtop) maxtop = top;      }      if (storage_mode_from_filename(prefix) != STORAGE_BY_BRICK) {         fac[0] = (maxtop > MRI_TYPE_maxval[datum]) ?                         maxtop/MRI_TYPE_maxval[datum] : 0.0 ;         for (ii=0; ii<ntout; ++ii)             ofac[ii] = fac[0];         if (verb) INFO_message("Forcing global scaling, Max = %f, fac = %f/n",                         maxtop, fac[0]);      } else {         if (verb) INFO_message("Reusing scaling factors of input dset/n");         upsample_1( nup, DSET_NVALS(inset), fac, ofac);      }

/*! Replace a voxel's value by the value's rank in the entire set of input datasets */int main( int argc , char * argv[] ){   THD_3dim_dataset ** dsets_in = NULL, *dset=NULL; /*input and output datasets*/   int nopt=0, nbriks=0, nsubbriks=0, ib=0, isb=0;   byte *cmask=NULL;   int *all_uniques=NULL, **uniques=NULL, *final_unq=NULL, *N_uniques=NULL;   int N_final_unq=0, iun=0, total_unq=0;   INT_HASH_DATUM *rmap=NULL, *hd=NULL;   int imax=0, iunq=0, ii=0, id = 0;   long int off=0;   char *prefix=NULL;   char stmp[THD_MAX_PREFIX+1]={""};    FILE *fout=NULL;   /*----- Read command line -----*/   if( argc < 2 || strcmp(argv[1],"-help") == 0 ){      Rank_help ();      exit(0) ;   }   mainENTRY("3dRank main"); machdep(); AFNI_logger("3dRank",argc,argv);   nopt = 1 ;      while( nopt < argc && argv[nopt][0] == '-' ){      if( strcmp(argv[nopt],"-ver") == 0 ){         PRINT_VERSION("3dRank"); AUTHOR("Ziad Saad");         nopt++; continue;      }      if( strcmp(argv[nopt],"-help") == 0 ){         Rank_help();         exit(0) ;      }      if( strcmp(argv[nopt],"-prefix") == 0 ){         ++nopt;         if (nopt>=argc) {            fprintf(stderr,"**ERROR: Need string after -prefix/n");            exit(1);         }         prefix = argv[nopt] ;         ++nopt; continue;      }      if( strcmp(argv[nopt],"-input") == 0 ){         dsets_in = (THD_3dim_dataset**)                        calloc(argc-nopt+1, sizeof(THD_3dim_dataset*));         ++nopt; nbriks=0;         while (nopt < argc ) {            dsets_in[nbriks] = THD_open_dataset( argv[nopt] );            if( !ISVALID_DSET(dsets_in[nbriks]) ){              fprintf(stderr,"**ERROR: can't open dataset %s/n",argv[nopt]) ;              exit(1);            }            ++nopt; ++nbriks;          }         continue;      }            ERROR_exit( " Error - unknown option %s", argv[nopt]);   }    if (nopt < argc) {      ERROR_exit( " Error unexplained trailing option: %s/n", argv[nopt]);   }   if (!nbriks) {      ERROR_exit( " Error no volumes entered on command line?");   }      /* some checks and inits*/   nsubbriks = 0;   for (ib = 0; ib<nbriks; ++ib) {      if (!is_integral_dset(dsets_in[ib], 0)) {         ERROR_exit("Dset %s is not of an integral data type.",                         DSET_PREFIX(dsets_in[ib]));      }      nsubbriks += DSET_NVALS(dsets_in[ib]);   }      /* Now get unique arrays */   uniques = (int **)calloc(nsubbriks, sizeof(int*));   N_uniques = (int *)calloc(nsubbriks, sizeof(int));   total_unq = 0;   iun = 0;   for (ib = 0; ib<nbriks; ++ib) {      DSET_mallocize(dsets_in[ib]); DSET_load(dsets_in[ib]);      for (isb=0; isb<DSET_NVALS(dsets_in[ib]); ++isb) {         uniques[iun] = THD_unique_vals(dsets_in[ib], isb,                                        &(N_uniques[iun]), cmask);         total_unq += N_uniques[iun];          ++iun;      }   }      /* put all the arrays together and get the unique of the uniques */   all_uniques = (int *)calloc(total_unq, sizeof(int));   off=0;   for (iun=0; iun<nsubbriks; ++iun) {      memcpy(all_uniques+off, uniques[iun], N_uniques[iun]*sizeof(int));      off += N_uniques[iun];   }//.........这里部分代码省略.........

int main( int argc , char * argv[] ){   THD_3dim_dataset *dset ;   int iarg=1 ;   char *cc1="x",*cc2="y",*cc3="z" ;   float th1=0.0, th2=0.0, th3=0.0 ;   float thx,thy,thz ;   int   axx,ayy,azz ;   char *fname="testcox.ppm" , fn[128] ;   void * rhand ;   int bot=1 , ii , nim=0 ;   float omap[128] , bfac ;   MRI_IMAGE * im , * brim ;   int hbr[256] , nperc,ibot,itop,sum ;   byte * bar ;   double ctim ;   int imode=CREN_TWOSTEP ;   int pmode=CREN_SUM_VOX ;   if( argc < 2 || strcmp(argv[1],"-help") == 0 ){      printf("Usage: testcox [-rotate a b c] [-mip|-MIP] [-out f] [-bot b] [-nn|-ts|-li] dset/n") ;      exit(0) ;   }   enable_mcw_malloc() ;   while( iarg < argc && argv[iarg][0] == '-' ){      if( strcmp(argv[iarg],"-MIP") == 0 ){        pmode = CREN_MIP_VOX ; iarg++ ; continue ;      }      if( strcmp(argv[iarg],"-mip") == 0 ){        pmode = CREN_MINIP_VOX ; iarg++ ; continue ;      }      if( strcmp(argv[iarg],"-nn") == 0 ){         imode = CREN_NN ;         iarg++ ; continue ;      }      if( strcmp(argv[iarg],"-ts") == 0 ){         imode = CREN_TWOSTEP ;         iarg++ ; continue ;      }      if( strcmp(argv[iarg],"-li") == 0 ){         imode = CREN_LINEAR ;         iarg++ ; continue ;      }      if( strcmp(argv[iarg],"-bot") == 0 ){         bot = strtod( argv[++iarg] , NULL ) ;         iarg++ ; continue ;      }      if( strcmp(argv[iarg],"-rotate") == 0 ){         th1 = (PI/180.0) * strtod( argv[++iarg] , &cc1 ) ;         th2 = (PI/180.0) * strtod( argv[++iarg] , &cc2 ) ;         th3 = (PI/180.0) * strtod( argv[++iarg] , &cc3 ) ;         iarg++ ; continue ;      }      if( strcmp(argv[iarg],"-out") == 0 ){         fname = argv[++iarg] ;         iarg++ ; continue ;      }      fprintf(stderr,"Illegal option: %s/n",argv[iarg]); exit(1);   }   if( iarg >= argc ){fprintf(stderr,"No dataset?/n"); exit(1); }   dset = THD_open_dataset( argv[iarg] ) ;   if( dset == NULL ){fprintf(stderr,"Can't open dataset!/n");exit(1);}   if( DSET_BRICK_TYPE(dset,0) != MRI_byte ){      fprintf(stderr,"Non-byte dataset input!/n");exit(1);   }   DSET_mallocize(dset) ; DSET_load(dset) ;   if( !DSET_LOADED(dset) ){      fprintf(stderr,"Can't load dataset!/n");exit(1);   }   rhand = new_CREN_renderer() ;#if 0   THD_rotangle_user_to_dset( dset ,                              th1,*cc1  , th2,*cc2  , th3,*cc3 ,                              &thx,&axx , &thy,&ayy , &thz,&azz ) ;   CREN_set_viewpoint( rhand , axx,thx,ayy,thy,azz,thz ) ;#else   CREN_set_angles( rhand , th1,th2,th3 ) ;#endif   for( ii=0 ; ii < 128 ; ii++ )      omap[ii] = (ii <= bot) ? 0.0                             : (ii-bot)/(127.0-bot) ;   CREN_set_opamap( rhand , omap , 1.0 ) ;//.........这里部分代码省略.........

THD_3dim_dataset * THD_localhistog( int nsar , THD_3dim_dataset **insar ,                                    int numval , int *rlist , MCW_cluster *nbhd ,                                    int do_prob , int verb ){   THD_3dim_dataset *outset=NULL , *inset ;   int nvox=DSET_NVOX(insar[0]) ;   int ids, iv, bb, nnpt=nbhd->num_pt ;   MRI_IMAGE *bbim ; int btyp ;   float **outar , **listar ;ENTRY("THD_localhistog") ;   /*---- create output dataset ----*/   outset = EDIT_empty_copy(insar[0]) ;   EDIT_dset_items( outset ,                      ADN_nvals     , numval    ,                      ADN_datum_all , MRI_float ,                      ADN_nsl       , 0         ,                      ADN_brick_fac , NULL      ,                    ADN_none ) ;   outar = (float **)malloc(sizeof(float *)*numval) ;   for( bb=0 ; bb < numval ; bb++ ){     EDIT_substitute_brick( outset , bb , MRI_float , NULL ) ;     outar[bb] = DSET_BRICK_ARRAY(outset,bb) ;   }   /*---- make mapping between values and arrays to get those values ----*/   listar = (float **)malloc(sizeof(float *)*TWO16) ;   for( bb=0 ; bb < TWO16 ; bb++ ) listar[bb] = outar[0] ;   for( bb=1 ; bb < numval ; bb++ ){     listar[ rlist[bb] + TWO15 ] = outar[bb] ;   }   /*----------- loop over datasets, add in counts for all voxels -----------*/   for( ids=0 ; ids < nsar ; ids++ ){              /* dataset loop */     inset = insar[ids] ; DSET_load(inset) ;     for( iv=0 ; iv < DSET_NVALS(inset) ; iv++ ){  /* sub-brick loop */       if( verb ) fprintf(stderr,".") ;       bbim = DSET_BRICK(inset,iv) ; btyp = bbim->kind ;       if( nnpt == 1 ){                            /* only 1 voxel in nbhd */         int qq,ii,jj,kk,ib,nb ;         switch( bbim->kind ){           case MRI_short:{             short *sar = MRI_SHORT_PTR(bbim) ;             for( qq=0 ; qq < nvox ; qq++ ) listar[sar[qq]+TWO15][qq]++ ;           }           break ;           case MRI_byte:{             byte *bar = MRI_BYTE_PTR(bbim) ;             for( qq=0 ; qq < nvox ; qq++ ) listar[bar[qq]+TWO15][qq]++ ;           }           break ;           case MRI_float:{             float *far = MRI_FLOAT_PTR(bbim) ; short ss ;             for( qq=0 ; qq < nvox ; qq++ ){ ss = SHORTIZE(far[qq]); listar[ss+TWO15][qq]++; }           }           break ;         }       } else {                                    /* multiple voxels in nbhd */ AFNI_OMP_START ;#pragma omp parallel { int qq,ii,jj,kk,ib,nb ; void *nar ; short *sar,ss ; byte *bar ; float *far ;   nar = malloc(sizeof(float)*nnpt) ;   sar = (short *)nar ; bar = (byte *)nar ; far = (float *)nar ;#pragma omp for         for( qq=0 ; qq < nvox ; qq++ ){           /* qq=voxel index */           ii = DSET_index_to_ix(inset,qq) ;           jj = DSET_index_to_jy(inset,qq) ;           kk = DSET_index_to_kz(inset,qq) ;           nb = mri_get_nbhd_array( bbim , NULL , ii,jj,kk , nbhd , nar ) ;           if( nb == 0 ) continue ;           switch( btyp ){             case MRI_short:               for( ib=0 ; ib < nb ; ib++ ) listar[sar[ib]+TWO15][qq]++ ;             break ;             case MRI_byte:               for( ib=0 ; ib < nb ; ib++ ) listar[bar[ib]+TWO15][qq]++ ;             break ;             case MRI_float:               for( ib=0 ; ib < nb ; ib++ ){ ss = SHORTIZE(far[ib]); listar[ss+TWO15][qq]++; }             break ;           }         } /* end of voxel loop */   free(nar) ; } /* end of OpenMP */ AFNI_OMP_END ;       }     } /* end of sub-brick loop */     DSET_unload(inset) ;   } /* end of dataset loop */   if( verb ) fprintf(stderr,"/n") ;   free(listar) ;   /*---- post-process output ---*///.........这里部分代码省略.........

//.........这里部分代码省略.........   PLUTO_next_option(plint) ;                 /* skip to next line */   ignore = PLUTO_get_number(plint) ;         /* get number item (ignore) */   ninp = DSET_NUM_TIMES(old_dset) ;   /* number of values in input */   nuse = ninp;              /* number of values to actually use */   nfreq=nuse;   nfft=nuse;   str  = PLUTO_get_string(plint) ;              /* get string item (the datum type) */   istr = PLUTO_string_index( str ,              /* find it in the list it came from */                              NUM_TYPE_STRINGSX ,                              type_stringsx ) ;   switch( istr ){      default:      case 0: image_type = 0; break;           }  PLUTO_next_option(plint) ;                 /* skip to next line */  scale = PLUTO_get_number(plint) ;         /* get number item (scale) */   /*------------------------------------------------------*/   /*---------- At this point, the inputs are OK ----------*/   PLUTO_popup_meter( plint ) ;  /* popup a progress meter */   /*--------- set up pointers to each sub-brick in the input dataset ---------*/   DSET_load( old_dset ) ;  /* must be in memory before we get pointers to it */   old_datum = DSET_BRICK_TYPE( old_dset , 0 ) ; /* get old dataset datum type */   switch( old_datum ){  /* pointer type depends on input datum type */      default:         return "******************************/n"                "Illegal datum in Input Dataset/n"                "******************************"  ;      /** create array of pointers into old dataset sub-bricks **/      /** Note that we skip the first 'ignore' sub-bricks here **/      /*--------- input is bytes ----------*/      /* voxel #i at time #k is bptr[k][i] */      /* for i=0..nvox-1 and k=0..nuse-1.  */      case MRI_byte:         bptr = (byte **) malloc( sizeof(byte *) * nuse ) ;         if( bptr == NULL ) return "Malloc/nFailure!/n [bptr]" ;         for( kk=0 ; kk < nuse ; kk++ )            bptr[kk] = (byte *) DSET_ARRAY(old_dset,kk) ;      break ;      /*--------- input is shorts ---------*/      /* voxel #i at time #k is sptr[k][i] */      /* for i=0..nvox-1 and k=0..nuse-1.  */      case MRI_short:         sptr = (short **) malloc( sizeof(short *) * nuse ) ;         if( sptr == NULL ) return "Malloc/nFailure!/n [sptr]" ;         for( kk=0 ; kk < nuse ; kk++ )

//.........这里部分代码省略.........   /*---- check for stupid user inputs ----*/   if( iarg >= argc ) ERROR_exit("No datasets on command line?") ;   if( ohist_name == NULL && strcmp(prefix,"NULL") == 0 )     ERROR_exit("-prefix NULL is only meaningful if you also use -hsave :-(") ;   /*------------ scan input datasets, built overall histogram ------------*/   nsar  = argc - iarg ;   insar = (THD_3dim_dataset **)malloc(sizeof(THD_3dim_dataset *)*nsar) ;   if( verb ) fprintf(stderr,"Scanning %d datasets ",nsar) ;   ohist = (UINT32 *)calloc(sizeof(UINT32),TWO16) ;   for( ids=iarg ; ids < argc ; ids++ ){                      /* dataset loop */     insar[ids-iarg] = inset = THD_open_dataset(argv[ids]) ;     CHECK_OPEN_ERROR(inset,argv[ids]) ;     if( ids == iarg ){       nx = DSET_NX(inset); ny = DSET_NY(inset); nz = DSET_NZ(inset); nvox = nx*ny*nz;     } else if( nx != DSET_NX(inset) ||                ny != DSET_NY(inset) || nz != DSET_NZ(inset) ){       ERROR_exit("Dataset %s grid doesn't match!",argv[ids]) ;     }     if( !THD_datum_constant(inset->dblk) )       ERROR_exit("Dataset %s doesn't have a fixed data type! :-(",argv[ids]) ;     if( THD_need_brick_factor(inset) )       ERROR_exit("Dataset %s has scale factors! :-(",argv[ids]) ;     if( DSET_BRICK_TYPE(inset,0) != MRI_byte  &&         DSET_BRICK_TYPE(inset,0) != MRI_short &&         DSET_BRICK_TYPE(inset,0) != MRI_float    )       ERROR_exit("Dataset %s is not byte- or short-valued! :-(",argv[ids]) ;     DSET_load(inset) ; CHECK_LOAD_ERROR(inset) ;     for( ii=0 ; ii < DSET_NVALS(inset) ; ii++ ){ /* add to overall histogram */       if( verb ) fprintf(stderr,".") ;       switch( DSET_BRICK_TYPE(inset,ii) ){         case MRI_short:{           short *sar = (short *)DSET_BRICK_ARRAY(inset,ii) ;           for( kk=0 ; kk < nvox ; kk++ ) ohist[ sar[kk]+TWO15 ]++ ;         }         break ;         case MRI_byte:{           byte *bar = (byte *)DSET_BRICK_ARRAY(inset,ii) ;           for( kk=0 ; kk < nvox ; kk++ ) ohist[ bar[kk]+TWO15 ]++ ;         }         break ;         case MRI_float:{           float *far = (float *)DSET_BRICK_ARRAY(inset,ii) ; short ss ;           for( kk=0 ; kk < nvox ; kk++ ){ ss = SHORTIZE(far[kk]); ohist[ss+TWO15]++; }         }         break ;       }     } /* end of sub-brick loop */     DSET_unload(inset) ;  /* will re-load later, as needed */   } /* end of dataset loop */   if( verb ) fprintf(stderr,"/n") ;   /*-------------- process overall histogram for fun and profit -------------*/   /* if we didn't actually find 0, put it in the histogram now */

static int * PLUTO_4D_to_nothing (THD_3dim_dataset * old_dset , int ignore , int detrend ,                         generic_func * user_func, void * user_data ){   byte    ** bptr = NULL ;  /* one of these will be the array of */   short   ** sptr = NULL ;  /* pointers to input dataset sub-bricks */   float   ** fptr = NULL ;  /* (depending on input datum type) */   complex ** cptr = NULL ;   float * fxar = NULL ;  /* array loaded from input dataset */   float * fac  = NULL ;  /* array of brick scaling factors */   float * dtr  = NULL ;  /* will be array of detrending coeff */   float val , d0fac , d1fac , x0,x1;   double tzero=0.0 , tdelta , ts_mean , ts_slope ;   int   ii , old_datum , nuse , use_fac , iz,izold, nxy,nvox ;   static int retval;	register int kk ;   /*----------------------------------------------------------*/   /*----- Check inputs to see if they are reasonable-ish -----*/   if( ! ISVALID_3DIM_DATASET(old_dset) ) return NULL ;   if( user_func == NULL ) return NULL ;   if( ignore < 0 ) ignore = 0 ;   /*--------- set up pointers to each sub-brick in the input dataset ---------*/   old_datum = DSET_BRICK_TYPE( old_dset , 0 ) ;   /* get old dataset datum */   nuse      = DSET_NUM_TIMES(old_dset) - ignore ; /* # of points on time axis */   if( nuse < 2 ) return NULL ;   DSET_load( old_dset ) ;  /* must be in memory before we get pointers to it */   kk = THD_count_databricks( old_dset->dblk ) ;  /* check if it was */   if( kk < DSET_NVALS(old_dset) ){               /* loaded correctly */      DSET_unload( old_dset ) ;      return NULL ;   }   switch( old_datum ){  /* pointer type depends on input datum type */      default:                      /** don't know what to do **/         DSET_unload( old_dset ) ;         return NULL ;      /** create array of pointers into old dataset sub-bricks **/      /*--------- input is bytes ----------*/      /* voxel #i at time #k is bptr[k][i] */      /* for i=0..nvox-1 and k=0..nuse-1.  */      case MRI_byte:         bptr = (byte **) malloc( sizeof(byte *) * nuse ) ;         if( bptr == NULL ) return NULL ;         for( kk=0 ; kk < nuse ; kk++ )            bptr[kk] = (byte *) DSET_ARRAY(old_dset,kk+ignore) ;      break ;      /*--------- input is shorts ---------*/      /* voxel #i at time #k is sptr[k][i] */      /* for i=0..nvox-1 and k=0..nuse-1.  */      case MRI_short:         sptr = (short **) malloc( sizeof(short *) * nuse ) ;         if( sptr == NULL ) return NULL ;         for( kk=0 ; kk < nuse ; kk++ )            sptr[kk] = (short *) DSET_ARRAY(old_dset,kk+ignore) ;      break ;      /*--------- input is floats ---------*/      /* voxel #i at time #k is fptr[k][i] */      /* for i=0..nvox-1 and k=0..nuse-1.  */      case MRI_float:         fptr = (float **) malloc( sizeof(float *) * nuse ) ;         if( fptr == NULL ) return NULL ;         for( kk=0 ; kk < nuse ; kk++ )            fptr[kk] = (float *) DSET_ARRAY(old_dset,kk+ignore) ;      break ;      /*--------- input is complex ---------*/      /* voxel #i at time #k is cptr[k][i]  */      /* for i=0..nvox-1 and k=0..nuse-1.   */      case MRI_complex:         cptr = (complex **) malloc( sizeof(complex *) * nuse ) ;         if( cptr == NULL ) return NULL ;         for( kk=0 ; kk < nuse ; kk++ )            cptr[kk] = (complex *) DSET_ARRAY(old_dset,kk+ignore) ;      break ;   } /* end of switch on input type */	nvox = old_dset->daxes->nxx * old_dset->daxes->nyy * old_dset->daxes->nzz ;      /*---- allocate space for 1 voxel timeseries ----*///.........这里部分代码省略.........

//.........这里部分代码省略.........		if( strcmp(argv[nopt],"-blur") == 0 ){			if( ++nopt >= argc ) ERROR_exit("need an argument after -blur!") ;			blur = strtod(argv[nopt],NULL) ;			if( blur <= 0.0f ) WARNING_message("non-positive blur?!") ;			nopt++ ; continue ;		}		if( strcmp(argv[nopt],"-localPV") == 0 ){			if( ++nopt >= argc ) ERROR_exit("need an argument after -localpv!") ;			pvrad = strtod(argv[nopt],NULL) ;			if( pvrad <= 0.0f ) WARNING_message("non-positive -localpv?!") ;			nopt++ ; continue ;		}		if( strcmp(argv[nopt],"-prefix") == 0 ){			if( ++nopt >= argc ) ERROR_exit("need an argument after -prefix!") ;			prefix = strdup(argv[nopt]) ;			if( !THD_filename_ok(prefix) ) ERROR_exit("bad -prefix option!") ;			nopt++ ; continue ;		}		if( strcmp(argv[nopt],"-automask") == 0 ){			if( mask != NULL ) ERROR_exit("Can't use -mask AND -automask!") ;			do_automask = 1 ; nopt++ ; continue ;		}		if( strcmp(argv[nopt],"-mask") == 0 ){			THD_3dim_dataset *mset ;			if( ++nopt >= argc ) ERROR_exit("Need argument after '-mask'") ;			if( mask != NULL || do_automask ) ERROR_exit("Can't have two mask inputs") ;			mset = THD_open_dataset( argv[nopt] ) ;			CHECK_OPEN_ERROR(mset,argv[nopt]) ;			DSET_load(mset) ; CHECK_LOAD_ERROR(mset) ;			mask_nx = DSET_NX(mset); mask_ny = DSET_NY(mset); mask_nz = DSET_NZ(mset);			mask = THD_makemask( mset , 0 , 0.5f, 0.0f ) ; DSET_delete(mset) ;			if( mask == NULL ) ERROR_exit("Can't make mask from dataset '%s'",argv[nopt]) ;			nmask = THD_countmask( mask_nx*mask_ny*mask_nz , mask ) ;			if( verb ) INFO_message("Number of voxels in mask = %d",nmask) ;			if( nmask < 1 ) ERROR_exit("Mask is too small to process") ;			nopt++ ; continue ;		}		if( strcmp(argv[nopt],"-norm") == 0 ){			do_norm = 1 ; nopt++ ; continue ;		}		if( strcmp(argv[nopt],"-quiet") == 0 ){			verb = 0 ; nopt++ ; continue ;		}		if( strcmp(argv[nopt],"-no_rs_out") == 0 ){ // @@			SERIES_OUT = 0 ; nopt++ ; continue ;		}		if( strcmp(argv[nopt],"-un_bp_out") == 0 ){ // @@			UNBP_OUT = 1 ; nopt++ ; continue ;		}		if( strcmp(argv[nopt],"-no_rsfa") == 0 ){ // @@			DO_RSFA = 0 ; nopt++ ; continue ;		}		if( strcmp(argv[nopt],"-bp_at_end") == 0 ){ // @@			BP_LAST = 1 ; nopt++ ; continue ;		}

//.........这里部分代码省略.........         int val = (int)strtod(argv[++nopt],&cpt) ;         if( *cpt != '/0' || val < -1 || val > 3 ){            ERROR_exit("Illegal value after -polort!") ;         }         polort = val ; nopt++ ; continue ;      }      if( strcmp(argv[nopt],"-mem_stat") == 0 ){         MEM_STAT = 1 ; nopt++ ; continue ;      }      if( strncmp(argv[nopt],"-mem_profile",8) == 0 ){         MEM_PROF = 1 ; nopt++ ; continue ;      }      /* check for 1d argument */      if ( strcmp(argv[nopt],"-out1D") == 0 ){          if (!(fout1D = fopen(argv[++nopt], "w"))) {             ERROR_message("Failed to open %s for writing", argv[nopt]);             exit(1);          }          nopt++ ; continue ;      }      ERROR_exit("Illegal option: %s",argv[nopt]) ;   }   /*-- open dataset, check for legality --*/   if( nopt >= argc ) ERROR_exit("Need a dataset on command line!?") ;   xset = THD_open_dataset(argv[nopt]); CHECK_OPEN_ERROR(xset,argv[nopt]);   if( DSET_NVALS(xset) < 3 )     ERROR_exit("Input dataset %s does not have 3 or more sub-bricks!",argv[nopt]) ;   DSET_load(xset) ; CHECK_LOAD_ERROR(xset) ;   /*-- compute mask array, if desired --*/   nvox = DSET_NVOX(xset) ; nvals = DSET_NVALS(xset) ;   INC_MEM_STATS((nvox * nvals * sizeof(double)), "input dset");   PRINT_MEM_STATS("inset");   /* if a mask was specified make sure it is appropriate */   if( mset ){      if( DSET_NVOX(mset) != nvox )         ERROR_exit("Input and mask dataset differ in number of voxels!") ;      mask  = THD_makemask(mset, 0, 1.0, 0.0) ;      /* update running memory statistics to reflect loading the image */      INC_MEM_STATS( mset->dblk->total_bytes, "mask dset" );      PRINT_MEM_STATS( "mset load" );      nmask = THD_countmask( nvox , mask ) ;      INC_MEM_STATS( nmask * sizeof(byte), "mask array" );      PRINT_MEM_STATS( "mask" );      INFO_message("%d voxels in -mask dataset",nmask) ;      if( nmask < 2 ) ERROR_exit("Only %d voxels in -mask, exiting...",nmask);      /* update running memory statistics to reflect loading the image */      DEC_MEM_STATS( mset->dblk->total_bytes, "mask dset" );      DSET_unload(mset) ;      PRINT_MEM_STATS( "mset unload" );   }    /* if automasking is requested, handle that now */   else if( do_autoclip ){      mask  = THD_automask( xset ) ;

//.........这里部分代码省略.........      }      if( strcmp(argv[iarg],"-prefix") == 0 ){         prefix = argv[++iarg] ;         if( !THD_filename_ok(prefix) ){            fprintf(stderr,"*** Illegal string after -prefix!/n"); exit(1) ;         }         iarg++ ; continue ;      }      if( strcmp(argv[iarg],"-maxstep") == 0 ){         maxgap = strtol( argv[++iarg] , NULL , 10 ) ;         if( maxgap < 1 ){            fprintf(stderr,"*** Illegal value after -maxgap!/n"); exit(1);         }         iarg++ ; continue ;      }      if( strcmp(argv[iarg],"-dir") == 0 ){         dstr = argv[++iarg] ;         iarg++ ; continue ;      }      fprintf(stderr,"*** Illegal option: %s/n",argv[iarg]) ; exit(1) ;   }   if( dstr == NULL ){      fprintf(stderr,"*** No -dir option on command line!/n"); exit(1);   }   if( iarg >= argc ){      fprintf(stderr,"*** No input dataset on command line!/n"); exit(1);   }   inset = THD_open_dataset( argv[iarg] ) ;   if( inset == NULL ){      fprintf(stderr,"*** Can't open dataset %s/n",argv[iarg]); exit(1);   }   outset = EDIT_empty_copy( inset ) ;   EDIT_dset_items( outset , ADN_prefix , prefix , ADN_none ) ;   if( THD_deathcon() && THD_is_file( DSET_HEADNAME(outset) ) ){      fprintf(stderr,"** Output file %s exists -- cannot overwrite!/n",              DSET_HEADNAME(outset) ) ;      exit(1) ;   }   tross_Copy_History( inset , outset ) ;   tross_Make_History( "3dZFillin" , argc,argv , outset ) ;   if( DSET_NVALS(inset) > 1 ){      fprintf(stderr,"++ WARNING: input dataset has more than one sub-brick!/n");      EDIT_dset_items( outset ,                         ADN_ntt   , 0 ,                         ADN_nvals , 1 ,                       ADN_none ) ;   }   if( DSET_BRICK_TYPE(outset,0) != MRI_byte ){      fprintf(stderr,"*** This program only works on byte datasets!/n");      exit(1) ;   }   switch( *dstr ){      case 'x': dcode = 1 ; break ;      case 'y': dcode = 2 ; break ;      case 'z': dcode = 3 ; break ;      default:        if( *dstr == ORIENT_tinystr[outset->daxes->xxorient][0] ||            *dstr == ORIENT_tinystr[outset->daxes->xxorient][1]   ) dcode = 1 ;        if( *dstr == ORIENT_tinystr[outset->daxes->yyorient][0] ||            *dstr == ORIENT_tinystr[outset->daxes->yyorient][1]   ) dcode = 2 ;        if( *dstr == ORIENT_tinystr[outset->daxes->zzorient][0] ||            *dstr == ORIENT_tinystr[outset->daxes->zzorient][1]   ) dcode = 3 ;      break ;   }   if( dcode == 0 ){      fprintf(stderr,"*** Illegal -dir direction!/n") ; exit(1) ;   }   if( verb )      fprintf(stderr,"++ Direction = axis %d in dataset/n",dcode) ;   DSET_load(inset) ; CHECK_LOAD_ERROR(inset) ;   brim = mri_copy( DSET_BRICK(inset,0) ) ;   DSET_unload(inset) ;   EDIT_substitute_brick( outset , 0 , brim->kind , mri_data_pointer(brim) ) ;   nftot = THD_dataset_zfillin( outset , 0 , dcode , maxgap ) ;   fprintf(stderr,"++ Number of voxels filled = %d/n",nftot) ;   if (DSET_write(outset) != False) {      fprintf(stderr,"++ output dataset: %s/n",DSET_BRIKNAME(outset)) ;      exit(0) ;   } else {      fprintf(stderr,         "** 3dZFillin: Failed to write output!/n" ) ;      exit(1) ;   }   }

int main( int argc , char *argv[] ){   THD_3dim_dataset *inset=NULL , *outset=NULL ;   MCW_cluster *nbhd=NULL ;   byte *mask=NULL ; int mask_nx,mask_ny,mask_nz , automask=0 ;   char *prefix="./LocalCormat" ;   int iarg=1 , verb=1 , ntype=0 , kk,nx,ny,nz,nxy,nxyz,nt , xx,yy,zz, vstep ;   float na,nb,nc , dx,dy,dz ;   MRI_IMARR *imar=NULL ; MRI_IMAGE *pim=NULL ;   int mmlag=10 , ii,jj , do_arma=0 , nvout ;   MRI_IMAGE *concim=NULL ; float *concar=NULL ;   if( argc < 2 || strcmp(argv[1],"-help") == 0 ){     printf(       "Usage: 3dLocalCORMAT [options] inputdataset/n"       "/n"       "Compute the correlation matrix (in time) of the input dataset,/n"       "up to lag given by -maxlag.  The matrix is averaged over the/n"       "neighborhood specified by the -nbhd option, and then the entries/n"       "are output at each voxel in a new dataset./n"       "/n"       "Normally, the input to this program would be the -errts output/n"       "from 3dDeconvolve, or the equivalent residuals from some other/n"       "analysis.  If you input a non-residual time series file, you at/n"       "least should use an appropriate -polort level for detrending!/n"       "/n"       "Options:/n"       "  -input inputdataset/n"       "  -prefix ppp/n"       "  -mask mset    {these 2 options are}/n"       "  -automask     {mutually exclusive.}/n"       "  -nbhd nnn     [e.g., 'SPHERE(9)' for 9 mm radius]/n"       "  -polort ppp   [default = 0, which is reasonable for -errts output]/n"       "  -concat ccc   [as in 3dDeconvolve]/n"       "  -maxlag mmm   [default = 10]/n"       "  -ARMA         [estimate ARMA(1,1) parameters into last 2 sub-bricks]/n"       "/n"       "A quick hack for my own benignant purposes -- RWCox -- June 2008/n"     ) ;     PRINT_COMPILE_DATE ; exit(0) ;   }   /*---- official startup ---*/   PRINT_VERSION("3dLocalCormat"); mainENTRY("3dLocalCormat main"); machdep();   AFNI_logger("3dLocalCormat",argc,argv); AUTHOR("Zhark the Toeplitzer");   /*---- loop over options ----*/   while( iarg < argc && argv[iarg][0] == '-' ){#if 0fprintf(stderr,"argv[%d] = %s/n",iarg,argv[iarg]) ;#endif     if( strcmp(argv[iarg],"-ARMA") == 0 ){       do_arma = 1 ; iarg++ ; continue ;     }     if( strcmp(argv[iarg],"-polort") == 0 ){       char *cpt ;       if( ++iarg >= argc )         ERROR_exit("Need argument after option %s",argv[iarg-1]) ;       pport = (int)strtod(argv[iarg],&cpt) ;       if( *cpt != '/0' )         WARNING_message("Illegal non-numeric value after -polort") ;       if( pport > 3 ){         pport = 3 ; WARNING_message("-polort set to 3 == max implemented") ;       } else if( pport < 0 ){         pport = 0 ; WARNING_message("-polort set to 0 == min implemented") ;       }       iarg++ ; continue ;     }     if( strcmp(argv[iarg],"-input") == 0 ){       if( inset != NULL  ) ERROR_exit("Can't have two -input options") ;       if( ++iarg >= argc ) ERROR_exit("Need argument after '-input'") ;       inset = THD_open_dataset( argv[iarg] ) ;       CHECK_OPEN_ERROR(inset,argv[iarg]) ;       iarg++ ; continue ;     }     if( strcmp(argv[iarg],"-prefix") == 0 ){       if( ++iarg >= argc ) ERROR_exit("Need argument after '-prefix'") ;       prefix = strdup(argv[iarg]) ;       iarg++ ; continue ;     }     if( strcmp(argv[iarg],"-mask") == 0 ){       THD_3dim_dataset *mset ; int mmm ;       if( ++iarg >= argc ) ERROR_exit("Need argument after '-mask'") ;       if( mask != NULL || automask ) ERROR_exit("Can't have two mask inputs") ;       mset = THD_open_dataset( argv[iarg] ) ;       CHECK_OPEN_ERROR(mset,argv[iarg]) ;       DSET_load(mset) ; CHECK_LOAD_ERROR(mset) ;       mask_nx = DSET_NX(mset); mask_ny = DSET_NY(mset); mask_nz = DSET_NZ(mset);       mask = THD_makemask( mset , 0 , 0.5f, 0.0f ) ; DSET_delete(mset) ;       if( mask == NULL ) ERROR_exit("Can't make mask from dataset '%s'",argv[iarg]) ;       mmm = THD_countmask( mask_nx*mask_ny*mask_nz , mask ) ;       INFO_message("Number of voxels in mask = %d",mmm) ;//.........这里部分代码省略.........