/*! \file DiffusionGP.h \brief Contains declaration of class for making Gaussian process based predictions about DWI data. \author Jesper Andersson \version 1.0b, Sep., 2012. */ // Declarations of class to make Gaussian-Process // based predictions about diffusion data. // // DiffusionGP.h // // Jesper Andersson, FMRIB Image Analysis Group // // Copyright (C) 2011 University of Oxford // #ifndef DiffusionGP_h #define DiffusionGP_h #include #include #include #include #include #include "armawrap/newmat.h" #include "newimage/newimageall.h" #include "miscmaths/miscmaths.h" #include "EddyHelperClasses.h" #include "DWIPredictionMaker.h" #include "KMatrix.h" #include "HyParEstimator.h" namespace EDDY { /****************************************************************//** * * \brief Class used to make Gaussian process based predictions * about diffusion data. * * Will make predictions about observed (smoothing) and unobserved * (interpolation) scans using Gaussian Processes. It works by first * creating an object, setting the # of scans one wants it to contain * and then populate it with scans. Once that is done it is ready * to start making predictions. If you think of the signal (in a * given voxel) as points on a surface the Gaussian process will make * predictions from the assumption that that surface should be smooth. * ********************************************************************/ class DiffusionGP : public DWIPredictionMaker { public: /// Default constructor DiffusionGP(const std::shared_ptr& Kmat, const std::shared_ptr& hpe) EddyTry : _Kmat(Kmat->Clone()), _hpe(hpe->Clone()), _pop(true), _mc(false) {} EddyCatch /// Constructor that takes filenames from which to load data DiffusionGP(const std::shared_ptr& Kmat, const std::shared_ptr& hpe, const std::string& scans_fname, const std::string& var_mask_fname, const std::vector& dpars, float fwhm=0.0, bool verbose=false); /// Returns prediction for point given by indx. virtual NEWIMAGE::volume Predict(unsigned int indx, bool exclude=false) const; /// Returns prediction for point given by indx. virtual NEWIMAGE::volume Predict(unsigned int indx, bool exclude=false); /// Returns prediction for point given by indx. This is used only as a means of directly comparing CPU and GPU outputs. virtual NEWIMAGE::volume PredictCPU(unsigned int indx, bool exclude=false); /// Returns predictions for points given by indicies virtual std::vector > Predict(const std::vector& indicies, bool exclude=false); /// Returns input data for point given by indx virtual NEWIMAGE::volume InputData(unsigned int indx) const; /// Returns input data for points given by indicies virtual std::vector > InputData(const std::vector& indicies) const; /// Returns variance of prediction for point given by indx. virtual double PredictionVariance(unsigned int indx, bool exclude=false); /// Returns measurement-error variance for point given by indx. virtual double ErrorVariance(unsigned int indx) const; /// Returns true if all data has been loaded virtual bool IsPopulated() const { return(_pop); } /// Indicates if it is ready to make predictions. virtual bool IsValid() const EddyTry { return(IsPopulated() && _mc && _Kmat->IsValid()); } EddyCatch /// Specify the # of points we plan to put into the predictor. virtual void SetNoOfScans(unsigned int n); /// Set a point given by indx. This function is thread safe as long as different threads set different points. virtual void SetScan(const NEWIMAGE::volume& scan, // May be thread safe if used "sensibly" const DiffPara& dp, unsigned int indx, unsigned int sess) { SetScan(scan,dp,indx); } /// Set a point given by indx. This function is thread safe as long as different threads set different points. virtual void SetScan(const NEWIMAGE::volume& scan, // May be thread safe if used "sensibly" const DiffPara& dp, unsigned int indx); /// Returns the number of hyperparameters for the model virtual unsigned int NoOfHyperPar() const EddyTry { return(_Kmat->NoOfHyperPar()); } EddyCatch /// Returns the hyperparameters for the model virtual std::vector GetHyperPar() const EddyTry { return(_Kmat->GetHyperPar()); } EddyCatch /// Evaluates the model so as to make the predictor ready to make predictions. virtual void EvaluateModel(const NEWIMAGE::volume& mask, float fwhm, bool verbose=false); /// Evaluates the model so as to make the predictor ready to make predictions. virtual void EvaluateModel(const NEWIMAGE::volume& mask, bool verbose=false) EddyTry { EvaluateModel(mask,0.0,verbose); } EddyCatch /// Writes internal content to disk for debug purposes virtual void WriteImageData(const std::string& fname) const; virtual void WriteMetaData(const std::string& fname) const EddyTry { if (!IsPopulated()) throw EddyException("DiffusionGP::WriteMetaData: Not yet fully populated"); _Kmat->Write(fname); } EddyCatch virtual void Write(const std::string& fname) const EddyTry { WriteImageData(fname); WriteMetaData(fname); } EddyCatch private: static const unsigned int nvoxhp = 500; std::vector > > _sptrs; // Pointers to the scans std::shared_ptr _Kmat; // K-matrix for Gaussian Process std::vector _dpars; // Diffusion parameters std::shared_ptr _hpe; // Hyperparameter estimator bool _pop; // True if all data is present bool _mc; // True if data has been mean-corrected std::vector > > _mptrs; // Pointers to mean images void mean_correct(const std::vector >& mi); bool is_populated() const; void predict_image_cpu(unsigned int indx, bool excl, const NEWMAT::RowVector& pv, NEWIMAGE::volume& ima) const; void predict_images_cpu(// Input const std::vector& indicies, bool exclude, const std::vector& pvecs, // Output std::vector >& pi) const; #ifdef COMPILE_GPU void predict_image_gpu(unsigned int indx, bool excl, const NEWMAT::RowVector& pv, NEWIMAGE::volume& ima) const; void predict_images_gpu(// Input const std::vector& indicies, bool exclude, const std::vector& pvecs, // Output std::vector >& pi) const; #endif unsigned int which_mean(unsigned int indx) const; bool get_y(// Input unsigned int i, unsigned int j, unsigned int k, unsigned int indx, bool exclude, // Output NEWMAT::ColumnVector& y) const; }; } // End namespace EDDY #endif // End #ifndef DiffusionGP_h