/*----------------------------------------------------------------------
* Exported functions - these are available for use in other functions
*
* r_new_resam_dset : reorient and/or resamle a dataset
* r_fill_resampled_data_brick : warp a dataset into memory
* r_dxyz_mod_dataxes : modify THD_dataxes struct per dxyz
*----------------------------------------------------------------------
*/
/*----------------------------------------------------------------------
* history:
*
* 2002.07.29
* - if we master, be sure output has same view type
*
* 2003.01.15
* - do not add any history in this library routine
*
* 2004.07.26
* - added sublist parameter to r_new_resam_dset, as a list of
* sub-bricks to resample (all, if NULL)
*
* 2005.03.22
* - removed all tabs
* - must touch file once per year
*
* 2005.08.03
* - allow dxyz to override those from master
*----------------------------------------------------------------------
*/
#include "mrilib.h"
#include "r_misc.h"
#include "r_new_resam_dset.h"
/* local #defines and function definitions */
static int apply_dataxes ( THD_3dim_dataset * dset, THD_dataxes * dax );
static int apply_orientation ( THD_3dim_dataset * dset, FD_brick * fdb,
char orient[] );
static int valid_resam_inputs( THD_3dim_dataset * , THD_3dim_dataset *,
double, double, double, char [], int );
#define LR_RESAM_IN_FAIL -1
#define LR_RESAM_IN_NONE 0
#define LR_RESAM_IN_REORIENT 1
#define LR_RESAM_IN_RESAM 2
static char * this_file = "r_new_resam_dset.c";
/*----------------------------------------------------------------------
* r_new_resam_dset - create a new dataset by resampling an existing one
*
* inputs:
* - din : input dataset
* - min : master dataset - overrides dxyz and orient input
* - dx,dy,dz : new deltas
* - orient : new orientation string
* - resam : resampling mode
* - sublist : list of sub-bricks to resample (all if NULL)
* {first element of sublist is the number of sub-bricks}
*
* output : a pointer to the resulting THD_3dim_dataset
* (or NULL, on failure)
*
* notes:
* - if (min != NULL) , dx, dy, dz and orient are ignored
* - if (dx == dy == dz == 0.0), no resampling will be done
* - if (orient == NULL), no reorienting will be done
*----------------------------------------------------------------------
*/
THD_3dim_dataset * r_new_resam_dset
(
THD_3dim_dataset * din, /* original dataset to resample */
THD_3dim_dataset * min, /* master dataset to align to */
double dx, /* delta for new x-axis */
double dy, /* delta for new y-axis */
double dz, /* delta for new z-axis */
char orient [], /* new orientation code */
int resam, /* mode to resample with */
int * sublist /* list of sub-bricks to resam */
)
{
THD_3dim_dataset * dout;
THD_3dim_dataset * dtmp;
THD_dataxes new_daxes;
FD_brick * fdb = NULL;
int work;
char * l_orient, l_orient_s[4] = "";
/* if we have a master, do not rely on orient - use local memory */
if ( min ) l_orient = l_orient_s;
else l_orient = orient;
work = valid_resam_inputs( din, min, dx, dy, dz, l_orient, resam );
if ( work == LR_RESAM_IN_FAIL )
return NULL;
if ( sublist )
dtmp = THD_copy_dset_subs(din, sublist);
else
dtmp = din;
dout = EDIT_empty_copy( dtmp ); /* create new brick */
if( ! ISVALID_3DIM_DATASET( dout ) )
{
fprintf( stderr, "ERROR: <%s> - failed to duplicate datset at %p\n",
this_file, din );
return NULL;
}
/* give junk name */
EDIT_dset_items(dout, ADN_prefix, "junk_resample", ADN_none );
/* possibly update view type to that of the master */
if ( min &&
(dout->view_type != min->view_type) &&
(dout->dblk && dout->dblk->diskptr) )
{
dout->view_type = min->view_type;
THD_init_diskptr_names( dout->dblk->diskptr, NULL, NULL, NULL,
min->view_type, True );
}
/* re-orient the brick? */
if ( work & LR_RESAM_IN_REORIENT )
{
if ( ( fdb = THD_oriented_brick( dout, l_orient ) ) == NULL )
{
fprintf( stderr, "<%s>: failed to create THD_brick with orient "
"string <%.6s>\n", this_file, l_orient );
THD_delete_3dim_dataset( dout, FALSE );
return NULL;
}
if ( apply_orientation( dout, fdb, l_orient ) != 0 )
{
THD_delete_3dim_dataset( dout, FALSE );
return NULL;
}
}
new_daxes = *dout->daxes; /* in case we don't resample */
/* resample the brick? */
if ( work & LR_RESAM_IN_RESAM )
{
/*-- given dx, dy and dz, create a new THD_dataxes structure ---- */
new_daxes.type = DATAXES_TYPE;
/* fill new_daxes, from the master, new dxyz, or both */
if ( min && dx == 0.0 ) new_daxes = *min->daxes;
else
{
/* possibly start with master (must have master's orient) */
if ( min ) *dout->daxes = *min->daxes;
if ( r_dxyz_mod_dataxes(dx, dy, dz, dout->daxes, &new_daxes) != 0 )
{
THD_delete_3dim_dataset( dout, FALSE );
return NULL;
}
}
if ( apply_dataxes( dout, &new_daxes ) != 0 )
{
THD_delete_3dim_dataset( dout, FALSE );
return NULL;
}
}
/* needed by THD_load_datablock() */
dout->warp_parent = dtmp;
dout->warp_parent_idcode = dtmp->idcode; /* needed for HEAD write */
/* needed to warp from parent */
dout->wod_flag = True; /* mark for WOD */
dout->vox_warp->type = ILLEGAL_TYPE; /* mark fo recomputation */
*dout->wod_daxes = new_daxes; /* used for actual warp */
dout->daxes->parent = (XtPointer)dout; /* parent is new dset */
if ( r_fill_resampled_data_brick( dout, resam ) != 0 )
{
THD_delete_3dim_dataset( dout, FALSE );
dout = NULL;
}
if ( sublist ) DSET_delete(dtmp);
return dout;
}
/*----------------------------------------------------------------------
* r_fill_resampled_data_brick - get the data brick into memory
*
* return FAIL or OK (0)
* - see adwarp_refashion_dataset
*----------------------------------------------------------------------
*/
int r_fill_resampled_data_brick( THD_3dim_dataset * dset, int resam )
{
THD_diskptr * diskptr = dset->dblk->diskptr;
MRI_IMAGE * im;
char * newdata, * dptr;
float bfac;
int ival, dsize;
int nx, ny, nz, nxy, nxyz, nv;
int slice;
if ( DSET_LOADED(dset) )
{
fprintf( stderr, "error <%s>: trying to fill pre-loaded dataset\n",
this_file );
return FAIL;
}
DSET_lock(dset); /* since it will just sit in memory for now */
/* a basic warp is needed if header is written out - PLUTO_add_dset() */
dset->warp = myXtNew( THD_warp );
*dset->warp = IDENTITY_WARP;
ADDTO_KILL( dset->kl, dset->warp );
diskptr->byte_order = mri_short_order();
diskptr->storage_mode = STORAGE_BY_BRICK;
/* note new dimensions */
nx = dset->daxes->nxx; ny = dset->daxes->nyy; nz = dset->daxes->nzz;
nv = diskptr->nvals;
nxy = nx * ny;
nxyz = nx * ny * nz;
/* recompute sub-bricks */
for ( ival = 0; ival < nv; ival++ ) /* for each sub-brick */
{
/* first create memory to deposit the slices into */
dsize = mri_datum_size( DSET_BRICK_TYPE(dset, ival) );
if ( (newdata = (char *)malloc( nxyz * dsize )) == NULL )
{
fprintf( stderr, "r frdb: alloc failure: %d bytes!\n",
nxyz * dsize );
return FAIL;
}
dptr = newdata; /* we will copy images at dptr */
/* force return of unscaled slices for output - reset fac at end */
bfac = DBLK_BRICK_FACTOR(dset->dblk,ival);
DBLK_BRICK_FACTOR(dset->dblk,ival) = 0.0;
/* for each slice, insert it into memory */
for ( slice = 0; slice < nz; slice++)
{
im = AFNI_dataset_slice( dset, 3, slice, ival, resam );
if ( im == NULL )
{
fprintf( stderr, "r_fill_resampled_data_brick: failure to "
"compute dataset slice %d\n", slice );
free( newdata );
return FAIL;
}
memcpy( (void *)dptr, mri_data_pointer(im), nxy*dsize );
mri_free( im );
dptr += nxy*dsize;
}
DBLK_BRICK_FACTOR(dset->dblk,ival) = bfac;
/* we now have the raw brick data, so insert it into the brick */
EDIT_substitute_brick(dset, ival, DSET_BRICK_TYPE(dset,ival),
(void *)newdata);
}
dset->dblk->malloc_type = DATABLOCK_MEM_MALLOC;
dset->wod_flag = False; /* since data is now in memory */
/* recompute statistics */
THD_load_statistics( dset );
return 0; /* OK */
}
/*----------------------------------------------------------------------
* Fill in a THD_dataxes structure to match the 3D box of an existing
* one, but where dx, dy and dz are modified.
*
* This is identical to THD_edit_dataxes(), except that the requirement
* for cubical voxels has been lifted.
*----------------------------------------------------------------------
*/
int r_dxyz_mod_dataxes( double dx, double dy, double dz,
THD_dataxes * daxin, THD_dataxes * daxout
)
{
double rex, rey, rez;
double lxx, lyy, lzz;
int ret_val;
if ( ! ISVALID_DATAXES( daxin ) || ! ISVALID_DATAXES( daxout ) )
return -1;
*daxout = *daxin; /* start with a copy */
if ( dx <= 0.0 || dy <= 0.0 || dz <= 0.0 )
return -1; /* having duplicated the structures */
rex = (daxout->xxdel > 0) ? dx : -dx; /* signed voxel sizes */
rey = (daxout->yydel > 0) ? dy : -dy;
rez = (daxout->zzdel > 0) ? dz : -dz;
lxx = daxin->nxx * daxin->xxdel; /* signed lengths of data box */
lyy = daxin->nyy * daxin->yydel;
lzz = daxin->nzz * daxin->zzdel;
daxout->nxx = (int)( lxx/rex + 0.499 ); /* so this is > 0 */
daxout->nyy = (int)( lyy/rey + 0.499 );
daxout->nzz = (int)( lzz/rez + 0.499 );
/* go from old edge to old center, then back out to new edge */
daxout->xxorg = daxin->xxorg + 0.5*(lxx - daxin->xxdel)
- 0.5*(daxout->nxx - 1)*rex;
daxout->yyorg = daxin->yyorg + 0.5*(lyy - daxin->yydel)
- 0.5*(daxout->nyy - 1)*rey;
daxout->zzorg = daxin->zzorg + 0.5*(lzz - daxin->zzdel)
- 0.5*(daxout->nzz - 1)*rez;
/* dave new dimensions */
daxout->xxdel = rex;
daxout->yydel = rey;
daxout->zzdel = rez;
/* create a new bounding box */
/* (note that xxdel<0 implies we must swap min/max) */
daxout->xxmin = daxout->xxorg;
daxout->xxmax = daxout->xxorg + (daxout->nxx-1)*daxout->xxdel;
if ( daxout->xxmin > daxout->xxmax )
{
double tmp = daxout->xxmin;
daxout->xxmin = daxout->xxmax;
daxout->xxmax = tmp;
}
daxout->yymin = daxout->yyorg;
daxout->yymax = daxout->yyorg + (daxout->nyy-1)*daxout->yydel;
if ( daxout->yymin > daxout->yymax )
{
double tmp = daxout->yymin;
daxout->yymin = daxout->yymax;
daxout->yymax = tmp;
}
daxout->zzmin = daxout->zzorg;
daxout->zzmax = daxout->zzorg + (daxout->nzz-1)*daxout->zzdel;
if ( daxout->zzmin > daxout->zzmax )
{
double tmp = daxout->zzmin;
daxout->zzmin = daxout->zzmax;
daxout->zzmax = tmp;
}
#ifdef EXTEND_BBOX
daxout->xxmin -= 0.5 * daxout->xxdel;
daxout->xxmax += 0.5 * daxout->xxdel;
daxout->yymin -= 0.5 * daxout->yydel;
daxout->yymax += 0.5 * daxout->yydel;
daxout->zzmin -= 0.5 * daxout->zzdel;
daxout->zzmax += 0.5 * daxout->zzdel;
#endif
return ret_val = 0;
}
/*----------------------------------------------------------------------
* validate an orientation string
*
* use explicit boolean return values
*----------------------------------------------------------------------
*/
Boolean r_is_valid_orient_str ( char ostr [] )
{
int o1, o2, o3;
if ( ostr == NULL )
return False;
o1 = ORCODE(toupper(ostr[0]));
o2 = ORCODE(toupper(ostr[1]));
o3 = ORCODE(toupper(ostr[2]));
if ( ( o1 != ILLEGAL_TYPE ) &&
( o2 != ILLEGAL_TYPE ) &&
( o3 != ILLEGAL_TYPE ) &&
OR3OK(o1,o2,o3) )
return True;
else
return False;
}
/*----------------------------------------------------------------------
* given an orientation string, fill an orientation code vector
*
* return 0 on success
*----------------------------------------------------------------------
*/
int r_orient_str2vec ( char ostr [], THD_ivec3 * ovec )
{
int o1, o2, o3;
if ( !ostr || !ovec )
{
fprintf( stderr, "%s: r_orient_str2vec - invalid parameter pair "
"(%p,%p)\n", this_file, ostr, ovec );
return -1;
}
ovec->ijk[0] = o1 = ORCODE(toupper(ostr[0]));
ovec->ijk[1] = o2 = ORCODE(toupper(ostr[1]));
ovec->ijk[2] = o3 = ORCODE(toupper(ostr[2]));
if ( ( o1 == ILLEGAL_TYPE ) ||
( o2 == ILLEGAL_TYPE ) ||
( o3 == ILLEGAL_TYPE ) ||
( !OR3OK(o1,o2,o3) ) )
{
fprintf( stderr, "%s: r_orient_str2vec - bad ostr <%.4s>\n",
this_file, ostr );
return -2;
}
return 0;
}
/*----------------------------------------------------------------------
* validate inputs to driving function
*
* return one of:
*
* LR_RESAM_IN_NONE : no changes from original
* LR_RESAM_IN_REORIENT : do a re-orientation
* LR_RESAM_IN_RESAM : resample dataset
* LR_RESAM_IN_FAIL : bad inputs
*----------------------------------------------------------------------
*/
static int valid_resam_inputs( THD_3dim_dataset * dset, THD_3dim_dataset * mset,
double dx, double dy, double dz,
char orient [], int resam )
{
int ret_val = LR_RESAM_IN_NONE;
if( ! ISVALID_3DIM_DATASET(dset) )
{
fprintf( stderr, "ERROR: <%s> - invalid input dataset\n", this_file );
return LR_RESAM_IN_FAIL;
}
/* validate the resampling mode - this is required in any case */
if ( resam < 0 || resam > LAST_RESAM_TYPE )
{
fprintf( stderr, "ERROR: <%s> - invalid resample mode of %d\n",
this_file, resam );
return LR_RESAM_IN_FAIL;
}
if( mset ) /* behold, the master! */
{
/* validate mset and its thd_dataxes structure */
if ( ! ISVALID_3DIM_DATASET(mset) || !ISVALID_DATAXES(mset->daxes) )
{
fprintf( stderr, "ERROR: <%s> - invalid master dataset\n",
this_file );
return LR_RESAM_IN_FAIL;
}
/* we should have orientation memory with a master */
if ( ! orient )
{
fprintf( stderr, "ERROR: <%s> - orientation memory should "
"come with master\n", this_file );
return LR_RESAM_IN_FAIL;
}
/* we have a master - get orientation now, dxyz will come later */
orient[0] = ORIENT_typestr[mset->daxes->xxorient][0];
orient[1] = ORIENT_typestr[mset->daxes->yyorient][0];
orient[2] = ORIENT_typestr[mset->daxes->zzorient][0];
orient[3] = '\0';
/* validate any new dxyz */
if ( dx != 0.0 && (dx < 0.0 || dy <= 0.0 || dz <= 0.0) )
{
fprintf( stderr, "ERROR: <%s> - invalid (dx,dy,dz) = (%f,%f,%f)\n",
this_file, dx, dy, dz );
return LR_RESAM_IN_FAIL;
}
return LR_RESAM_IN_RESAM | LR_RESAM_IN_REORIENT; /* we're good to go */
}
/* so no master dataset, validate dxyz and orientation */
/* validate resampling */
if ( dx == 0.0 && dy == 0.0 && dz == 0.0 ) /* no resampling */
ret_val &= ~LR_RESAM_IN_RESAM;
else if ( dx <= 0.0 || dy <= 0.0 || dz <= 0.0 )
{
fprintf( stderr, "ERROR: <%s> - invalid (dx,dy,dz) = (%f,%f,%f)\n",
this_file, dx, dy, dz );
return LR_RESAM_IN_FAIL;
}
else
ret_val |= LR_RESAM_IN_RESAM;
/* validate orientation */
if ( orient == NULL )
ret_val &= ~LR_RESAM_IN_REORIENT;
else if ( r_is_valid_orient_str ( orient ) == False )
{
fprintf( stderr, "ERROR: <%s> - invalid orientation string <%.6s>\n",
this_file, orient );
return LR_RESAM_IN_FAIL;
}
else
ret_val |= LR_RESAM_IN_REORIENT;
return ret_val;
}
/*----------------------------------------------------------------------
* Apply orientation and delta values to the dataset.
*
* return 0 on success
*----------------------------------------------------------------------
*/
static int apply_orientation( THD_3dim_dataset * dset, FD_brick * fdb,
char orient[] )
{
THD_dataxes * daxp = dset->daxes;
THD_ivec3 ivnxyz, ivorient;
THD_fvec3 fvdel, fvorg;
float org4[4], del4[4];
int aa1, aa2, aa3;
int a1, a2, a3;
int ret_val;
LOAD_IVEC3(ivnxyz, fdb->n1, fdb->n2, fdb->n3);
if ( (ret_val = r_orient_str2vec(orient, &ivorient)) != 0 )
return ret_val;
LOAD_FVEC3( fvdel,
ORIENT_sign[ivorient.ijk[0]]=='+' ? fdb->del1 : -fdb->del1,
ORIENT_sign[ivorient.ijk[1]]=='+' ? fdb->del2 : -fdb->del2,
ORIENT_sign[ivorient.ijk[2]]=='+' ? fdb->del3 : -fdb->del3 );
/* use original org and delta with permuted index to get new orgs */
org4[1] = daxp->xxorg; org4[2] = daxp->yyorg; org4[3] = daxp->zzorg;
del4[1] = daxp->xxdel; del4[2] = daxp->yydel; del4[3] = daxp->zzdel;
a1 = fdb->a123.ijk[0]; a2 = fdb->a123.ijk[1]; a3 = fdb->a123.ijk[2];
aa1 = abs(a1); aa2 = abs(a2); aa3 = abs(a3);
LOAD_FVEC3( fvorg,
(a1 > 0) ? org4[aa1] : org4[aa1]+(fdb->n1-1)*del4[aa1],
(a2 > 0) ? org4[aa2] : org4[aa2]+(fdb->n2-1)*del4[aa2],
(a3 > 0) ? org4[aa3] : org4[aa3]+(fdb->n3-1)*del4[aa3] );
/* now update the current brick with our new orientation data */
EDIT_dset_items( dset,
ADN_nxyz, ivnxyz,
ADN_xyzdel, fvdel,
ADN_xyzorg, fvorg,
ADN_xyzorient, ivorient,
ADN_none );
return 0;
}
/*----------------------------------------------------------------------
* Apply new dataxes structure to dataset and prepare for warp speed.
*
* return 0 on success
*----------------------------------------------------------------------
*/
static int apply_dataxes( THD_3dim_dataset * dset, THD_dataxes * dax )
{
THD_ivec3 ivnxyz;
THD_fvec3 fvdel, fvorg;
LOAD_IVEC3( ivnxyz, dax->nxx, dax->nyy, dax->nzz );
LOAD_FVEC3( fvorg, dax->xxorg, dax->yyorg, dax->zzorg );
LOAD_FVEC3( fvdel, dax->xxdel, dax->yydel, dax->zzdel );
EDIT_dset_items( dset,
ADN_nxyz, ivnxyz,
ADN_xyzdel, fvdel,
ADN_xyzorg, fvorg,
ADN_none );
/* prepare the dataset for warping */
*dset->wod_daxes = *dax; /* used for actual warp */
*dset->daxes = *dax;
return 0;
}
syntax highlighted by Code2HTML, v. 0.9.1