// Declarations of classes that implements useful // utility functions for the eddy current project. // They are collections of statically declared // functions that have been collected into classes // to make it explicit where they come from. There // will never be any instances of theses classes. // // EddyUtils.h // // Jesper Andersson, FMRIB Image Analysis Group // // Copyright (C) 2011 University of Oxford // #ifndef EddyUtils_h #define EddyUtils_h #include #include #include #include #include "newmat.h" #include "newimage/newimageall.h" #include "miscmaths/miscmaths.h" #include "EddyHelperClasses.h" #include "ECScanClasses.h" namespace EDDY { //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ // // Class EddyUtils // // Helper Class used to perform various useful tasks for // the eddy current correction project. // // {{{ @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ class EddyUtils { private: static const int b_range = 100; // b-values within this range considered equivalent static NEWIMAGE::volume4D get_partial_derivatives_in_scan_space(// Input const NEWIMAGE::volume& pred, const EDDY::ECScan& scan, boost::shared_ptr > susc, EDDY::Parameters whichp); static double param_update(// Input const NEWIMAGE::volume& pred, // Prediction in model space boost::shared_ptr > susc, // Susceptibility induced off-resonance field const NEWIMAGE::volume& pmask, EDDY::Parameters whichp, // What parameters do we want to update bool cbs, // Check (success of parameters) Before Set // Input/output EDDY::ECScan& scan, // Scan we want to register to pred // Output NEWMAT::ColumnVector *update); // Vector of updates, optional output static EDDY::ImageCoordinates transform_coordinates_from_model_to_scan_space(// Input const NEWIMAGE::volume& pred, const EDDY::ECScan& scan, boost::shared_ptr > susc, // Output NEWIMAGE::volume *omask, NEWIMAGE::volume *jac); // Returns coordinates into f transformed with // displacement field d and affine matrix M static void transform_coordinates(// Input const NEWIMAGE::volume& f, const NEWIMAGE::volume4D& d, const NEWMAT::Matrix& M, // Output ImageCoordinates& c, NEWIMAGE::volume *omask); // Calculates X.t()*X where X is a matrix where each column is one of the volumes in vols static NEWMAT::Matrix make_XtX(const NEWIMAGE::volume4D& vols, const NEWIMAGE::volume& mask); // Calculates X.t()*y where X is a matrix where each column is one of the volumes in Xvols // and where y is the volume in Yvol. static NEWMAT::ColumnVector make_Xty(const NEWIMAGE::volume4D& Xvols, const NEWIMAGE::volume& Yvol, const NEWIMAGE::volume& mask); public: // This function has been temporarily moved into public space. Should probably be // moved back to private space at some stage. static NEWIMAGE::volume transform_model_to_scan_space(// Input const NEWIMAGE::volume& pred, const EDDY::ECScan& scan, boost::shared_ptr > susc, bool jacmod, // Output NEWIMAGE::volume& omask, NEWIMAGE::volume *jac, NEWIMAGE::volume4D *grad); // Some functions for comparing diffusion parameters static bool AreInSameShell(const DiffPara& dp1, const DiffPara& dp2) { return(fabs(dp1.bVal()-dp2.bVal()) double(b_range)); } static bool Isb0(const DiffPara& dp) { return(!IsDiffusionWeighted(dp)); } static bool HaveSameDirection(const DiffPara& dp1, const DiffPara& dp2) { return(NEWMAT::DotProduct(dp1.bVec(),dp2.bVec())>0.999); } // Random functions to set extrapolation and interpolation // template static void SetTrilinearInterp(V& vol) { if (vol.getinterpolationmethod() != NEWIMAGE::trilinear) vol.setinterpolationmethod(NEWIMAGE::trilinear); if (vol.getextrapolationmethod() != NEWIMAGE::mirror) vol.setextrapolationmethod(NEWIMAGE::mirror); } template static void SetSplineInterp(V& vol) { if (vol.getinterpolationmethod() != NEWIMAGE::spline) vol.setinterpolationmethod(NEWIMAGE::spline); if (vol.getsplineorder() != 3) vol.setsplineorder(3); if (vol.getextrapolationmethod() != NEWIMAGE::mirror) vol.setextrapolationmethod(NEWIMAGE::mirror); } // Check if a pair of ECScans can potentially be used in an LSR reconstruction static bool AreMatchingPair(const ECScan& s1, const ECScan& s2); // Get indicies for non-zero b-values static std::vector GetIndiciesOfDWIs(const std::vector& dpars); // Removes bvecs associated with zero b-values. static std::vector GetDWIDiffParas(const std::vector& dpars); // Reads all diffusion weighted images from 4D volume static int read_DWI_volume4D(NEWIMAGE::volume4D& dwivols, const std::string& fname, const std::vector& dpars); // Converts char mask (from the general_transform functions) to a float mask static NEWIMAGE::volume ConvertMaskToFloat(const NEWIMAGE::volume& charmask); // Calculates slice-wise statistics from the difference between observation and predicton in observation space static DiffStats GetSliceWiseStats(// Input const NEWIMAGE::volume& pred, // Prediction in model space boost::shared_ptr > susc, // Susceptibility induced off-resonance field const NEWIMAGE::volume& pmask, // "Data valid" mask in model space const NEWIMAGE::volume& bmask, // Brain mask in model space EDDY::ECScan& scan); // Scan corresponding to pred // Performs an update of the movement parameters for one scan static double MovParamUpdate(// Input const NEWIMAGE::volume& pred, boost::shared_ptr > susc, const NEWIMAGE::volume& pmask, bool cbs, // Input/output EDDY::ECScan& scan) { return(param_update(pred,susc,pmask,MOVEMENT,cbs,scan,NULL)); } // Performs an update of the EC parameters for one scan static double ECParamUpdate(// Input const NEWIMAGE::volume& pred, boost::shared_ptr > susc, const NEWIMAGE::volume& pmask, bool cbs, // Input/output EDDY::ECScan& scan) { return(param_update(pred,susc,pmask,EC,cbs,scan,NULL)); } // Performs an update of the EC parameters for one scan static double MovAndECParamUpdate(// Input const NEWIMAGE::volume& pred, boost::shared_ptr > susc, const NEWIMAGE::volume& pmask, bool cbs, // Input/output EDDY::ECScan& scan) { return(param_update(pred,susc,pmask,ALL,cbs,scan,NULL)); } // Transforms an image from model/prediction space to observation space static NEWIMAGE::volume TransformModelToScanSpace(// Input const NEWIMAGE::volume& pred, const EDDY::ECScan& scan, boost::shared_ptr > susc) { NEWIMAGE::volume mask(pred.xsize(),pred.ysize(),pred.zsize()); NEWIMAGE::volume jac(pred.xsize(),pred.ysize(),pred.zsize()); return(transform_model_to_scan_space(pred,scan,susc,true,mask,&jac,NULL)); } static NEWIMAGE::volume TransformScanToModelSpace(// Input const EDDY::ECScan& scan, boost::shared_ptr > susc, // Output NEWIMAGE::volume& omask); // The next two are alternate transformation routines that // performs the transforms in several resampling steps. // They are intended for debugging. static NEWIMAGE::volume DirectTransformScanToModelSpace(// Input const EDDY::ECScan& scan, boost::shared_ptr > susc, // Output NEWIMAGE::volume& omask); static NEWIMAGE::volume DirectTransformModelToScanSpace(// Input const NEWIMAGE::volume& ima, const EDDY::ECScan& scan, boost::shared_ptr > susc, // Output NEWIMAGE::volume& omask); // Alternate version of updating routine for writing out debug information. static double param_update_debug(// Input const NEWIMAGE::volume& pred, // Prediction in model space boost::shared_ptr > susc, // Susceptibility induced off-resonance field const NEWIMAGE::volume& pmask, // Parameters whichp, // What parameters do we want to update bool cbs, // Check (success of parameters) Before Set unsigned int scindx, // Scan index unsigned int iter, // Iteration unsigned int level, // Determines how much gets written // Input/output EDDY::ECScan& scan, // Scan we want to register to pred // Output NEWMAT::ColumnVector *rupdate); // Vector of updates, optional output }; // }}} End of fold. //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ // // Class FieldUtils // // Helper Class used to perform various useful calculations // on off-resonance and displacement fields. // //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ class FieldUtils { public: // Rigid body transform of off-resonance field static NEWIMAGE::volume RigidBodyTransformHzField(const NEWIMAGE::volume& hzfield); // Some conversion routines off-resonance->displacements static NEWIMAGE::volume4D Hz2VoxelDisplacements(const NEWIMAGE::volume& hzfield, const AcqPara& acqp); static NEWIMAGE::volume4D Hz2MMDisplacements(const NEWIMAGE::volume& hzfield, const AcqPara& acqp); static NEWIMAGE::volume4D Voxel2MMDisplacements(const NEWIMAGE::volume4D& voxdisp) { NEWIMAGE::volume4D mmd=voxdisp; mmd[0] *= mmd.xdim(); mmd[1] *= mmd.ydim(); mmd[2] *= mmd.zdim(); return(mmd); } // Inverts 3D displacement field, ASSUMING it is really 1D (only non-zero displacements in one direction). static NEWIMAGE::volume4D InvertDisplacementField(// Input const NEWIMAGE::volume4D& dfield, const AcqPara& acqp, const NEWIMAGE::volume& inmask, // Output NEWIMAGE::volume& omask); // Inverts 1D displacement field. Input must be scaled to voxels (i.e. not mm). static NEWIMAGE::volume InvertDisplacementField(// Input const NEWIMAGE::volume& dfield, const AcqPara& acqp, const NEWIMAGE::volume& inmask, // Output NEWIMAGE::volume& omask); // Calculates Jacobian of a displacement field static NEWIMAGE::volume GetJacobian(const NEWIMAGE::volume4D& dfield, const AcqPara& acqp); // Calculates Jacobian of a 1D displacement field static NEWIMAGE::volume GetJacobianFrom1DField(const NEWIMAGE::volume& dfield, unsigned int dir); private: }; } // End namespace EDDY #endif // End #ifndef EddyUtils_h