00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006    
00007 #include "thd_maker.h"
00008 
00009 /*-------------------------------------------------------------------------
00010    Routine to create a 3D 'fbuc' dataset from a 3D+time dataset,
00011    using a user supplied function to process each time series.
00012    "user_func" should like
00013 
00014    void user_func( double tzero , double tdelta ,
00015                    int npts , float ts[] , double ts_mean , double ts_slope ,
00016                    void * ud , int nbrik , float * val )
00017 
00018    where tzero  =  time at ts[0]
00019          tdelta =  time at ts[1] (i.e., ts[k] is at tzero + k*tdelta);
00020                      tzero and tdelta will be in sec if this is truly "time"
00021          npts   =  number of points in ts array
00022          ts     =  one voxel time series array, ts[0] .. ts[npts-1];
00023                      note that this will always be a float array, and
00024                      that ts will start with the "ignore"-th point of
00025                      the actual voxel time series.
00026         ts_mean =  mean value of ts array
00027        ts_slope =  slope of ts array;
00028                      this will be inversely proportional to tdelta
00029                      (units of 1/sec);
00030                      if "detrend" is nonzero, then the mean and slope
00031                      will been removed from the ts array
00032          ud     =  the user_data pointer passed in here -- this can
00033                      contain whatever control information the user wants
00034          nbrik  =  number of output values that this function should return
00035                      for the voxel corresponding to input data 'ts'
00036          val    =  pointer to return values for this voxel;
00037                      note that this is a float array of length nbrik,
00038                      and that values that you don't fill in will be
00039                      set to zero (don't overrun this array!).
00040 
00041   Before the first timeseries is passed, user_func will be called with
00042   arguments
00043      ( 0.0 , 0.0 , nvox , NULL , 0.0 , 0.0 , user_data , nbrik , NULL )
00044   where nvox = number of voxels that will be processed.
00045   This is to allow for some setup (e.g., malloc, PLUTO_popup_meter, ...).
00046 
00047   After the last timeseries is passed, user_func will be called again with
00048   arguments
00049      ( 0.0 , 0.0 , 0 , NULL , 0.0 , 0.0 , user_data , nbrik , NULL )
00050   This is to allow for cleanup (e.g., free of malloc, ...).  Note that the
00051   only difference between these "notification" calls is the third argument.
00052 
00053   The inputs to the present routine are
00054     old_dset   = pointer to old dataset;
00055                    note that this dataset must not be warp-on-demand
00056     new_prefix = string to use as filename prefix
00057     new_datum  = type of data to store in output brick;
00058                    if negative, will use datum from old_dset
00059     ignore     = number of data points to ignore at the start
00060     detrend    = if nonzero, this routine will detrend (a+b*t)
00061                    each time series before passing it to user_func
00062     nbrik      = number of values (and sub-bricks) to create at each
00063                    voxel location
00064     user_func  = discussed above
00065     user_data  = discussed above
00066 
00067   The output is a pointer to a new dataset.  If NULL is returned,
00068   some error occurred.
00069 
00070   Note that sub-brick labels, statistical data, etc., must be added to
00071   the new dataset (using EDIT_dset_items) after this routine returns.
00072 ---------------------------------------------------------------------------*/
00073 
00074 /*------------------ macros to return workspace at exit -------------------*/
00075 
00076 #undef  FREE_FOUT
00077 #define FREE_FOUT                                            \
00078   do{ int jv ;                                               \
00079       if( fout != NULL ){                                    \
00080          for( jv=0 ; jv < nbrik ; jv++ ) FREEUP(fout[jv]) ;  \
00081          free(fout) ;                                        \
00082       } } while(0)
00083 
00084 #undef  FREE_WORKSPACE
00085 #define FREE_WORKSPACE                              \
00086   do{ FREEUP(bptr) ; FREEUP(sptr) ; FREEUP(fptr) ;  \
00087       FREEUP(cptr) ; FREEUP(fxar) ; FREEUP(fac)  ;  \
00088       FREEUP(dtr)  ; FREEUP(val)  ; FREE_FOUT    ;  \
00089     } while(0)
00090 
00091 /*-------------------------------------------------------------------------*/
00092 
00093 THD_3dim_dataset * MAKER_4D_to_typed_fbuc( THD_3dim_dataset * old_dset ,
00094                                            char * new_prefix , int new_datum ,
00095                                            int ignore , int detrend ,
00096                                            int nbrik , generic_func * user_func ,
00097                                            void * user_data )
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 }