/*****************************************************************************
Major portions of this software are copyrighted by the Medical College
of Wisconsin, 1994-2000, and are released under the Gnu General Public
License, Version 2. See the file README.Copyright for details.
******************************************************************************/
#include "mrilib.h"
#include "thd.h"
/*---------------------------------------------------------------
Routine to extract a time-series (fixed index, variable ival)
from a dataset.
ixyz = spatial index of desired voxel
= ix + jy * n1 + kz * n1*n2
raw = 0 if you always want floats
= 1 if you want the truly stored data type
05 Nov 2001: Split into 2 functions -
THD_extract_series() produces a new image each time
THD_extract_array() copies data into a user-supplied array
-----------------------------------------------------------------*/
MRI_IMAGE * THD_extract_series( int ind , THD_3dim_dataset *dset , int raw )
{
int nv , typ , ii ;
MRI_IMAGE *im ;
void *imar ;
ENTRY("THD_extract_series") ;
if( !ISVALID_DSET(dset) ) RETURN(NULL) ;
nv = dset->dblk->nvals ;
if( raw ) typ = DSET_BRICK_TYPE(dset,0) ; /* type of output array */
else typ = MRI_float ;
im = mri_new( nv , 1 , typ ) ; /* output image */
imar = mri_data_pointer(im) ;
ii = THD_extract_array( ind , dset , raw , imar ) ; /* get data */
if( ii != 0 ){ mri_free(im) ; RETURN(NULL) ; } /* bad */
if( dset->taxis != NULL ){ /* 21 Oct 1996 */
float zz , tt ;
int kz = ind / ( dset->daxes->nxx * dset->daxes->nyy ) ;
zz = dset->daxes->zzorg + kz * dset->daxes->zzdel ;
tt = THD_timeof( 0 , zz , dset->taxis ) ;
im->xo = tt ; im->dx = dset->taxis->ttdel ; /* origin and delta */
if( dset->taxis->units_type == UNITS_MSEC_TYPE ){ /* convert to sec */
im->xo *= 0.001 ; im->dx *= 0.001 ;
}
} else {
im->xo = 0.0 ; im->dx = 1.0 ; /* 08 Nov 1996 */
}
RETURN(im) ;
}
/*---------------------------------------------------------------------------
Return value is 0 for all is good, -1 for all is bad.
If raw == 0, uar is of type float.
If raw != 0, uar is of the same type as the dataset brick
(the user need to know the type ahead of time if raw != 0).
Data goes into a user-supplied array.
-----------------------------------------------------------------------------*/
int THD_extract_array( int ind, THD_3dim_dataset *dset, int raw, void *uar )
{
MRI_TYPE typ ;
int nv , ival , nb , nb1 ;
char *iar ; /* brick in the input */
float *far=NULL ; /* non-raw output */
static void *tar=NULL ; static int ntar=0 ;
ENTRY("THD_extract_array") ;
if( ind < 0 || uar == NULL ||
!ISVALID_DSET(dset) || ind >= DSET_NVOX(dset) ) RETURN(-1) ;
nv = dset->dblk->nvals ;
iar = DSET_ARRAY(dset,0) ;
if( iar == NULL ){ /* load data from disk? */
DSET_load(dset) ;
iar = DSET_ARRAY(dset,0); if( iar == NULL ) RETURN(-1) ;
}
typ = DSET_BRICK_TYPE(dset,0) ; /* raw data type */
/* will extract nb bytes of raw data into array tar */
nb1 = mri_datum_size(typ); nb = nb1 * (nv+1); nb1 = nb1 * nv;
if( nb > ntar ){ tar = AFREALL(tar,void *,nb) ; ntar = nb ; }
memset(tar,0,nb) ;
if( !raw ) far = (float *) uar ; /* non-raw output */
switch( typ ){
default: /* don't know what to do --> return nada */
RETURN(-1);
break ;
case MRI_byte:{
byte *ar = (byte *)tar , *bar ;
for( ival=0 ; ival < nv ; ival++ ){
bar = (byte *) DSET_ARRAY(dset,ival) ;
if( bar != NULL ) ar[ival] = bar[ind] ;
}
if( !raw ){
for( ival=0 ; ival < nv ; ival++ ) far[ival] = ar[ival] ;
}
}
break ;
case MRI_short:{
short *ar = (short *)tar , *bar ;
for( ival=0 ; ival < nv ; ival++ ){
bar = (short *) DSET_ARRAY(dset,ival) ;
if( bar != NULL ) ar[ival] = bar[ind] ;
}
if( !raw ){
for( ival=0 ; ival < nv ; ival++ ) far[ival] = ar[ival] ;
}
}
break ;
case MRI_float:{
float *ar = (float *)tar , *bar ;
for( ival=0 ; ival < nv ; ival++ ){
bar = (float *) DSET_ARRAY(dset,ival) ;
if( bar != NULL ) ar[ival] = bar[ind] ;
}
if( !raw ){
for( ival=0 ; ival < nv ; ival++ ) far[ival] = ar[ival] ;
}
}
break ;
case MRI_int:{
int *ar = (int *)tar , *bar ;
for( ival=0 ; ival < nv ; ival++ ){
bar = (int *) DSET_ARRAY(dset,ival) ;
if( bar != NULL ) ar[ival] = bar[ind] ;
}
if( !raw ){
for( ival=0 ; ival < nv ; ival++ ) far[ival] = ar[ival] ;
}
}
break ;
case MRI_double:{
double *ar = (double *)tar , *bar ;
for( ival=0 ; ival < nv ; ival++ ){
bar = (double *) DSET_ARRAY(dset,ival) ;
if( bar != NULL ) ar[ival] = bar[ind] ;
}
if( !raw ){
for( ival=0 ; ival < nv ; ival++ ) far[ival] = ar[ival] ;
}
}
break ;
case MRI_complex:{
complex *ar = (complex *)tar , *bar ;
for( ival=0 ; ival < nv ; ival++ ){
bar = (complex *) DSET_ARRAY(dset,ival) ;
if( bar != NULL ) ar[ival] = bar[ind] ;
}
if( !raw ){
for( ival=0 ; ival < nv ; ival++ ) far[ival] = CABS(ar[ival]) ;
}
}
break ;
}
if( raw ){ memcpy(uar,tar,nb1); RETURN(0); }
if( THD_need_brick_factor(dset) ){
for( ival=0 ; ival < nv ; ival++ )
if( DSET_BRICK_FACTOR(dset,ival) > 0.0 )
far[ival] *= DSET_BRICK_FACTOR(dset,ival) ;
}
RETURN(0);
}
/*----------------------------------------------------------------------------*/
int THD_voxel_is_constant( int ind , THD_3dim_dataset *dset )
{
float *far ; int ii,nvox,nvals ;
if( !ISVALID_DSET(dset) ) return 1 ;
if( ind < 0 || ind >= DSET_NVOX(dset) ) return 1 ;
nvals = DSET_NVALS(dset) ; if( nvals == 1 ) return 1 ;
far = (float *)malloc(sizeof(float)*nvals) ;
ii = THD_extract_array( ind , dset , 0 , far ) ;
if( ii < 0 ){ free(far); return 1; }
for( ii=1 ; ii < nvals && far[ii]==far[0]; ii++ ) ; /*nada*/
free(far) ; return (ii==nvals) ;
}
/*----------------------------------------------------------------------------
04 Feb 2000: do a bunch of timeseries at once (for efficiency)
27 Feb 2003: rearranged slightly for more efficiency
------------------------------------------------------------------------------*/
MRI_IMARR * THD_extract_many_series( int ns, int *ind, THD_3dim_dataset *dset )
{
MRI_IMARR *imar ; /* output */
MRI_IMAGE *im ;
int nv , ival , kk ;
char *iar ; /* brick in the input */
float **far ; /* 27 Feb 2003: ptrs to output */
ENTRY("THD_extract_many_series") ;
if( ns <= 0 || ind == NULL | dset == NULL ) RETURN( NULL );
/* try to load dataset */
nv = dset->dblk->nvals ;
iar = DSET_ARRAY(dset,0) ;
if( iar == NULL ){ /* if data needs to be loaded from disk */
(void) THD_load_datablock( dset->dblk ) ;
iar = DSET_ARRAY(dset,0) ;
if( iar == NULL ) RETURN( NULL );
}
/* create output */
far = (float **) malloc(sizeof(float *)*ns) ; /* 27 Feb 2003 */
INIT_IMARR(imar) ;
for( kk=0 ; kk < ns ; kk++ ){
im = mri_new( nv , 1 , MRI_float ) ; /* N.B.: now does 0 fill */
far[kk] = MRI_FLOAT_PTR(im) ; /* ptr to kk-th output series */
ADDTO_IMARR(imar,im) ;
}
/* fill the output */
switch( DSET_BRICK_TYPE(dset,0) ){
default: /* don't know what to do --> return nada */
DESTROY_IMARR(imar) ; free(far) ;
RETURN( NULL );
case MRI_byte:{
byte * bar ;
for( ival=0 ; ival < nv ; ival++ ){
bar = (byte *) DSET_ARRAY(dset,ival) ;
if( bar != NULL ){
for( kk=0 ; kk < ns ; kk++ ){
far[kk][ival] = (float)bar[ind[kk]] ;
}
}
}
}
break ;
case MRI_short:{
short * bar ;
for( ival=0 ; ival < nv ; ival++ ){
bar = (short *) DSET_ARRAY(dset,ival) ;
if( bar != NULL ){
for( kk=0 ; kk < ns ; kk++ ){
far[kk][ival] = (float)bar[ind[kk]] ;
}
}
}
}
break ;
case MRI_float:{
float * bar ;
for( ival=0 ; ival < nv ; ival++ ){
bar = (float *) DSET_ARRAY(dset,ival) ;
if( bar != NULL ){
for( kk=0 ; kk < ns ; kk++ ){
far[kk][ival] = bar[ind[kk]] ;
}
}
}
}
break ;
case MRI_int:{
int * bar ;
for( ival=0 ; ival < nv ; ival++ ){
bar = (int *) DSET_ARRAY(dset,ival) ;
if( bar != NULL ){
for( kk=0 ; kk < ns ; kk++ ){
far[kk][ival] = bar[ind[kk]] ;
}
}
}
}
break ;
case MRI_double:{
double * bar ;
for( ival=0 ; ival < nv ; ival++ ){
bar = (double *) DSET_ARRAY(dset,ival) ;
if( bar != NULL ){
for( kk=0 ; kk < ns ; kk++ ){
far[kk][ival] = (float)bar[ind[kk]] ;
}
}
}
}
break ;
case MRI_complex:{
complex * bar ;
for( ival=0 ; ival < nv ; ival++ ){
bar = (complex *) DSET_ARRAY(dset,ival) ;
if( bar != NULL ){
for( kk=0 ; kk < ns ; kk++ ){
far[kk][ival] = bar[ind[kk]].r ;
}
}
}
}
break ;
}
/* scale outputs, if needed */
if( THD_need_brick_factor(dset) ){
MRI_IMAGE * qim ;
for( kk=0 ; kk < ns ; kk++ ){
im = IMARR_SUBIMAGE(imar,kk) ;
qim = mri_mult_to_float( dset->dblk->brick_fac , im ) ;
mri_free(im) ;
IMARR_SUBIMAGE(imar,kk) = qim ;
}
}
#if 0 /* 27 Feb 2003 */
/* convert to floats, if needed */
if( IMARR_SUBIMAGE(imar,0)->kind != MRI_float ){
MRI_IMAGE * qim ;
for( kk=0 ; kk < ns ; kk++ ){
im = IMARR_SUBIMAGE(imar,kk) ;
qim = mri_to_float( im ) ;
mri_free(im) ;
IMARR_SUBIMAGE(imar,kk) = qim ;
}
}
#endif
/* add time axis stuff to output images, if present */
if( dset->taxis != NULL ){
float zz , tt ;
int kz ;
for( kk=0 ; kk < ns ; kk++ ){
kz = ind[kk] / ( dset->daxes->nxx * dset->daxes->nyy ) ;
zz = dset->daxes->zzorg + kz * dset->daxes->zzdel ;
tt = THD_timeof( 0 , zz , dset->taxis ) ;
im = IMARR_SUBIMAGE(imar,kk) ;
im->xo = tt ; im->dx = dset->taxis->ttdel ; /* origin and delta */
if( dset->taxis->units_type == UNITS_MSEC_TYPE ){ /* convert to sec */
im->xo *= 0.001 ; im->dx *= 0.001 ;
}
}
} else {
for( kk=0 ; kk < ns ; kk++ ){
im = IMARR_SUBIMAGE(imar,kk) ;
im->xo = 0.0 ; im->dx = 1.0 ;
}
}
free(far) ; RETURN(imar);
}
syntax highlighted by Code2HTML, v. 0.9.1