/*---------------------------------------------------------------------- * 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 *---------------------------------------------------------------------- */ #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; if ( min ) /* if we have a master, just copy the dataxes struct */ new_daxes = *min->daxes; else 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; } 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'; 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 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; } /* 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; }