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