do{ int jv ;                                               \
      if( fout != NULL ){                                    \
         for( jv=0 ; jv < nbrik ; jv++ ) FREEUP(fout[jv]) ;  \
         free(fout) ;                                        \
      } } while(0)

do{ FREEUP(bptr) ; FREEUP(sptr) ; FREEUP(fptr) ;  \
      FREEUP(cptr) ; FREEUP(fxar) ; FREEUP(fac)  ;  \
      FREEUP(dtr)  ; FREEUP(val)  ; FREE_FOUT    ;  \
    } while(0)

00098 {
00099    THD_3dim_dataset * new_dset ;  /* output dataset */
00100 
00101    byte    ** bptr = NULL ;  /* one of these will be the array of */
00102    short   ** sptr = NULL ;  /* pointers to input dataset sub-bricks */
00103    float   ** fptr = NULL ;  /* (depending on input datum type) */
00104    complex ** cptr = NULL ;
00105 
00106    float *  fxar = NULL ;  /* array loaded from input dataset */
00107    float *  fac  = NULL ;  /* array of input brick scaling factors */
00108    float ** fout = NULL ;  /* will be arrays of output floats */
00109    float *  dtr  = NULL ;  /* will be array of detrending coeff */
00110    float *  val  = NULL ;  /* will be array of output values */
00111 
00112    float d0fac , d1fac , x0,x1;
00113    double tzero , tdelta , ts_mean , ts_slope ;
00114    int   ii , old_datum , nuse , use_fac , iz,izold, nxy,nvox , iv ;
00115    register int kk ;
00116    int nbad=0 ;        /* 08 Aug 2000 */
00117 
00118    void (*ufunc)(double,double,int,float *,double,double,void *,int,float *)
00119      = (void (*)(double,double,int,float *,double,double,void *,int,float *)) user_func;
00120 
00121    /*----------------------------------------------------------*/
00122    /*----- Check inputs to see if they are reasonable-ish -----*/
00123 
00124    if( ! ISVALID_3DIM_DATASET(old_dset) ) return NULL ;
00125 
00126    if( new_datum >= 0         &&
00127        new_datum != MRI_byte  &&
00128        new_datum != MRI_short &&
00129        new_datum != MRI_float   ) return NULL ;
00130 
00131    if( user_func == NULL ) return NULL ;
00132 
00133    if( nbrik <= 0 ) return NULL ;
00134 
00135    if( ignore < 0 ) ignore = 0 ;
00136 
00137    /*--------- set up pointers to each sub-brick in the input dataset ---------*/
00138 
00139    old_datum = DSET_BRICK_TYPE( old_dset , 0 ) ;   /* get old dataset datum */
00140    nuse      = DSET_NUM_TIMES(old_dset) - ignore ; /* # of points on time axis */
00141    if( nuse < 2 ) return NULL ;
00142 
00143    if( new_datum < 0 ) new_datum = old_datum ;   /* output datum = input */
00144    if( new_datum == MRI_complex ) return NULL ;  /* but complex = bad news */
00145 
00146    DSET_load( old_dset ) ;  /* must be in memory before we get pointers to it */
00147 
00148    kk = THD_count_databricks( old_dset->dblk ) ;  /* check if it was */
00149    if( kk < DSET_NVALS(old_dset) ){               /* loaded correctly */
00150       DSET_unload( old_dset ) ;
00151       return NULL ;
00152    }
00153 
00154    switch( old_datum ){  /* pointer type depends on input datum type */
00155 
00156       default:                      /** don't know what to do **/
00157          DSET_unload( old_dset ) ;
00158          return NULL ;
00159 
00160       /** create array of pointers into old dataset sub-bricks **/
00161 
00162       /*--------- input is bytes ----------*/
00163       /* voxel #i at time #k is bptr[k][i] */
00164       /* for i=0..nvox-1 and k=0..nuse-1.  */
00165 
00166       case MRI_byte:
00167          bptr = (byte **) malloc( sizeof(byte *) * nuse ) ;
00168          if( bptr == NULL ) return NULL ;
00169          for( kk=0 ; kk < nuse ; kk++ )
00170             bptr[kk] = (byte *) DSET_ARRAY(old_dset,kk+ignore) ;
00171       break ;
00172 
00173       /*--------- input is shorts ---------*/
00174       /* voxel #i at time #k is sptr[k][i] */
00175       /* for i=0..nvox-1 and k=0..nuse-1.  */
00176 
00177       case MRI_short:
00178          sptr = (short **) malloc( sizeof(short *) * nuse ) ;
00179          if( sptr == NULL ) return NULL ;
00180          for( kk=0 ; kk < nuse ; kk++ )
00181             sptr[kk] = (short *) DSET_ARRAY(old_dset,kk+ignore) ;
00182       break ;
00183 
00184       /*--------- input is floats ---------*/
00185       /* voxel #i at time #k is fptr[k][i] */
00186       /* for i=0..nvox-1 and k=0..nuse-1.  */
00187 
00188       case MRI_float:
00189          fptr = (float **) malloc( sizeof(float *) * nuse ) ;
00190          if( fptr == NULL ) return NULL ;
00191          for( kk=0 ; kk < nuse ; kk++ )
00192             fptr[kk] = (float *) DSET_ARRAY(old_dset,kk+ignore) ;
00193       break ;
00194 
00195       /*--------- input is complex ---------*/
00196       /* voxel #i at time #k is cptr[k][i]  */
00197       /* for i=0..nvox-1 and k=0..nuse-1.   */
00198 
00199       case MRI_complex:
00200          cptr = (complex **) malloc( sizeof(complex *) * nuse ) ;
00201          if( cptr == NULL ) return NULL ;
00202          for( kk=0 ; kk < nuse ; kk++ )
00203             cptr[kk] = (complex *) DSET_ARRAY(old_dset,kk+ignore) ;
00204       break ;
00205 
00206    } /* end of switch on input type */
00207 
00208    /*---- allocate space for 1 voxel timeseries ----*/
00209 
00210    fxar = (float *) malloc( sizeof(float) * nuse ) ;   /* voxel timeseries */
00211    if( fxar == NULL ){ FREE_WORKSPACE ; return NULL ; }
00212 
00213    /*--- get scaling factors for input sub-bricks ---*/
00214 
00215    fac = (float *) malloc( sizeof(float) * nuse ) ;   /* factors */
00216    if( fac == NULL ){ FREE_WORKSPACE ; return NULL ; }
00217 
00218    use_fac = 0 ;
00219    for( kk=0 ; kk < nuse ; kk++ ){
00220       fac[kk] = DSET_BRICK_FACTOR(old_dset,kk+ignore) ;
00221       if( fac[kk] != 0.0 ) use_fac++ ;
00222       else                 fac[kk] = 1.0 ;
00223    }
00224    if( !use_fac ) FREEUP(fac) ;
00225 
00226    /*--- setup for detrending ---*/
00227 
00228    dtr = (float *) malloc( sizeof(float) * nuse ) ;
00229    if( dtr == NULL ){ FREE_WORKSPACE ; return NULL ; }
00230 
00231    d0fac = 1.0 / nuse ;
00232    d1fac = 12.0 / nuse / (nuse*nuse - 1.0) ;
00233    for( kk=0 ; kk < nuse ; kk++ )
00234       dtr[kk] = kk - 0.5 * (nuse-1) ;  /* linear trend, orthogonal to 1 */
00235 
00236    /*---------------------- make a new dataset ----------------------*/
00237 
00238    new_dset = EDIT_empty_copy( old_dset ) ; /* start with copy of old one */
00239 
00240    /*-- edit some of its internal parameters --*/
00241 
00242    ii = EDIT_dset_items(
00243            new_dset ,
00244               ADN_prefix      , new_prefix ,           /* filename prefix */
00245               ADN_malloc_type , DATABLOCK_MEM_MALLOC , /* store in memory */
00246               ADN_datum_all   , new_datum ,            /* atomic datum */
00247               ADN_nvals       , nbrik ,                /* # sub-bricks */
00248               ADN_ntt         , 0 ,                    /* # time points */
00249               ADN_type        , ISHEAD(old_dset)       /* dataset type */
00250                                  ? HEAD_FUNC_TYPE
00251                                  : GEN_FUNC_TYPE ,
00252               ADN_func_type   , FUNC_BUCK_TYPE ,        /* function type */
00253            ADN_none ) ;
00254 
00255    if( ii != 0 ){
00256       THD_delete_3dim_dataset( new_dset , False ) ;  /* some error above */
00257       FREE_WORKSPACE ; return NULL ;
00258    }
00259 
00260    THD_init_datablock_labels( new_dset->dblk ) ;
00261    THD_init_datablock_keywords( new_dset->dblk ) ;
00262    THD_init_datablock_stataux( new_dset->dblk ) ;
00263 
00264    /*------ make floating point output bricks
00265             (only at the end will scale to byte or shorts) ------*/
00266 
00267    nvox = old_dset->daxes->nxx * old_dset->daxes->nyy * old_dset->daxes->nzz ;
00268 
00269    fout = (float **) malloc( sizeof(float *) * nbrik ) ;
00270 
00271    if( fout == NULL ){
00272       THD_delete_3dim_dataset( new_dset , False ) ;
00273       FREE_WORKSPACE ; return NULL ;
00274    }
00275 
00276    for( iv=0 ; iv < nbrik ; iv++ ) fout[iv] = NULL ;
00277 
00278    for( iv=0 ; iv < nbrik ; iv++ ){
00279       fout[iv] = (float *) malloc( sizeof(float) * nvox ) ;
00280       if( fout[iv] == NULL ){
00281          THD_delete_3dim_dataset( new_dset , False ) ;
00282          FREE_WORKSPACE ; return NULL ;
00283       }
00284    }
00285 
00286    /*-- floating point storage for output from 1 voxel --*/
00287 
00288    val = (float *) malloc( sizeof(float) * nbrik ) ;
00289    if( val == NULL ){
00290       THD_delete_3dim_dataset( new_dset , False ) ;
00291       FREE_WORKSPACE ; return NULL ;
00292    }
00293 
00294    /*----- set up to find time at each voxel -----*/
00295 
00296    tdelta = old_dset->taxis->ttdel ;
00297    if( DSET_TIMEUNITS(old_dset) == UNITS_MSEC_TYPE ) tdelta *= 0.001 ;
00298    if( tdelta == 0.0 ) tdelta = 1.0 ;
00299 
00300    izold  = -666 ;
00301    nxy    = old_dset->daxes->nxx * old_dset->daxes->nyy ;
00302 
00303    /*----------------------------------------------------*/
00304    /*----- Setup has ended.  Now do some real work. -----*/
00305 
00306    /* start notification */
00307 
00308 #if 0
00309    user_func(  0.0 , 0.0 , nvox , NULL,0.0,0.0 , user_data , nbrik , NULL ) ;
00310 #else
00311    ufunc(  0.0 , 0.0 , nvox , NULL,0.0,0.0 , user_data , nbrik , NULL ) ;
00312 #endif
00313 
00314    /***** loop over voxels *****/
00315 
00316    for( ii=0 ; ii < nvox ; ii++  ){  /* 1 time series at a time */
00317 
00318       /*** load data from input dataset, depending on type ***/
00319 
00320       switch( old_datum ){
00321 
00322          /*** input = bytes ***/
00323 
00324          case MRI_byte:
00325             for( kk=0 ; kk < nuse ; kk++ ) fxar[kk] = bptr[kk][ii] ;
00326          break ;
00327 
00328          /*** input = shorts ***/
00329 
00330          case MRI_short:
00331             for( kk=0 ; kk < nuse ; kk++ ) fxar[kk] = sptr[kk][ii] ;
00332          break ;
00333 
00334          /*** input = floats ***/
00335 
00336          case MRI_float:
00337             for( kk=0 ; kk < nuse ; kk++ ) fxar[kk] = fptr[kk][ii] ;
00338          break ;
00339 
00340          /*** input = complex (note we use absolute value) ***/
00341 
00342          case MRI_complex:
00343             for( kk=0 ; kk < nuse ; kk++ ) fxar[kk] = CABS(cptr[kk][ii]) ;
00344          break ;
00345 
00346       } /* end of switch over input type */
00347 
00348       /*** scale? ***/
00349 
00350       if( use_fac )
00351          for( kk=0 ; kk < nuse ; kk++ ) fxar[kk] *= fac[kk] ;
00352 
00353       /** compute mean and slope **/
00354 
00355       x0 = x1 = 0.0 ;
00356       for( kk=0 ; kk < nuse ; kk++ ){
00357          x0 += fxar[kk] ; x1 += fxar[kk] * dtr[kk] ;
00358       }
00359 
00360       x0 *= d0fac ; x1 *= d1fac ;  /* factors to remove mean and trend */
00361 
00362       ts_mean  = x0 ;
00363       ts_slope = x1 / tdelta ;
00364 
00365       /** detrend? **/
00366 
00367       if( detrend )
00368          for( kk=0 ; kk < nuse ; kk++ ) fxar[kk] -= (x0 + x1 * dtr[kk]) ;
00369 
00370       /** compute start time of this timeseries **/
00371 
00372       iz = ii / nxy ;    /* which slice am I in? */
00373 
00374       if( iz != izold ){          /* in a new slice? */
00375          tzero = THD_timeof( ignore ,
00376                              old_dset->daxes->zzorg
00377                            + iz*old_dset->daxes->zzdel , old_dset->taxis ) ;
00378          izold = iz ;
00379 
00380          if( DSET_TIMEUNITS(old_dset) == UNITS_MSEC_TYPE ) tzero *= 0.001 ;
00381       }
00382 
00383       /*** compute output ***/
00384 
00385       for( iv=0 ; iv < nbrik ; iv++ ) val[iv] = 0.0 ;
00386 
00387 #if 0
00388       user_func( tzero,tdelta, nuse,fxar,ts_mean,ts_slope, user_data, nbrik,val );
00389 #else
00390       ufunc( tzero,tdelta, nuse,fxar,ts_mean,ts_slope, user_data, nbrik,val );
00391 #endif
00392 
00393       for( iv=0 ; iv < nbrik ; iv++ ) fout[iv][ii] = val[iv] ;
00394 
00395    } /* end of outer loop over 1 voxels at a time */
00396 
00397    DSET_unload( old_dset ) ;  /* don't need this no more */
00398 
00399    /* end notification */
00400 
00401 #if 0
00402    user_func( 0.0 , 0.0 , 0 , NULL,0.0,0.0 , user_data , nbrik , NULL ) ;
00403 #else
00404    ufunc( 0.0 , 0.0 , 0 , NULL,0.0,0.0 , user_data , nbrik , NULL ) ;
00405 #endif
00406 
00407    /*---- Count and correct float errors ----*/
00408 
00409    for( iv=0 ; iv < nbrik ; iv++ )
00410       nbad += thd_floatscan( nvox, fout[iv] ) ;
00411 
00412    if( nbad > 0 )
00413       fprintf(stderr,
00414               "++ Warning: %d bad floats computed in MAKER_4D_to_typed_fbuc\n\a",
00415               nbad ) ;
00416 
00417    /*------------------------------------------------------------------------*/
00418    /*------- The output is now in fout[iv][ii], iv=0..nbrik-1, ii=0..nvox-1.
00419              We must now put this into the output dataset. ------------------*/
00420 
00421    switch( new_datum ){
00422 
00423       /*** output is floats is the simplest:
00424            we just have to attach the fout brick to the dataset ***/
00425 
00426       case MRI_float:
00427          for( iv=0 ; iv < nbrik ; iv++ ){
00428             EDIT_substitute_brick( new_dset , iv , MRI_float , fout[iv] ) ;
00429             fout[iv] = NULL ;  /* so it won't be freed later */
00430          }
00431       break ;
00432 
00433       /*** output is shorts:
00434            we have to create scaled sub-bricks from fout ***/
00435 
00436       case MRI_short:{
00437          short * bout ;
00438          float sfac ;
00439 
00440          for( iv=0 ; iv < nbrik ; iv++ ){
00441             bout = (short *) malloc( sizeof(short) * nvox ) ;
00442             if( bout == NULL ){
00443                fprintf(stderr,
00444                 "\nFinal malloc error in MAKER_4D_to_fbuc - is memory exhausted?\n\a");
00445                EXIT(1) ;
00446             }
00447             sfac = MCW_vol_amax( nvox,1,1 , MRI_float , fout[iv] ) ;
00448             if( sfac > 0.0 ){
00449                sfac = 32767.0 / sfac ;
00450                EDIT_coerce_scale_type( nvox,sfac ,
00451                                        MRI_float,fout[iv] , MRI_short,bout ) ;
00452                sfac = 1.0 / sfac ;
00453             }
00454             val[iv] = sfac ;
00455             EDIT_substitute_brick( new_dset , iv , MRI_short , bout ) ;
00456          }
00457          EDIT_dset_items( new_dset , ADN_brick_fac , val , ADN_none ) ;
00458       }
00459       break ;
00460 
00461       /*** output is bytes (byte = unsigned char)
00462            we have to create a scaled sub-brick from fout ***/
00463 
00464       case MRI_byte:{
00465          byte * bout ;
00466          float sfac ;
00467 
00468          for( iv=0 ; iv < nbrik ; iv++ ){
00469             bout = (byte *) malloc( sizeof(byte) * nvox ) ;
00470             if( bout == NULL ){
00471                fprintf(stderr,
00472                 "\nFinal malloc error in MAKER_4D_to_fbuc - is memory exhausted?\n\a");
00473                EXIT(1) ;
00474             }
00475             sfac = MCW_vol_amax( nvox,1,1 , MRI_float , fout[iv] ) ;
00476             if( sfac > 0.0 ){
00477                sfac = 255.0 / sfac ;
00478                EDIT_coerce_scale_type( nvox,sfac ,
00479                                        MRI_float,fout[iv] , MRI_byte,bout ) ;
00480                sfac = 1.0 / sfac ;
00481             }
00482             val[iv] = sfac ;
00483             EDIT_substitute_brick( new_dset , iv , MRI_byte , bout ) ;
00484          }
00485          EDIT_dset_items( new_dset , ADN_brick_fac , val , ADN_none ) ;
00486       }
00487       break ;
00488 
00489    } /* end of switch on output data type */
00490 
00491    /*-------------- Cleanup and go home ----------------*/
00492 
00493    FREE_WORKSPACE ;
00494    return new_dset ;
00495 }