EddyUtils.h

// 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 <cstdlib>
#include <string>
#include <vector>
#include <cmath>
#include <memory>
#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<float> get_partial_derivatives_in_scan_space(// Input
									 const NEWIMAGE::volume<float>&                    pred,      // Prediction in model space
									 const EDDY::ECScan&                               scan,      // Scan space
									 std::shared_ptr<const NEWIMAGE::volume<float> >   susc,      // Susceptibility off-resonance field
                                                                         EDDY::ParametersType                              whichp);

  static NEWIMAGE::volume4D<float> get_direct_partial_derivatives_in_scan_space(// Input
										const NEWIMAGE::volume<float>&                    pred,     // Prediction in model space
										const EDDY::ECScan&                               scan,     // Scan space
										std::shared_ptr<const NEWIMAGE::volume<float> >   susc,     // Susceptibility off-resonance field
										EDDY::ParametersType                              whichp);

  static NEWIMAGE::volume4D<float> get_long_ec_partial_derivatives_in_scan_space(// Input
										 const NEWIMAGE::volume<float>&                  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<const NEWIMAGE::volume<float> > susc); // Susceptibility off-resonance field

  static double param_update(// Input
			     const NEWIMAGE::volume<float>&                      pred,      // Prediction in model space
			     std::shared_ptr<const NEWIMAGE::volume<float> >     susc,      // Susceptibility off-resonance field
			     std::shared_ptr<const NEWIMAGE::volume<float> >     bias,      // Recieve bias field
			     const NEWIMAGE::volume<float>&                      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<float>&                   derivs,
						const NEWIMAGE::volume<float>&                     mask,
						const NEWIMAGE::volume<float>&                     mios,
						const NEWIMAGE::volume<float>&                     pios,
						const NEWIMAGE::volume<float>&                     jac,
						std::shared_ptr<const NEWIMAGE::volume<float> >    susc,
						std::shared_ptr<const NEWIMAGE::volume<float> >    bias,
						const NEWIMAGE::volume<float>&                     pred,
						const NEWIMAGE::volume<float>&                     dima,
						const NEWIMAGE::volume<float>&                     sims,
						const NEWIMAGE::volume<float>&                     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<double> long_ec_update(// Input
					    const std::vector<NEWIMAGE::volume<float> >&      pred,         // Predictions in model space
					    const NEWIMAGE::volume<float>&                    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<unsigned int>&                  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<NEWIMAGE::volume<float> >& pred,           // Predictions in model space
					const NEWIMAGE::volume<float>&               pmask,          // "Data valid" mask in model space
					float                                        fwhm,
					const std::vector<arma::mat>&                XtX,            // Hessian for esimation of updates
					const std::vector<arma::colvec>&             Xty,            // Gradient for estimation of updates
					const std::vector<double>&                   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<unsigned int>&             dbg_indx,       // Indicies of scans to write debug info for
					// Input/output
					EDDY::ECScanManager&                         sm,             // Scan manager
					// Output
					std::vector<arma::colvec>&                   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<NEWIMAGE::volume<float> >&    pred,              // Predictions in model space
					 const NEWIMAGE::volume<float>&                  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<unsigned int>&                debug_indx,        // Indicies of scans to write debug info for
					 // Output
					 std::vector<double>&                            mss,               // Vector of sum-of-mean-square of difference
					 std::vector<arma::mat>&                         XtX,               // Vector of XtX matrices
					 std::vector<arma::colvec>&                      Xty);              // Vector of Xty vectors

  static void write_debug_info_long_ec_calculate_matrices(const EDDY::ECScan&                                scan,
							  unsigned int                                       scindx,
							  const std::vector<unsigned int>&                   debug_indx,
							  unsigned int                                       iter,
							  unsigned int                                       level,
							  const NEWIMAGE::volume4D<float>&                   derivs,
							  const NEWIMAGE::volume<float>&                     mask,
							  const NEWIMAGE::volume<float>&                     mios,
							  const NEWIMAGE::volume<float>&                     pios,
							  const NEWIMAGE::volume<float>&                     pred,
							  const NEWIMAGE::volume<float>&                     dima,
							  const NEWIMAGE::volume<float>&                     sims,
							  const NEWIMAGE::volume<float>&                     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<NEWIMAGE::volume<float> >& pred,               // Predictions in model space
						const NEWIMAGE::volume<float>&               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<unsigned int>&             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<NEWIMAGE::volume<float> >& pred,           // Predictions in model space
				       const NEWIMAGE::volume<float>&               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<unsigned int>&             debug_indx,     // Indicies of scans to write debug info for
				       // Output
				       std::vector<double>                          mss);           // Vector of mean-sum-of-square-differences 
  
  static double calculate_new_mss(const NEWIMAGE::volume<float>&                  pred,               // Predictions in model space
				  std::shared_ptr<const NEWIMAGE::volume<float> > susc,               // Susceptibility-induced field
				  const NEWIMAGE::volume<float>&                  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<unsigned int>&                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<unsigned int>&                   debug_indx,
					   unsigned int                                       iter,
					   unsigned int                                       level,
					   const NEWIMAGE::volume<float>&                     pred,
					   const NEWIMAGE::volume<float>&                     mask,
					   const NEWIMAGE::volume<float>&                     dima,
					   const NEWIMAGE::volume<float>&                     pios,
					   std::shared_ptr<const NEWIMAGE::volume<float> >    susc);

  static EDDY::ImageCoordinates transform_coordinates_from_model_to_scan_space(// Input
									       const NEWIMAGE::volume<float>&                    pred,
									       const EDDY::ECScan&                               scan,
									       std::shared_ptr<const NEWIMAGE::volume<float> >   susc,
									       // Output
									       NEWIMAGE::volume<float>                           *omask,
									       NEWIMAGE::volume<float>                           *jac);
  // Returns coordinates into f transformed with
  // displacement field d and affine matrix M
  static void transform_coordinates(// Input
				    const NEWIMAGE::volume<float>&    f,
				    const NEWIMAGE::volume4D<float>&  d,
				    const NEWMAT::Matrix&             M,
				    std::vector<unsigned int>         slices,
				    // Output
				    ImageCoordinates&                 c,
				    NEWIMAGE::volume<float>           *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<float>& vols,
				 const NEWIMAGE::volume<float>&   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<float>& Xvols,
				       const NEWIMAGE::volume<float>&   Yvol,
				       const NEWIMAGE::volume<float>&   mask);
  static bool get_groups(// Input
			 const std::vector<DiffPara>&             dpv,
			 // Output
			 std::vector<std::vector<unsigned int> >& grps,
			 std::vector<unsigned int>&               grpi,
			 std::vector<double>&                     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<float> transform_model_to_scan_space(// Input
							       const NEWIMAGE::volume<float>&                    pred,
							       const EDDY::ECScan&                               scan,
							       std::shared_ptr<const NEWIMAGE::volume<float> >   susc,
							       bool                                              jacmod,
							       // Output
							       NEWIMAGE::volume<float>&                          omask,
							       NEWIMAGE::volume<float>                           *jac,
							       NEWIMAGE::volume4D<float>                         *grad);

  /*
  static NEWIMAGE::volume<float> transform_model_to_scan_space(// Input
							       const NEWIMAGE::volume<float>&                    pred,
							       const EDDY::ECScan&                               scan,
							       std::shared_ptr<const NEWIMAGE::volume<float> >   susc,
							       bool                                              jacmod,
							       // Output
							       NEWIMAGE::volume<float>&                          omask,
							       NEWIMAGE::volume<float>                           *jac,
							       NEWIMAGE::volume4D<float>                         *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<DiffPara>& dpv);
  /// Returns true if a vector of DiffPara objects implies a multi-shell design
  static bool IsMultiShell(const std::vector<DiffPara>& 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<DiffPara>&             dpv,
			// Output
			std::vector<unsigned int>&               grpi,
			std::vector<double>&                     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<DiffPara>&             dpv,
			// Output
			std::vector<std::vector<unsigned int> >& grps,
			std::vector<double>&                     grpb);
  /// Returns group info in both the formats of the routines above.
  static bool GetGroups(// Input
			const std::vector<DiffPara>&             dpv,
			// Output
			std::vector<std::vector<unsigned int> >& grps,
			std::vector<unsigned int>&               grpi,
			std::vector<double>&                     grpb);
  // Random functions to set extrapolation and interpolation //
  template <class V>
  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 <class V>
  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<unsigned int> GetIndiciesOfDWIs(const std::vector<DiffPara>& dpars);

  // Get vector of forward movement matrices, one per slice.
  static std::vector<NEWMAT::Matrix> GetSliceWiseForwardMovementMatrices(const EDDY::ECScan& scan);

  // Get vector of inverse movement matrices, one per slice.
  static std::vector<NEWMAT::Matrix> GetSliceWiseInverseMovementMatrices(const EDDY::ECScan& scan);

  // Removes bvecs associated with zero b-values.
  static std::vector<DiffPara> GetDWIDiffParas(const std::vector<DiffPara>&   dpars);

  // Reads all diffusion weighted images from 4D volume
  static int read_DWI_volume4D(NEWIMAGE::volume4D<float>&     dwivols,
			       const std::string&             fname,
			       const std::vector<DiffPara>&   dpars);

  // Converts char mask (from the general_transform functions) to a float mask
  static NEWIMAGE::volume<float> ConvertMaskToFloat(const NEWIMAGE::volume<char>& charmask);

  // 3D Smooth 3D/4D volume within mask
  static NEWIMAGE::volume<float> Smooth(const NEWIMAGE::volume<float>& ima,   // Image to smooth
					float                          fwhm,  // FWHM of Gaussian kernel
					const NEWIMAGE::volume<float>& 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<arma::colvec> Negate(const std::vector<arma::colvec>& 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<float> MakeNoiseIma(const NEWIMAGE::volume<float>&   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<float>&                  pred,     // Prediction in model space
				     std::shared_ptr<const NEWIMAGE::volume<float> > susc,     // Susceptibility induced off-resonance field
				     const NEWIMAGE::volume<float>&                  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<unsigned int> 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<float>&                    pred,      // Prediction in model space
			       std::shared_ptr<const NEWIMAGE::volume<float> >   susc,      // Susceptibility off-resonance field
			       std::shared_ptr<const NEWIMAGE::volume<float> >   bias,      // Recieve bias field
			       const NEWIMAGE::volume<float>&                    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<float>&                    pred,      // Prediction in model space
			      std::shared_ptr<const NEWIMAGE::volume<float> >   susc,      // Susceptibility off-resonance field
			      std::shared_ptr<const NEWIMAGE::volume<float> >   bias,      // Recieve bias field
			      const NEWIMAGE::volume<float>&                    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<float>&                    pred,      // Prediction in model space
				    std::shared_ptr<const NEWIMAGE::volume<float> >   susc,      // Susceptibility off-resonance field
				    std::shared_ptr<const NEWIMAGE::volume<float> >   bias,      // Recieve bias field
				    const NEWIMAGE::volume<float>&                    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<float>&                    pred,      // Prediction in model space
				    std::shared_ptr<const NEWIMAGE::volume<float> >   susc,      // Susceptibility off-resonance field
				    std::shared_ptr<const NEWIMAGE::volume<float> >   bias,      // Recieve bias field
				    const NEWIMAGE::volume<float>&                    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<double> LongECParamUpdate(// Input
					       const std::vector<NEWIMAGE::volume<float> >&      pred,         // Predictions in model space
					       const NEWIMAGE::volume<float>&                    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<unsigned int>(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<double> LongECParamUpdate(// Input
					       const std::vector<NEWIMAGE::volume<float> >&      pred,         // Predictions in model space
					       const NEWIMAGE::volume<float>&                    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<unsigned int>&                  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<float> TransformModelToScanSpace(// Input
							   const NEWIMAGE::volume<float>&                    pred,
							   const EDDY::ECScan&                               scan,
							   std::shared_ptr<const NEWIMAGE::volume<float> >   susc) EddyTry {
    NEWIMAGE::volume<float> mask(pred.xsize(),pred.ysize(),pred.zsize());
    NEWIMAGE::volume<float> jac(pred.xsize(),pred.ysize(),pred.zsize());
    return(transform_model_to_scan_space(pred,scan,susc,true,mask,&jac,NULL));
  } EddyCatch
  static NEWIMAGE::volume<float> TransformScanToModelSpace(// Input
							   const EDDY::ECScan&                             scan,
							   std::shared_ptr<const NEWIMAGE::volume<float> > susc,
							   // Output
							   NEWIMAGE::volume<float>&                        omask);

  // The next two are alternate transformation routines that
  // performs the transforms in several resampling steps.
  // They are intended for debugging.
  static NEWIMAGE::volume<float> DirectTransformScanToModelSpace(// Input
								 const EDDY::ECScan&                             scan,
								 std::shared_ptr<const NEWIMAGE::volume<float> > susc,
								 // Output
								 NEWIMAGE::volume<float>&                        omask);
  /*
  static NEWIMAGE::volume<float> DirectTransformModelToScanSpace(// Input
								 const NEWIMAGE::volume<float>&                    ima,
								 const EDDY::ECScan&                               scan,
								 const EDDY::MultiBandGroups&                      mbg,
								 std::shared_ptr<const NEWIMAGE::volume<float> >   susc,
								 // Output
								 NEWIMAGE::volume<float>&                          omask);
  */
  static NEWIMAGE::volume4D<float> DerivativesForModelToScanSpaceTransform(const EDDY::ECScan&                             scan,
									   const NEWIMAGE::volume<float>&                  mima,
									   std::shared_ptr<const NEWIMAGE::volume<float> > susc)
  {
    return(EddyUtils::get_partial_derivatives_in_scan_space(mima,scan,susc,EDDY::ParametersType::All));
  }

  static NEWIMAGE::volume4D<float> DirectDerivativesForModelToScanSpaceTransform(const EDDY::ECScan&                               scan,
										 const NEWIMAGE::volume<float>&                    mima,
										 std::shared_ptr<const NEWIMAGE::volume<float> >   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<float> RigidBodyTransformHzField(const NEWIMAGE::volume<float>& hzfield);

  // Some conversion routines off-resonance->displacements
  static NEWIMAGE::volume4D<float> Hz2VoxelDisplacements(const NEWIMAGE::volume<float>& hzfield,
                                                         const AcqPara&                 acqp);
  static NEWIMAGE::volume4D<float> Hz2MMDisplacements(const NEWIMAGE::volume<float>& hzfield,
                                                      const AcqPara&                 acqp);
  static NEWIMAGE::volume4D<float> Voxel2MMDisplacements(const NEWIMAGE::volume4D<float>& voxdisp) EddyTry {
    NEWIMAGE::volume4D<float> 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<float> Invert3DDisplacementField(// Input
							     const NEWIMAGE::volume4D<float>& dfield,
							     const AcqPara&                   acqp,
							     const NEWIMAGE::volume<float>& inmask,
							     // Output
							     NEWIMAGE::volume<float>&       omask);

  // Inverts 1D displacement field. Input must be scaled to voxels (i.e. not mm).
  static NEWIMAGE::volume<float> Invert1DDisplacementField(// Input
							   const NEWIMAGE::volume<float>& dfield,
							   const AcqPara&                 acqp,
							   const NEWIMAGE::volume<float>& inmask,
							   // Output
							   NEWIMAGE::volume<float>&       omask);

  // Calculates Jacobian of a displacement field
  static NEWIMAGE::volume<float> GetJacobian(const NEWIMAGE::volume4D<float>& dfield,
                                             const AcqPara&                   acqp);

  // Calculates Jacobian of a 1D displacement field
  static NEWIMAGE::volume<float> GetJacobianFrom1DField(const NEWIMAGE::volume<float>& 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<unsigned int> FindStillVolumes(ScanType st, const std::vector<unsigned int>& 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