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    BDW, 24 Feb 1997
00011 
00012    Routine to create a 3D 'fith' dataset from a 3D+time dataset,
00013    using a user supplied function to process each time series.
00014    "user_func" should like
00015 
00016    void user_func( double tzero , double tdelta ,
00017                    int npts , float ts[] , double ts_mean , double ts_slope ,
00018                    void * ud , float * val , float * tval )
00019 
00020    where tzero  =  time at ts[0]
00021          tdelta =  time at ts[1] (i.e., ts[k] is at tzero + k*tdelta);
00022                      tzero and tdelta will be in sec if this is truly "time"
00023          npts   =  number of points in ts array
00024          ts     =  one voxel time series array, ts[0] .. ts[npts-1];
00025                      note that this will always be a float array, and
00026                      that ts will start with the "ignore"-th point of
00027                      the actual voxel time series.
00028         ts_mean =  mean value of ts array
00029        ts_slope =  slope of ts array;
00030                      this will be inversely proportional to tdelta
00031                      (units of 1/sec);
00032                      if "detrend" is nonzero, then the mean and slope
00033                      will been removed from the ts array
00034          ud     =  the user_data pointer passed in here -- this can
00035                      contain whatever control information the user wants
00036          val    =  pointer to return "intensity" value for this voxel;
00037                      note that this must be a float
00038          tval   =  pointer to return "threshold" value for this voxel;
00039                      note that this must be a float between -1 and 1 (inclusive)
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 , NULL , 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 , NULL , 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   Note also that the only difference between the user_func for this function
00053   and the "to_typed_fim" function is the presence of the last argument.
00054 
00055   The inputs to the present routine are
00056     old_dset   = pointer to old dataset;
00057                    note that this dataset must not be warp-on-demand
00058     new_prefix = string to use as filename prefix
00059     new_datum  = type of data to store in output brick;
00060                    if negative, will use datum from old_dset
00061     ignore     = number of data points to ignore at the start
00062     detrend    = if nonzero, this routine will detrend (a+b*t)
00063                    each time series before passing it to user_func
00064     user_func  = discussed above
00065     user_data  = discussed above
00066 
00067   The output consists of pointers to the new dataset.  If NULL is returned,
00068   some error occurred.
00069 ---------------------------------------------------------------------------*/
00070 /* BDW, 24 Feb 1997 */
00071 
00072 /*------------------ macro to return workspace at exit -------------------*/
00073 
00074 #undef  FREE_WORKSPACE
00075 #define FREE_WORKSPACE                              \
00076   do{ FREEUP(bptr) ; FREEUP(sptr) ; FREEUP(fptr) ;  \
00077       FREEUP(cptr) ; FREEUP(fxar) ; FREEUP(fac)  ;  \
00078       FREEUP(fout) ; FREEUP(dtr)  ; FREEUP(tout) ;  \
00079     } while(0)
00080 
00081 /*-------------------------------------------------------------------------*/
00082 /* BDW, 24 Feb 1997 */
00083 
00084 THD_3dim_dataset * MAKER_4D_to_typed_fith( THD_3dim_dataset * old_dset ,
00085                                            char * new_prefix , int new_datum ,
00086                                            int ignore , int detrend ,
00087                                            generic_func * user_func ,
00088                                            void * user_data )
00089 {
00090    THD_3dim_dataset * new_dset ;  /* output dataset */
00091 
00092    byte    ** bptr = NULL ;  /* one of these will be the array of */
00093    short   ** sptr = NULL ;  /* pointers to input dataset sub-bricks */
00094    float   ** fptr = NULL ;  /* (depending on input datum type) */
00095    complex ** cptr = NULL ;
00096 
00097    float * fxar = NULL ;  /* array loaded from input dataset */
00098    float * fac  = NULL ;  /* array of brick scaling factors */
00099    float * fout = NULL ;  /* will be array of output floats (intensity) */
00100    float * tout = NULL ;  /* will be array of output floats (threshold) */
00101    float * dtr  = NULL ;  /* will be array of detrending coeff */
00102 
00103    float val , d0fac , d1fac , x0,x1;
00104    double tzero , tdelta , ts_mean , ts_slope ;
00105    int   ii , old_datum , nuse , use_fac , iz,izold, nxy,nvox , nbad ;
00106    register int kk ;
00107 
00108    void (*ufunc)(double,double,int,float *,double,double,void *,float *,float *)
00109      = (void (*)(double,double,int,float *,double,double,void *,float *,float *))user_func ;
00110 
00111    /*----------------------------------------------------------*/
00112    /*----- Check inputs to see if they are reasonable-ish -----*/
00113 
00114    if( ! ISVALID_3DIM_DATASET(old_dset) ) return NULL ;
00115 
00116    if( new_datum >= 0         &&
00117        new_datum != MRI_short &&
00118        new_datum != MRI_float   ) return NULL ;
00119 
00120    if( user_func == NULL ) return NULL ;
00121 
00122    if( ignore < 0 ) ignore = 0 ;
00123 
00124    /*--------- set up pointers to each sub-brick in the input dataset ---------*/
00125 
00126    old_datum = DSET_BRICK_TYPE( old_dset , 0 ) ;   /* get old dataset datum */
00127    nuse      = DSET_NUM_TIMES(old_dset) - ignore ; /* # of points on time axis */
00128    if( nuse < 2 ) return NULL ;
00129 
00130    if( new_datum < 0 ) new_datum = old_datum ;   /* output datum = input */
00131    if( new_datum == MRI_complex ) return NULL ;  /* but complex = bad news */
00132 
00133    DSET_load( old_dset ) ;  /* must be in memory before we get pointers to it */
00134 
00135    kk = THD_count_databricks( old_dset->dblk ) ;  /* check if it was */
00136    if( kk < DSET_NVALS(old_dset) ){               /* loaded correctly */
00137       DSET_unload( old_dset ) ;
00138       return NULL ;
00139    }
00140 
00141    switch( old_datum ){  /* pointer type depends on input datum type */
00142 
00143       default:                      /** don't know what to do **/
00144          DSET_unload( old_dset ) ;
00145          return NULL ;
00146 
00147       /** create array of pointers into old dataset sub-bricks **/
00148 
00149       /*--------- input is bytes ----------*/
00150       /* voxel #i at time #k is bptr[k][i] */
00151       /* for i=0..nvox-1 and k=0..nuse-1.  */
00152 
00153       case MRI_byte:
00154          bptr = (byte **) malloc( sizeof(byte *) * nuse ) ;
00155          if( bptr == NULL ) return NULL ;
00156          for( kk=0 ; kk < nuse ; kk++ )
00157             bptr[kk] = (byte *) DSET_ARRAY(old_dset,kk+ignore) ;
00158       break ;
00159 
00160       /*--------- input is shorts ---------*/
00161       /* voxel #i at time #k is sptr[k][i] */
00162       /* for i=0..nvox-1 and k=0..nuse-1.  */
00163 
00164       case MRI_short:
00165          sptr = (short **) malloc( sizeof(short *) * nuse ) ;
00166          if( sptr == NULL ) return NULL ;
00167          for( kk=0 ; kk < nuse ; kk++ )
00168             sptr[kk] = (short *) DSET_ARRAY(old_dset,kk+ignore) ;
00169       break ;
00170 
00171       /*--------- input is floats ---------*/
00172       /* voxel #i at time #k is fptr[k][i] */
00173       /* for i=0..nvox-1 and k=0..nuse-1.  */
00174 
00175       case MRI_float:
00176          fptr = (float **) malloc( sizeof(float *) * nuse ) ;
00177          if( fptr == NULL ) return NULL ;
00178          for( kk=0 ; kk < nuse ; kk++ )
00179             fptr[kk] = (float *) DSET_ARRAY(old_dset,kk+ignore) ;
00180       break ;
00181 
00182       /*--------- input is complex ---------*/
00183       /* voxel #i at time #k is cptr[k][i]  */
00184       /* for i=0..nvox-1 and k=0..nuse-1.   */
00185 
00186       case MRI_complex:
00187          cptr = (complex **) malloc( sizeof(complex *) * nuse ) ;
00188          if( cptr == NULL ) return NULL ;
00189          for( kk=0 ; kk < nuse ; kk++ )
00190             cptr[kk] = (complex *) DSET_ARRAY(old_dset,kk+ignore) ;
00191       break ;
00192 
00193    } /* end of switch on input type */
00194 
00195    /*---- allocate space for 1 voxel timeseries ----*/
00196 
00197    fxar = (float *) malloc( sizeof(float) * nuse ) ;   /* voxel timeseries */
00198    if( fxar == NULL ){ FREE_WORKSPACE ; return NULL ; }
00199 
00200    /*--- get scaling factors for sub-bricks ---*/
00201 
00202    fac = (float *) malloc( sizeof(float) * nuse ) ;   /* factors */
00203    if( fac == NULL ){ FREE_WORKSPACE ; return NULL ; }
00204 
00205    use_fac = 0 ;
00206    for( kk=0 ; kk < nuse ; kk++ ){
00207       fac[kk] = DSET_BRICK_FACTOR(old_dset,kk+ignore) ;
00208       if( fac[kk] != 0.0 ) use_fac++ ;
00209       else                 fac[kk] = 1.0 ;
00210    }
00211    if( !use_fac ) FREEUP(fac) ;
00212 
00213    /*--- setup for detrending ---*/
00214 
00215    dtr = (float *) malloc( sizeof(float) * nuse ) ;
00216    if( dtr == NULL ){ FREE_WORKSPACE ; return NULL ; }
00217 
00218    d0fac = 1.0 / nuse ;
00219    d1fac = 12.0 / nuse / (nuse*nuse - 1.0) ;
00220    for( kk=0 ; kk < nuse ; kk++ )
00221       dtr[kk] = kk - 0.5 * (nuse-1) ;  /* linear trend, orthogonal to 1 */
00222 
00223    /*---------------------- make a new dataset ----------------------*/
00224 
00225    new_dset = EDIT_empty_copy( old_dset ) ; /* start with copy of old one */
00226 
00227    /*-- edit some of its internal parameters --*/
00228 
00229    ii = EDIT_dset_items(
00230            new_dset ,
00231               ADN_prefix      , new_prefix ,           /* filename prefix */
00232               ADN_malloc_type , DATABLOCK_MEM_MALLOC , /* store in memory */
00233               ADN_datum_all   , new_datum ,            /* atomic datum */
00234               ADN_nvals       , 2 ,                    /* # sub-bricks */
00235               ADN_ntt         , 0 ,                    /* # time points */
00236               ADN_type        , ISHEAD(old_dset)       /* dataset type */
00237                                  ? HEAD_FUNC_TYPE
00238                                  : GEN_FUNC_TYPE ,
00239               ADN_func_type   , FUNC_THR_TYPE ,        /* function type */
00240            ADN_none ) ;
00241 
00242    if( ii != 0 ){
00243       THD_delete_3dim_dataset( new_dset , False ) ;  /* some error above */
00244       FREE_WORKSPACE ; return NULL ;
00245    }
00246 
00247    /*------ make floating point output brick
00248             (only at the end will scale to byte or shorts) ------*/
00249 
00250    nvox = old_dset->daxes->nxx * old_dset->daxes->nyy * old_dset->daxes->nzz ;
00251 
00252    fout = (float *) malloc( sizeof(float) * nvox ); /* ptr to intensity brick */
00253    tout = (float *) malloc( sizeof(float) * nvox ); /* ptr to threshold brick */
00254 
00255    if( ( fout == NULL ) || (tout == NULL) ){
00256       THD_delete_3dim_dataset( new_dset , False ) ;
00257       FREE_WORKSPACE ; return NULL ;
00258    }
00259 
00260    /*----- set up to find time at each voxel -----*/
00261 
00262    tdelta = old_dset->taxis->ttdel ;
00263    if( DSET_TIMEUNITS(old_dset) == UNITS_MSEC_TYPE ) tdelta *= 0.001 ;
00264    if( tdelta == 0.0 ) tdelta = 1.0 ;
00265 
00266    izold  = -666 ;
00267    nxy    = old_dset->daxes->nxx * old_dset->daxes->nyy ;
00268 
00269    /*----------------------------------------------------*/
00270    /*----- Setup has ended.  Now do some real work. -----*/
00271 
00272    /* start notification */
00273 
00274 #if 0
00275    user_func(  0.0 , 0.0 , nvox , NULL,0.0,0.0 , user_data , NULL , NULL ) ;
00276 #else
00277    ufunc(  0.0 , 0.0 , nvox , NULL,0.0,0.0 , user_data , NULL , NULL ) ;
00278 #endif
00279 
00280    /***** loop over voxels *****/
00281 
00282    for( ii=0 ; ii < nvox ; ii++  ){  /* 1 time series at a time */
00283 
00284       /*** load data from input dataset, depending on type ***/
00285 
00286       switch( old_datum ){
00287 
00288          /*** input = bytes ***/
00289 
00290          case MRI_byte:
00291             for( kk=0 ; kk < nuse ; kk++ ) fxar[kk] = bptr[kk][ii] ;
00292          break ;
00293 
00294          /*** input = shorts ***/
00295 
00296          case MRI_short:
00297             for( kk=0 ; kk < nuse ; kk++ ) fxar[kk] = sptr[kk][ii] ;
00298          break ;
00299 
00300          /*** input = floats ***/
00301 
00302          case MRI_float:
00303             for( kk=0 ; kk < nuse ; kk++ ) fxar[kk] = fptr[kk][ii] ;
00304          break ;
00305 
00306          /*** input = complex (note we use absolute value) ***/
00307 
00308          case MRI_complex:
00309             for( kk=0 ; kk < nuse ; kk++ ) fxar[kk] = CABS(cptr[kk][ii]) ;
00310          break ;
00311 
00312       } /* end of switch over input type */
00313 
00314       /*** scale? ***/
00315 
00316       if( use_fac )
00317          for( kk=0 ; kk < nuse ; kk++ ) fxar[kk] *= fac[kk] ;
00318 
00319       /** compute mean and slope **/
00320 
00321       x0 = x1 = 0.0 ;
00322       for( kk=0 ; kk < nuse ; kk++ ){
00323          x0 += fxar[kk] ; x1 += fxar[kk] * dtr[kk] ;
00324       }
00325 
00326       x0 *= d0fac ; x1 *= d1fac ;  /* factors to remove mean and trend */
00327 
00328       ts_mean  = x0 ;
00329       ts_slope = x1 / tdelta ;
00330 
00331       /** detrend? **/
00332 
00333       if( detrend )
00334          for( kk=0 ; kk < nuse ; kk++ ) fxar[kk] -= (x0 + x1 * dtr[kk]) ;
00335 
00336       /** compute start time of this timeseries **/
00337 
00338       iz = ii / nxy ;    /* which slice am I in? */
00339 
00340       if( iz != izold ){          /* in a new slice? */
00341          tzero = THD_timeof( ignore ,
00342                              old_dset->daxes->zzorg
00343                            + iz*old_dset->daxes->zzdel , old_dset->taxis ) ;
00344          izold = iz ;
00345 
00346          if( DSET_TIMEUNITS(old_dset) == UNITS_MSEC_TYPE ) tzero *= 0.001 ;
00347       }
00348 
00349       /*** compute output ***/
00350 
00351 #if 0
00352       user_func( tzero,tdelta , nuse,fxar,ts_mean,ts_slope , user_data ,
00353                  fout+ii , tout+ii ) ;
00354 #else
00355       ufunc( tzero,tdelta , nuse,fxar,ts_mean,ts_slope , user_data ,
00356                  fout+ii , tout+ii ) ;
00357 #endif
00358 
00359       /***  limit threshold data to [-1,+1]  ***/
00360       if (tout[ii] >  1.0)  tout[ii] =  1.0;
00361       if (tout[ii] < -1.0)  tout[ii] = -1.0;
00362 
00363    } /* end of outer loop over 1 voxels at a time */
00364 
00365    DSET_unload( old_dset ) ;  /* don't need this no more */
00366 
00367    /* end notification */
00368 
00369 #if 0
00370    user_func( 0.0 , 0.0 , 0 , NULL,0.0,0.0 , user_data , NULL , NULL ) ;
00371 #else
00372    ufunc( 0.0 , 0.0 , 0 , NULL,0.0,0.0 , user_data , NULL , NULL ) ;
00373 #endif
00374 
00375    nbad = thd_floatscan(nvox,fout) + thd_floatscan(nvox,tout) ;
00376    if( nbad > 0 )
00377       fprintf(stderr,
00378               "++ Warning: %d bad floats computed in MAKER_4D_to_typed_fith\n\a",
00379               nbad ) ;
00380 
00381    /*-------------------------------------------------------------------*/
00382    /*------- The output is now in fout[ii] and tout[ii],  ii=0..nvox-1.
00383              We must now put this into the output dataset --------------*/
00384 
00385    switch( new_datum ){
00386 
00387       /*** output is floats is the simplest:
00388            we just have to attach the fout and tout bricks to the dataset ***/
00389 
00390       case MRI_float:
00391          EDIT_substitute_brick( new_dset , 0 , MRI_float , fout ) ;
00392          EDIT_substitute_brick( new_dset , 1 , MRI_float , tout ) ;
00393          fout = NULL ;   tout = NULL;    /* so it won't be freed later */
00394       break ;
00395 
00396       /*** output is shorts:
00397            we have to create scaled sub-bricks from fout and tout ***/
00398 
00399       case MRI_short:{
00400          short * bout ;
00401          short * tbout ;
00402          float sfac[2] ;
00403 
00404          /*-- get output sub-bricks --*/
00405 
00406          bout = (short *) malloc( sizeof(short) * nvox ) ;
00407          tbout = (short *) malloc( sizeof(short) * nvox ) ;
00408          if( (bout == NULL) || (tbout == NULL) ){
00409             fprintf(stderr,
00410              "\nFinal malloc error in MAKER_4D_to_fith - is memory exhausted?\n\a") ;
00411             EXIT(1) ;
00412          }
00413 
00414          /*-- find scaling and then scale --*/
00415 
00416          sfac[0] = MCW_vol_amax( nvox,1,1 , MRI_float , fout ) ;
00417          if( sfac[0] > 0.0 ){
00418             sfac[0] = 32767.0 / sfac[0] ;
00419             EDIT_coerce_scale_type( nvox,sfac[0] ,
00420                                     MRI_float,fout , MRI_short,bout ) ;
00421             sfac[0] = 1.0 / sfac[0] ;
00422          }
00423          sfac[1] = FUNC_THR_SCALE_SHORT;
00424          EDIT_coerce_scale_type( nvox,sfac[1] ,
00425                                  MRI_float,tout , MRI_short,tbout ) ;
00426          sfac[1] = 1.0 / sfac[1];
00427 
00428          /*-- put output bricks into dataset, and store scale factor --*/
00429 
00430          EDIT_substitute_brick( new_dset , 0 , MRI_short , bout ) ;
00431          EDIT_substitute_brick( new_dset , 1 , MRI_short , tbout ) ;
00432          EDIT_dset_items( new_dset , ADN_brick_fac , sfac , ADN_none ) ;
00433       }
00434       break ;
00435 
00436 
00437    } /* end of switch on output data type */
00438 
00439    /*-------------- Cleanup and go home ----------------*/
00440 
00441    FREE_WORKSPACE ;
00442    return new_dset ;
00443 }