// 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 #include "armawrap/newmat.h" #include "newimage/newimageall.h" #include "miscmaths/miscmaths.h" #include "EddyHelperClasses.h" #include "ECScanClasses.h" #include "EddyCommandLineOptions.h" namespace EDDY { //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ // // Class EddyUtils // // Helper Class used to perform various useful tasks for // the eddy current correction project. // //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ class EddyUtils { private: /// b-values within this range considered equivalent static double b_range; /// Uncomment next line and delete line above when we move NVCC to C++17 // inline static double b_range = 100; /// bladibla static NEWIMAGE::volume4D get_partial_derivatives_in_scan_space(// Input const NEWIMAGE::volume& pred, // Prediction in model space const EDDY::ECScan& scan, // Scan space std::shared_ptr > susc, // Susceptibility off-resonance field EDDY::ParametersType whichp); static NEWIMAGE::volume4D get_direct_partial_derivatives_in_scan_space(// Input const NEWIMAGE::volume& pred, // Prediction in model space const EDDY::ECScan& scan, // Scan space std::shared_ptr > susc, // Susceptibility off-resonance field EDDY::ParametersType whichp); static NEWIMAGE::volume4D get_long_ec_partial_derivatives_in_scan_space(// Input const NEWIMAGE::volume& pred, // Prediction in model space const EDDY::ECScan& scan, // Scan space unsigned int si, // Scan index const EDDY::LongECModel& lecm, // Model for long EC std::shared_ptr > susc); // Susceptibility off-resonance field static double param_update(// Input const NEWIMAGE::volume& pred, // Prediction in model space std::shared_ptr > susc, // Susceptibility off-resonance field std::shared_ptr > bias, // Recieve bias field const NEWIMAGE::volume& pmask, // "Data valid" mask in model space EDDY::ParametersType whichp, // Which parameters to update float fwhm, // FWHM for Gaussian smoothing */ bool very_verbose, // These input parameters are for debugging only 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 static void write_debug_info_for_param_update(const EDDY::ECScan& scan, unsigned int scindx, unsigned int iter, unsigned int level, float fwhm, const NEWIMAGE::volume4D& derivs, const NEWIMAGE::volume& mask, const NEWIMAGE::volume& mios, const NEWIMAGE::volume& pios, const NEWIMAGE::volume& jac, std::shared_ptr > susc, std::shared_ptr > bias, const NEWIMAGE::volume& pred, const NEWIMAGE::volume& dima, const NEWIMAGE::volume& sims, const NEWIMAGE::volume& pmask, const NEWMAT::Matrix& XtX, const NEWMAT::ColumnVector& Xty, const NEWMAT::ColumnVector& update); /// Updates the parameters pertaining to modelling of time-varying EC static std::vector long_ec_update(// Input const std::vector >& pred, // Predictions in model space const NEWIMAGE::volume& pmask, // "Data valid" mask in model space float fwhm, // FWHM for Gaussian smoothing unsigned int nthr, // Number of threads to run it on bool very_verbose, // Detailed output to screen? // These input parameters are for debugging only unsigned int iter, // Iteration unsigned int level, // Determines how much gets written const std::vector& debug_index, // Indicies of scans to write debug info for // Input/output EDDY::ECScanManager& sm); // Scans we want to register to predictions /// Performs the levenberg updates of the long EC parameters. static void long_ec_levenberg_updates(// Input unsigned int first_vol, // Index of first vol to work on unsigned int last_vol, // One past index of last vol to work on const std::vector >& pred, // Predictions in model space const NEWIMAGE::volume& pmask, // "Data valid" mask in model space float fwhm, const std::vector& XtX, // Hessian for esimation of updates const std::vector& Xty, // Gradient for estimation of updates const std::vector& mss, // mss for current long EC parameters bool very_verbose, // Should we write info to screen? // Input used for debugging output only unsigned int iter, // Iteration unsigned int level, // Determines how much gets written const std::vector& dbg_indx, // Indicies of scans to write debug info for // Input/output EDDY::ECScanManager& sm, // Scan manager // Output std::vector& updates); // Applied updates /// Called by the function above to calculate matrices XtX and Xty. Allows for thread library parallelisation. static void long_ec_calculate_matrices(// Input unsigned int first_vol, // Index of first vol to work on unsigned int last_vol, // One past index of last vol to work on const std::vector >& pred, // Predictions in model space const NEWIMAGE::volume& pmask, // "Data valid" mask in model space float fwhm, // FWHM for Gaussian smoothing const EDDY::ECScanManager& sm, // Scan manager // These input parameters are for debugging only unsigned int iter, // Iteration unsigned int level, // Determines how much gets written const std::vector& debug_indx, // Indicies of scans to write debug info for // Output std::vector& mss, // Vector of sum-of-mean-square of difference std::vector& XtX, // Vector of XtX matrices std::vector& Xty); // Vector of Xty vectors static void write_debug_info_long_ec_calculate_matrices(const EDDY::ECScan& scan, unsigned int scindx, const std::vector& debug_indx, unsigned int iter, unsigned int level, const NEWIMAGE::volume4D& derivs, const NEWIMAGE::volume& mask, const NEWIMAGE::volume& mios, const NEWIMAGE::volume& pios, const NEWIMAGE::volume& pred, const NEWIMAGE::volume& dima, const NEWIMAGE::volume& sims, const NEWIMAGE::volume& pmask, const arma::mat& XtX, const arma::colvec& Xty); /// Used by EddyUtils::long_ec_update when checking if a joint update step was successful static double calculate_long_ec_total_new_mss(const std::vector >& pred, // Predictions in model space const NEWIMAGE::volume& pmask, // "Data valid" mask in model space float fwhm, // FWHM for Gaussian smoothing const EDDY::ECScanManager& sm, // Scan manager unsigned int nthr, // No. of threads // These input parameters are for debugging only unsigned int iter, // Iteration unsigned int level, // Determines how much gets written const std::vector& debug_indx); // Indicies of scans to write debug info for /// Helper function to parallelise calls to EddyUtils::calculate_long_ec_new_mss static void long_ec_total_mss_helper(// Input unsigned int first_vol, // Index of first vol to work on unsigned int last_vol, // One past index of last vol to work on const std::vector >& pred, // Predictions in model space const NEWIMAGE::volume& pmask, // "Data valid" mask in model space float fwhm, // FWHM for Gaussian smoothing const EDDY::ECScanManager& sm, // Scan manager // These input parameters are for debugging only unsigned int iter, // Iteration unsigned int level, // Determines how much gets written const std::vector& debug_indx, // Indicies of scans to write debug info for // Output std::vector mss); // Vector of mean-sum-of-square-differences static double calculate_new_mss(const NEWIMAGE::volume& pred, // Predictions in model space std::shared_ptr > susc, // Susceptibility-induced field const NEWIMAGE::volume& pmask, // "Data valid" mask in model space const EDDY::ECScan& scan, // Scan float fwhm, // FWHM for Gaussian smoothing EDDY::ParametersType whichp, // What parameters have been updated // These input parameters are for debugging only unsigned int scindx, // Scan index unsigned int iter, // Iteration unsigned int level, // Determines how much gets written const std::vector& debug_indx); // Indicies of scans to write debug info for /// Called by the EddyUtils::calculate_new_mss function above. It writes out debug info. static void write_debug_info_updated_mss(const EDDY::ECScan& scan, unsigned int scindx, EDDY::ParametersType whichp, const std::vector& debug_indx, unsigned int iter, unsigned int level, const NEWIMAGE::volume& pred, const NEWIMAGE::volume& mask, const NEWIMAGE::volume& dima, const NEWIMAGE::volume& pios, std::shared_ptr > susc); static EDDY::ImageCoordinates transform_coordinates_from_model_to_scan_space(// Input const NEWIMAGE::volume& pred, const EDDY::ECScan& scan, std::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, std::vector slices, // 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); static bool get_groups(// Input const std::vector& dpv, // Output std::vector >& grps, std::vector& grpi, std::vector& grpb); 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, std::shared_ptr > susc, bool jacmod, // Output NEWIMAGE::volume& omask, NEWIMAGE::volume *jac, NEWIMAGE::volume4D *grad); /* static NEWIMAGE::volume transform_model_to_scan_space(// Input const NEWIMAGE::volume& pred, const EDDY::ECScan& scan, std::shared_ptr > susc, bool jacmod, // Output NEWIMAGE::volume& omask, NEWIMAGE::volume *jac, NEWIMAGE::volume4D *grad, // Tmp unsigned int scindx, unsigned int iter, unsigned int level); */ /// Returns a string representation of the EC-model static std::string ModelString(ECModelType mt); /// Functions for setting/getting b_range static void SetbRange(double val) { b_range = val; } static double GetbRange() { return(b_range); } // Some functions for comparing diffusion parameters /// Returns true if the difference in b-value is less than EddyUtils::b_range static bool AreInSameShell(const DiffPara& dp1, const DiffPara& dp2) EddyTry { return(fabs(dp1.bVal()-dp2.bVal()) < EddyUtils::b_range); } EddyCatch static bool IsDiffusionWeighted(const DiffPara& dp) EddyTry { return(dp.bVal() > EddyUtils::b_range); } EddyCatch static bool Isb0(const DiffPara& dp) EddyTry { return(!IsDiffusionWeighted(dp)); } EddyCatch /// Returns true if the inner product of the b-vectors is greater than 0.999 static bool HaveSameDirection(const DiffPara& dp1, const DiffPara& dp2) EddyTry { return(NEWMAT::DotProduct(dp1.bVec(),dp2.bVec())>0.999); } EddyCatch /// Returns true if a vector of DiffPara objects implies a "shelled" (i.e. non-DSI) design static bool IsShelled(const std::vector& dpv); /// Returns true if a vector of DiffPara objects implies a multi-shell design static bool IsMultiShell(const std::vector& dpv); /// Returns a vector of group indices (one for each element in dpv) and b-values for the different groups static bool GetGroups(// Input const std::vector& dpv, // Output std::vector& grpi, std::vector& grpb); /// Returns n vectors of indicies into dpv, where n is the number of groups. Also b-values for each group. static bool GetGroups(// Input const std::vector& dpv, // Output std::vector >& grps, std::vector& grpb); /// Returns group info in both the formats of the routines above. static bool GetGroups(// Input const std::vector& dpv, // Output std::vector >& grps, std::vector& grpi, std::vector& grpb); // Random functions to set extrapolation and interpolation // template static void SetTrilinearInterp(V& vol) EddyTry { if (vol.getinterpolationmethod() != NEWIMAGE::trilinear) vol.setinterpolationmethod(NEWIMAGE::trilinear); if (vol.getextrapolationmethod() != NEWIMAGE::mirror) vol.setextrapolationmethod(NEWIMAGE::mirror); } EddyCatch template static void SetSplineInterp(V& vol) EddyTry { 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); } EddyCatch // 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); // Get vector of forward movement matrices, one per slice. static std::vector GetSliceWiseForwardMovementMatrices(const EDDY::ECScan& scan); // Get vector of inverse movement matrices, one per slice. static std::vector GetSliceWiseInverseMovementMatrices(const EDDY::ECScan& scan); // 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); // 3D Smooth 3D/4D volume within mask static NEWIMAGE::volume Smooth(const NEWIMAGE::volume& ima, // Image to smooth float fwhm, // FWHM of Gaussian kernel const NEWIMAGE::volume& mask); // Mask within which to smooth /// Negate the values in an std::vector of arma::colvec /// Used to set back parameters when they fail to improve the SSD. static std::vector Negate(const std::vector& vec); /// Checks for reasonable values in update. Based on empirical assessments. static bool UpdateMakesSense(const EDDY::ECScan& scan, const NEWMAT::ColumnVector& update); // Make image with normal distributed noise static NEWIMAGE::volume MakeNoiseIma(const NEWIMAGE::volume& ima, // Template image float mu, // Mean of noise float stdev); // Stdev of noise // Calculates slice-wise statistics from the difference between observation and predicton in observation space static DiffStats GetSliceWiseStats(// Input const EddyCommandLineOptions& clo, ScanType st, const NEWIMAGE::volume& pred, // Prediction in model space std::shared_ptr > susc, // Susceptibility induced off-resonance field const NEWIMAGE::volume& pmask, // "Data valid" mask in model space const ECScanManager& sm, unsigned int s, // Scan index unsigned int iter, // Iteration unsigned int dl); // Debug level // Returns a vector nvols such that thread i is responsible for the range of volumes nvols[i]--nvols[i+1]-1 static std::vector ScansPerThread(unsigned int nscans, // No. of scans/volumes unsigned int nthreads); // No. of threads // Performs an update of the movement parameters for one scan static double MovParamUpdate(// Input const NEWIMAGE::volume& pred, // Prediction in model space std::shared_ptr > susc, // Susceptibility off-resonance field std::shared_ptr > bias, // Recieve bias field const NEWIMAGE::volume& pmask, // "Data valid" mask in model space float fwhm, // FWHM for Gaussian smoothing bool very_verbose, // Input/output EDDY::ECScan& scan) EddyTry { return(param_update(pred,susc,bias,pmask,ParametersType::Movement,fwhm,very_verbose,0,0,0,scan)); } EddyCatch // Performs an update of the EC parameters for one scan static double ECParamUpdate(// Input const NEWIMAGE::volume& pred, // Prediction in model space std::shared_ptr > susc, // Susceptibility off-resonance field std::shared_ptr > bias, // Recieve bias field const NEWIMAGE::volume& pmask, // "Data valid" mask in model space float fwhm, // FWHM for Gaussian smoothing bool very_verbose, // Input/output EDDY::ECScan& scan) EddyTry { return(param_update(pred,susc,bias,pmask,ParametersType::EC,fwhm,very_verbose,0,0,0,scan)); } EddyCatch // Performs an update of the EC parameters for one scan // Does currently not use the bias parameter static double MovAndECParamUpdate(// Input const NEWIMAGE::volume& pred, // Prediction in model space std::shared_ptr > susc, // Susceptibility off-resonance field std::shared_ptr > bias, // Recieve bias field const NEWIMAGE::volume& pmask, // "Data valid" mask in model space float fwhm, // FWHM for Gaussian smoothing bool very_verbose, // Input/output EDDY::ECScan& scan) EddyTry { return(param_update(pred,susc,bias,pmask,ParametersType::All,fwhm,very_verbose,0,0,0,scan)); } EddyCatch // Performs an update of the EC parameters and writes debug into for one scan // Does currently not use the bias parameter static double MovAndECParamUpdate(// Input const NEWIMAGE::volume& pred, // Prediction in model space std::shared_ptr > susc, // Susceptibility off-resonance field std::shared_ptr > bias, // Recieve bias field const NEWIMAGE::volume& pmask, // "Data valid" mask in model space float fwhm, // FWHM for Gaussian smoothing bool very_verbose, // Parameters for debugging unsigned int scindx, // Scan index unsigned int iter, // Iteration unsigned int level, // Determines how much gets written // Input/output EDDY::ECScan& scan) EddyTry { return(param_update(pred,susc,bias,pmask,ParametersType::All,fwhm,very_verbose,scindx,iter,level,scan)); } EddyCatch // Performs an update of the long (time varying) EC parameters for all scans // Does currently not use the bias parameter static std::vector LongECParamUpdate(// Input const std::vector >& pred, // Predictions in model space const NEWIMAGE::volume& pmask, // "Data valid" mask in model space float fwhm, // FWHM for Gaussian smoothing unsigned int nthr, // Number of threads to run it on bool very_verbose, // Detailed output to screen? // Input/output EDDY::ECScanManager& sm) EddyTry { // Scans we want to register to predictions return(long_ec_update(pred,pmask,fwhm,nthr,very_verbose,0,0,std::vector(0),sm)); } EddyCatch // Performs an update of the long (time varying) EC parameters for all scans, and writes debug into // Does currently not use the bias parameter static std::vector LongECParamUpdate(// Input const std::vector >& pred, // Predictions in model space const NEWIMAGE::volume& pmask, // "Data valid" mask in model space float fwhm, // FWHM for Gaussian smoothing unsigned int nthr, // Number of threads to run it on bool very_verbose, // Detailed output to screen? // These input parameters are for debugging only unsigned int iter, // Iteration unsigned int level, // Determines how much gets written const std::vector& debug_index, // Indicies of scans to write debug info for // Input/output EDDY::ECScanManager& sm) EddyTry { // Scans we want to register to predictions return(long_ec_update(pred,pmask,fwhm,nthr,very_verbose,iter,level,debug_index,sm)); } EddyCatch // Transforms an image from model/prediction space to observation space static NEWIMAGE::volume TransformModelToScanSpace(// Input const NEWIMAGE::volume& pred, const EDDY::ECScan& scan, std::shared_ptr > susc) EddyTry { 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)); } EddyCatch static NEWIMAGE::volume TransformScanToModelSpace(// Input const EDDY::ECScan& scan, std::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, std::shared_ptr > susc, // Output NEWIMAGE::volume& omask); /* static NEWIMAGE::volume DirectTransformModelToScanSpace(// Input const NEWIMAGE::volume& ima, const EDDY::ECScan& scan, const EDDY::MultiBandGroups& mbg, std::shared_ptr > susc, // Output NEWIMAGE::volume& omask); */ static NEWIMAGE::volume4D DerivativesForModelToScanSpaceTransform(const EDDY::ECScan& scan, const NEWIMAGE::volume& mima, std::shared_ptr > susc) { return(EddyUtils::get_partial_derivatives_in_scan_space(mima,scan,susc,EDDY::ParametersType::All)); } static NEWIMAGE::volume4D DirectDerivativesForModelToScanSpaceTransform(const EDDY::ECScan& scan, const NEWIMAGE::volume& mima, std::shared_ptr > susc) EddyTry { return(EddyUtils::get_direct_partial_derivatives_in_scan_space(mima,scan,susc,EDDY::ParametersType::All)); } EddyCatch }; //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ // // 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) EddyTry { NEWIMAGE::volume4D mmd=voxdisp; mmd[0] *= mmd.xdim(); mmd[1] *= mmd.ydim(); mmd[2] *= mmd.zdim(); return(mmd); } EddyCatch // Inverts 3D displacement field, ASSUMING it is really 1D (only non-zero displacements in one direction). static NEWIMAGE::volume4D Invert3DDisplacementField(// 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 Invert1DDisplacementField(// 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: }; /****************************************************************//** * * \brief This class estimates amount of s2v movement * * This class estimates amount of s2v movement * ********************************************************************/ class s2vQuant { public: s2vQuant(const ECScanManager& sm) EddyTry : _sm(sm), _trth(0.3), _rotth(0.3) { common_construction(); } EddyCatch s2vQuant(const ECScanManager& sm, double trth, double rotth) EddyTry : _sm(sm), _trth(trth), _rotth(rotth) { common_construction(); } EddyCatch /// Returns a vector of indicies of "still" volumes std::vector FindStillVolumes(ScanType st, const std::vector& mbsp) const; private: /// Performs common construction tasks void common_construction(); const ECScanManager& _sm; ///< Local copy of the scan manager NEWMAT::Matrix _tr; ///< Matrix with across-volume std of translations NEWMAT::Matrix _rot; ///< Matrix with across-volume std of rotations double _trth; ///< Threshold for mean translation std for being considered "still" double _rotth; ///< Threshold for mean rotation std for being considered "still" }; } // End namespace EDDY #endif // End #ifndef EddyUtils_h