// Definitions of a class that does the parsing and
// sanity checking of commnad line options for the
// "eddy" application.
//
// EddyCommandLineOptions.cpp
//
// Jesper Andersson, FMRIB Image Analysis Group
//
// Copyright (C) 2011 University of Oxford
//
/*  Part of FSL - FMRIB's Software Library
    http://www.fmrib.ox.ac.uk/fsl
    fsl@fmrib.ox.ac.uk

    Developed at FMRIB (Oxford Centre for Functional Magnetic Resonance
    Imaging of the Brain), Department of Clinical Neurology, Oxford
    University, Oxford, UK


    LICENCE

    FMRIB Software Library, Release 6.0 (c) 2018, The University of
    Oxford (the "Software")

    The Software remains the property of the Oxford University Innovation
    ("the University").

    The Software is distributed "AS IS" under this Licence solely for
    non-commercial use in the hope that it will be useful, but in order
    that the University as a charitable foundation protects its assets for
    the benefit of its educational and research purposes, the University
    makes clear that no condition is made or to be implied, nor is any
    warranty given or to be implied, as to the accuracy of the Software,
    or that it will be suitable for any particular purpose or for use
    under any specific conditions. Furthermore, the University disclaims
    all responsibility for the use which is made of the Software. It
    further disclaims any liability for the outcomes arising from using
    the Software.

    The Licensee agrees to indemnify the University and hold the
    University harmless from and against any and all claims, damages and
    liabilities asserted by third parties (including claims for
    negligence) which arise directly or indirectly from the use of the
    Software or the sale of any products based on the Software.

    No part of the Software may be reproduced, modified, transmitted or
    transferred in any form or by any means, electronic or mechanical,
    without the express permission of the University. The permission of
    the University is not required if the said reproduction, modification,
    transmission or transference is done without financial return, the
    conditions of this Licence are imposed upon the receiver of the
    product, and all original and amended source code is included in any
    transmitted product. You may be held legally responsible for any
    copyright infringement that is caused or encouraged by your failure to
    abide by these terms and conditions.

    You are not permitted under this Licence to use this Software
    commercially. Use for which any financial return is received shall be
    defined as commercial use, and includes (1) integration of all or part
    of the source code or the Software into a product for sale or license
    by or on behalf of Licensee to third parties or (2) use of the
    Software or any derivative of it for research with the final aim of
    developing software products for sale or license to a third party or
    (3) use of the Software or any derivative of it for research with the
    final aim of developing non-software products for sale or license to a
    third party, or (4) use of the Software to provide any service to an
    external organisation for which payment is received. If you are
    interested in using the Software commercially, please contact Oxford
    University Innovation ("OUI"), the technology transfer company of the
    University, to negotiate a licence. Contact details are:
    fsl@innovation.ox.ac.uk quoting Reference Project 9564, FSL.*/
#include <cstdlib>
#include <string>
#include <vector>
#include <algorithm>
#include <thread>
#include "utils/options.h"
#include "topup/topup_file_io.h"
#include "EddyHelperClasses.h"
#include "EddyUtils.h"
#include "EddyCommandLineOptions.h"

using namespace std;
using namespace EDDY;

EddyCommandLineOptions::diff::diff() try :
  // Parameters related to specification of diffusion weighting
  _bvecs(string("--bvecs"),string(""),string("\tFile containing the b-vectors for all volumes in --imain"),true,Utilities::requires_argument),
  _bvals(string("--bvals"),string(""),string("\tFile containing the b-values for all volumes in --imain"),true,Utilities::requires_argument),
  _bdeltas(string("--bdeltas"),string(""),string("File containing the b-delta values for all volumes in --imain"),false,Utilities::requires_argument),
  _data_is_shelled(string("--data_is_shelled"),false,string("Assume, don't check, that data is shelled (default false)"),false, Utilities::no_argument),
  _brange(string("--b_range"),100,string("Range of b-values considered to be the same shell (default 100)"),false, Utilities::requires_argument),
  // Parameters related to least squares reconstruction of the data
  _resamp(string("--resamp"),string("jac"),string("Final resampling method (jac/lsr, default jac)"),false, Utilities::requires_argument),
  _lsr_lambda(string("--lsr_lambda"),0.01,string("Regularisation weight for LSR-resampling."),false,Utilities::requires_argument),
  // Parameters related to eddy current modelling
  _flm(string("--flm"),string("quadratic"),string("\tFirst level EC model (movement/linear/quadratic/cubic, default quadratic)"),false,Utilities::requires_argument),
  _slm_str(string("--slm"),string("none"),string("\tSecond level EC model (none/linear/quadratic, default none)"),false,Utilities::requires_argument),
  _b0_flm(string("--b0_flm"),string("movement"),string("First level EC model for b0 scans (movement/linear/quadratic, default movement)"),false,Utilities::requires_argument),
  _b0_slm_str(string("--b0_slm"),string("none"),string("Second level EC model for b0 scans (none/linear/quadratic, default none)"),false,Utilities::requires_argument),
  _rep_time(string("--rep_time"),5.0,string("Repetition time (seconds)"),false,Utilities::requires_argument),
  _long_ec_str(string("--lecm"),string("none"),string("Model for long (time-varying) EC (none/weights/tc/jointweights/jointtc default none)"),false,Utilities::requires_argument),
  _long_ec_dont_reest(string("--lecm_dont_reest"),false,string("Don't reestimate EC when estimating long EC (default false)"),false, Utilities::no_argument),
  _long_ec_sep_offs_move(string("--lecm_sep_offs_move"),false,string("Attempt to separate field offset from subject movement when estimating long EC (default false)"),false, Utilities::no_argument),
  // Parameters related to outlier detection/replacement
  _rep_noise(string("--rep_noise"),false,string("Add noise to replaced outliers (default false)"),false, Utilities::no_argument),
  _ol_ss(string("--ol_ss"),string("sw"),string("Summary statistics for OL, shell-wise (sw) or pooled (pooled). (default sw)"),false,Utilities::requires_argument),
  _ol_pos(string("--ol_pos"),false,string("Consider both positive and negative outliers if set (default false)"),false,Utilities::no_argument),
  _ol_sqr(string("--ol_sqr"),false,string("Consider outliers among sums-of-squared differences if set (default false)"),false,Utilities::no_argument),
  // Parameters related to the Gaussian process
  _covfunc(string("--covfunc"),string("spheri"),string("Covariance function for GP (spheri/expo/old, default spheri)"),false, Utilities::requires_argument),
  // Should bvecs be rotated during estimations (probably not)
  _rbvde(string("--rbvde"),false,string("\tRotate b-vecs during estimation (default false)"),false,Utilities::no_argument),
  // Parameters related to "tricky" movements for diffusion
  _dont_sep_offs_move(string("--dont_sep_offs_move"),false,string("Do NOT attempt to separate field offset from subject movement (default false)"),false, Utilities::no_argument),
  _offset_model(string("--offset_model"),string("linear_lag"),string("Second level model for field offset (default linear_lag)"),false, Utilities::requires_argument),
  _dont_peas(string("--dont_peas"),false,string("Do NOT perform a post-eddy alignment of shells (default false)"),false, Utilities::no_argument),
  _use_b0s_for_peas(string("--b0_peas"),false,string("Use interspersed b0s to perform post-eddy alignment of shells (default false)"),false, Utilities::no_argument),
  _session(string("--session"),string(""),string("File containing session indices for all volumes in --imain"),false,Utilities::requires_argument),
  // Parameters related to the writing of various "extra" output
  _fields(string("--fields"),false,string("Write EC fields as images (default false)"),false, Utilities::no_argument),
  _write_cnr_maps(string("--cnr_maps"),false,string("Write shell-wise cnr-maps (default false)"),false, Utilities::no_argument),
  _write_range_cnr_maps(string("--range_cnr_maps"),false,string("Write shell-wise range-cnr-maps (default false)"),false, Utilities::no_argument),
  _write_scatter_brain_predictions(string("--write_scatter_brain_predictions"),false,string("Write predictions obtained with a scattered data approach (in addition to corrected, default false)"),false, Utilities::no_argument),
  // Parameters related to things that can be useful for debugging/development
  _dwi_only(string("--dwi_only"),false,string("Only register the dwi images (default false)"),false, Utilities::no_argument),
  _b0_only(string("--b0_only"),false,string("Only register the b0 images (default false)"),false, Utilities::no_argument),
  _test_rot(string("--test_rot"),std::vector<float>(0),string("Do a large rotation to test b-vecs"),false,Utilities::requires_argument),
  _print_mi_values(string("--pmiv"),false,string("Write text file of MI values between shells (default false)"),false,Utilities::no_argument),
  _print_mi_planes(string("--pmip"),false,string("Write text file of (2D) MI values between shells (default false)"),false,Utilities::no_argument),
  // --mb_offs will soon be deprecated, so we should warn about it
  _mb_offs(string("--mb_offs"),0,string("Multi-band offset (-1 if bottom slice removed, 1 if top slice removed) (DEPRECATED!)"),false,Utilities::requires_argument),
  // The remaining are deprecated parameters
  _sep_offs_move(string("--sep_offs_move"),true,string("Attempt to separate field offset from subject movement (DEPRECATED!)"),false, Utilities::no_argument),
  _peas(string("--peas"),true,string("Perform a post-eddy alignment of shells (DEPRECATED!)"),false, Utilities::no_argument),
  _rms(string("--rms"),true,string("Write movement induced RMS (DEPRECATED!)"),false, Utilities::no_argument),
  _slorder(string("--slorder"),string(""),string("Name of text file defining slice/group order (DEPRECATED!)"),false,Utilities::requires_argument)
{} catch (const EddyInputError& e)
{
  cout << e.what() << endl;
  cout << "Terminating program" << endl;
  exit(EXIT_FAILURE);
} EddyCatch

EddyCommandLineOptions::fmri::fmri() try :
  _session(string("--session"),string(""),string("File containing session indices for all volumes in --imain"),true,Utilities::requires_argument),
  _rep_time(string("--rep_time"),3.0,string("Repetition time (seconds)"),true,Utilities::requires_argument),
  _echo_time(string("--echo_time"),40.0,string("Echo time (milliseconds)"),false,Utilities::requires_argument),
  _covfunc(string("--covfunc"),string("sqexp"),string("Covariance function for GP (sqexp, default sqexp)"),false, Utilities::requires_argument),
  _write_snr_maps(string("--snr_maps"),false,string("Write session-wise SNR maps (default false)"),false, Utilities::no_argument),
  _fast_map(string("--fast_map"),string(""),string("FAST pve map for CSF"),false, Utilities::requires_argument)
{}
catch (const EddyInputError& e)
{
  cout << e.what() << endl;
  cout << "Terminating program" << endl;
  exit(EXIT_FAILURE);
} EddyCatch

EddyCommandLineOptions::EddyCommandLineOptions(int  argc, char *argv[]) try :
  _title("eddy \nCopyright(c) 2015, University of Oxford (Jesper Andersson)"),
  _examples("eddy --monsoon"),
  // Mandatory parameters
  _imain(string("--imain"),string(""),string("\tFile containing all the images to estimate distortions for"),true,Utilities::requires_argument),
  _mask(string("--mask"),string(""),string("\tMask to indicate brain"),true,Utilities::requires_argument),
  _acqp(string("--acqp"),string(""),string("\tFile containing acquisition parameters"),true,Utilities::requires_argument),
  _index(string("--index"),string(""),string("\tFile containing indices for all volumes in --imain into --acqp and --topup"),true,Utilities::requires_argument),
  _out(string("--out"),string(""),string("\tBasename for output"),true,Utilities::requires_argument),
  // Parameters related to input susceptibility field
  _topup(string("--topup"),string(""),string("\tBase name for output files from topup"),false,Utilities::requires_argument),
  _field(string("--field"),string(""),string("\tName of file with susceptibility field (in Hz)"),false,Utilities::requires_argument),
  _field_mat(string("--field_mat"),string(""),string("Name of rigid body transform for susceptibility field"),false,Utilities::requires_argument),
  // Paramaters related to the general running of vanilla eddy
  _niter_tmp(string("--niter"),5,string("\tNumber of iterations (default 5)"),false,Utilities::requires_argument),
  _fwhm_tmp(string("--fwhm"),std::vector<float>(1,0.0),string("\tFWHM for conditioning filter when estimating the parameters (default 0)"),false,Utilities::requires_argument),
  _ref_scan_no(string("--ref_scan_no"),0,string("Zero-offset # of ref (for location) volume (default 0)"),false,Utilities::requires_argument),
  _shape_ref_scan_nos(string("--shape_ref_scan_nos"),std::vector<int>(0),string("Zero-offset #s of refs (for shape) volumes"),false,Utilities::requires_argument),
  _interp(string("--interp"),string("spline"),string("Interpolation model for estimation step (spline/trilinear, default spline)"),false,Utilities::requires_argument),
  _extrap(string("--extrap"),string("periodic"),string("Extrapolation model for estimation step (periodic/mirror, default periodic)"),false,Utilities::requires_argument),
  _epvalid(string("--epvalid"),false,string("Indicates that extrapolation is valid in EP direction (default false)"),false, Utilities::no_argument),
  // Parameters related to outlier detection/replacement
  _rep_ol(string("--repol"),false,string("\tDetect and replace outlier slices (default false))"),false, Utilities::no_argument),
  _ol_nstd(string("--ol_nstd"),4.0,string("Number of std off to qualify as outlier (default 4)"),false,Utilities::requires_argument),
  _ol_nvox(string("--ol_nvox"),250,string("Min # of voxels in a slice for inclusion in outlier detection (default 250)"),false,Utilities::requires_argument),
  _ol_ec(string("--ol_ec"),1,string("\tError type (1 or 2) to keep constant for outlier detection (default 1)"),false,Utilities::requires_argument),
  _ol_type(string("--ol_type"),string("sw"),string("Type of outliers, slicewise (sw), groupwise (gw) or both (both). (default sw)"),false,Utilities::requires_argument),
  _ol_jacut(string("--ol_jacut"),3.0,string("Upper threshold for Jacobian outlier detection mask (default 3)"),false,Utilities::requires_argument),
  // Parameters related to slice-to-volume motion correction
  _mb(string("--mb"),1,string("\tMulti-band factor"),false,Utilities::requires_argument),
  _slspec(string("--slspec"),string(""),string("Name of text file completely specifying slice/group acuistion. N.B. --slspec and --json are mutually exclusive."),false,Utilities::requires_argument),
  _json(string("--json"),string(""),string("\tName of .json text file with information about slice timing. N.B. --json and --slspec are mutually exclusive."),false,Utilities::requires_argument),
  _mp_order(string("--mporder"),std::vector<int>(1,0),string("Order of slice-to-vol movement model (default 0, i.e. vol-to-vol"),false,Utilities::requires_argument),
  _s2v_lambda(string("--s2v_lambda"),std::vector<float>(1,1.0),string("Regularisation weight for slice-to-vol movement. (default 1, reasonable range 1--10"),false,Utilities::requires_argument),
  _s2v_niter(string("--s2v_niter"),std::vector<int>(1,5),string("Number of iterations for slice-to-vol (default 5)"),false,Utilities::requires_argument),
  _s2v_fwhm(string("--s2v_fwhm"),std::vector<float>(1,0.0),string("FWHM for conditioning filter when estimating slice-to-vol parameters (default 0)"),false,Utilities::requires_argument),
  _s2v_interp(string("--s2v_interp"),string("trilinear"),string("Slice-to-vol interpolation model for estimation step (spline/trilinear, default trilinear)"),false,Utilities::requires_argument),
  // Parameters related to susceptibility-by-movement interaction estimation
  _estimate_mbs(string("--estimate_move_by_susceptibility"),false,string("Estimate how susceptibility field changes with subject movement (default false)"),false,Utilities::no_argument),
  _mbs_niter(string("--mbs_niter"),10,string("Number of iterations for MBS estimation (default 10)"),false,Utilities::requires_argument),
  _mbs_lambda(string("--mbs_lambda"),10.0,string("Weighting of regularisation for MBS estimation (default 10)"),false,Utilities::requires_argument),
  _mbs_ksp(string("--mbs_ksp"),10.0,string("Knot-spacing for MBS field estimation (default 10mm)"),false,Utilities::requires_argument),
  // Parameters related to the Gaussian process
  _hyparcostfunc(string("--hpcf"),string("CV"),string("\tCost-function for GP hyperparameters (MML/CV/GPP/CC, default CV)"),false, Utilities::requires_argument),
  _nvoxhp(string("--nvoxhp"),1000,string("# of voxels used to estimate the hyperparameters (default 1000)"),false, Utilities::requires_argument),
  _initrand(string("--initrand"),0,string("Seeds rand for when selecting voxels (default 0=no seeding)"),false, Utilities::optional_argument),
  _hyparfudgefactor(string("--ff"),10.0,string("\tFudge factor for hyperparameter error variance (default 10.0)"),false, Utilities::requires_argument),
  _hypar(string("--hypar"),std::vector<float>(0),string("\tUser specified values for GP hyperparameters"),false,Utilities::requires_argument),
  // Miscallenous parameters related to missing planes/output masking
  _fep(string("--fep"),false,string("\tFill empty planes in x- or y-directions (default false)"),false, Utilities::no_argument),
  _dont_mask_output(string("--dont_mask_output"),false,string("Do not mask output to include only voxels present for all volumes (default false)"),false, Utilities::no_argument),
  // Parameters related to initialisation of movement/distortion parameters
  _init(string("--init"),string(""),string("\tText file with parameters for initialisation"),false,Utilities::requires_argument),
  _init_s2v(string("--init_s2v"),string(""),string("Text file with parameters for initialisation of slice-to-vol movement"),false,Utilities::requires_argument),
  _init_mbs(string("--init_mbs"),string(""),string("4D image file for initialisation of movement-by-susceptibility"),false,Utilities::requires_argument),
  // Parameters related to the writing of various "extra" output
  _dfields(string("--dfields"),false,string("Write total displacement fields (default false)"),false, Utilities::no_argument),
  _residuals(string("--residuals"),false,string("Write residuals (between GP and observations), (default false)"),false, Utilities::no_argument),
  _with_outliers(string("--with_outliers"),false,string("Write corrected data (additionally) with outliers retained (default false)"),false, Utilities::no_argument),
  _history(string("--history"),false,string("Write history of mss and parameter estimates (default false)"),false, Utilities::no_argument),
  _write_slice_stats(string("--wss"),false,string("\tWrite slice-wise stats for each iteration (default false)"),false,Utilities::no_argument),
  _write_predictions(string("--write_predictions"),false,string("Write predicted data (in addition to corrected, default false)"),false, Utilities::no_argument),
  _no_text_files(string("--no_textout"),false,string("Supress writing of separate text output files (default false)"),false, Utilities::no_argument),
  // Parameters related to things that can be useful for debugging/development
  _debug_tmp(string("--debug"),0,string("\tLevel of debug print-outs (default 0)"),false,Utilities::requires_argument),
  _dbg_indx_str(string("--dbgindx"),string(""),string("Indicies (zero offset) of volumes for debug print-outs"),false,Utilities::requires_argument),
  _log_timings(string("--log_timings"),false,string("Write timing information (defaut false)"),false, Utilities::no_argument),
  // Parameters related to architecture/threads/cores etc
  _nthr(string("--nthr"),1,string("# of threads used when running eddy (default 1)"),false, Utilities::requires_argument),
  // Finally some generally helpful parameters
  _very_verbose(string("--very_verbose"),false,string("Switch on detailed diagnostic messages (default false)"),false, Utilities::no_argument),
  _verbose(string("-v,--verbose"),false,string("switch on diagnostic messages"),false, Utilities::no_argument),
  _help(string("-h,--help"),false,string("display this message"),false, Utilities::no_argument),
  _is_fmri(false), _init_rand(false), _fixed_hpar(false), _shape_references_set(false)  
{
  // All the commmon parameters have been initialised. Now we need to initialise also the modality specific ones.
  if (argc > 1) {
    if (std::string(argv[1]) != "fmri") {
      if (std::string(argv[1]) != "diffusion") std::cout << std::endl << "Warning: In a future release the first argument will have to be \"diffusion\" when using eddy on diffusion data, i.e." << std::endl << "eddy diffusion --imain='my_ima' --acqp='my_acqp' ..." << std::endl << std::endl;
      _is_fmri = false;
    }
    else { // This means that the first argument is "fmri"
      _is_fmri = true;
      // Temporarily block the fMRI branch to allow for bugfix release.
      throw EddyInputError("fMRI functionality is currently disabled awaiting more testing");
    }
  }
  else { // This means that someone only typed eddy
    std::cout << std::endl << "Warning: In a future release the first argument will have to be \"diffusion\" when using eddy on diffusion data, i.e." << std::endl << "eddy diffusion --imain='my_ima' --acqp='my_acqp' ..." << std::endl << std::endl;
    std::cout << "If you want to see the help screen for eddy for fmri, just type \"eddy fmri\"" << std::endl << std::endl;
  }

  do_initial_parsing(argc,argv);

  // Do sanity checking and some initial reading of files
  // First warn about use of deprecated parameters
  if (IsDiffusion()) {
    if (_diff._sep_offs_move.set()) std::cout << "Warning: --sep_offs_move parameter is deprecated and will crash a future version" << std::endl << std::endl;
    if (_diff._peas.set()) std::cout << "Warning: --peas parameter is deprecated and will crash a future version" << std::endl << std::endl;
    if (_diff._rms.set()) std::cout << "Warning: --rms parameter is deprecated and will crash a future version" << std::endl << std::endl;
    if (_diff._mb_offs.set()) std::cout << "Warning: --mb_offs parameter is deprecated and will crash a future version" << std::endl << std::endl;
    if (_diff._slorder.set()) throw EddyInputError("--slorder has been deprecated");
  }
  if (!_no_text_files.value()) std::cout << "Warning: Writing of individual text files will be discontinued in favour of a single .json file in future versions" << std::endl << std::endl;

  // Next, check on all common (to diffusion and fmri) parameters.
  // I have tried to keep the order of the checks to be the same as the order in which they are added to options.

  // First mandatory parameters
  // Image file exists?
  NEWIMAGE::volume4D<float> imahdr;
  try { NEWIMAGE::read_volume_hdr_only(imahdr,_imain.value()); }  // Throws if there is a problem
  catch(...) { throw EddyInputError("Error when attempting to read --imain file " + _imain.value()); }
  NEWIMAGE::volume<float> maskhdr;
  try { NEWIMAGE::read_volume_hdr_only(maskhdr,_mask.value()); }  // Throws if there is a problem
  catch(...) { throw EddyInputError("Error when attempting to read --mask file " + _mask.value()); }
  if (!samesize(imahdr[0],maskhdr,3,true)) throw EddyInputError("--imain and --mask images must have the same dimensions");
  if (maskhdr.tsize() != 1) throw EddyInputError("--mask image must be 3D volume");
  // File with acquisition parameters exists and makes sense
  NEWMAT::Matrix acqpM;
  try {
    acqpM = MISCMATHS::read_ascii_matrix(_acqp.value());
    if (acqpM.Ncols() != 4) throw EddyInputError("Each row of the --acqp file must contain 4 values");
  }
  catch (...) { throw EddyInputError("Error when attempting to read --acqp file " + _acqp.value()); }
  // Index file exists and makes sense given image file
  NEWMAT::Matrix index;
  try {
    index = MISCMATHS::read_ascii_matrix(_index.value());
    if (std::min(index.Nrows(),index.Ncols()) != 1 || std::max(index.Nrows(),index.Ncols()) != imahdr.tsize()) {
      throw EddyInputError("--index must be an 1xN or Nx1 matrix where N is the number of volumes in --imain");
    }
  }
  catch (...) { throw EddyInputError("Error when attempting to read --index file " + _index.value()); }
  if (!this->indicies_kosher(index,acqpM)) throw EddyInputError("Mismatch between --index and --acqp");
  if (_out.value() == string("")) throw EddyInputError("--out must be a non-empty string");

  // Parameters related to input susceptibility field
  // Check that not both topup and field files have been specified
  if (_topup.set() && _topup.value() != string("") && _field.set() && _field.value() != string("")) {
    throw EddyInputError("One cannot specify both --field and --topup file");
  }
  // Topup file exists if one has been specified
  if (_topup.set() && _topup.value() != string("")) {
    try { TOPUP::TopupFileReader  tfr(_topup.value()); }
    catch (const TOPUP::TopupFileIOException& e) { cout << e.what() << endl; throw EddyInputError("Error when attempting to read --topup files " + _topup.value() + "_fieldcoef.nii.gz and " + _topup.value() + "_movpar.txt"); }
    catch (...) { throw EddyInputError("Error when attempting to read --topup files " + _topup.value() + "_fieldcoef.nii.gz and " + _topup.value() + "_movpar.txt"); }
  }
  // Field file exists if one has been specified
  if (_field.set() && _field.value() != string("")) {
    try { TOPUP::TopupFileReader  tfr(_field.value()); }
    catch (const TOPUP::TopupFileIOException& e) { cout << e.what() << endl; throw EddyInputError("Error when attempting to read --field file " + _field.value()); }
    catch (...) { throw EddyInputError("Error when attempting to read --field file " + _field.value()); }
  }
  // Field mat file exists if one has been specified
  if (_field_mat.set() && _field_mat.value() != string("")) {
    try {
      NEWMAT::Matrix fieldM = MISCMATHS::read_ascii_matrix(_field_mat.value());
      if (fieldM.Nrows() != 4 || fieldM.Ncols() != 4) throw EddyInputError("--field_mat must be a 4x4 matrix");
      NEWMAT::Matrix ul = fieldM.SubMatrix(1,3,1,3);
      float det = ul.Determinant();
      if (std::abs(det-1.0) > 1e-6) throw EddyInputError("--field_mat must be a rigid body transformation");
    }
    catch (...) { throw EddyInputError("Error when attempting to read --field_mat file" + _field_mat.value()); }
  }

  // Paramaters related to the general running of vanilla eddy
  // FWHM either once and for all or one per iteration
  _niter = static_cast<unsigned int>(_niter_tmp.value());
  _fwhm = _fwhm_tmp.value();
  if (_fwhm.size() != 1 && _fwhm.size() != _niter) throw EddyInputError("--fwhm must be one value or one per iteration");
  if (_fwhm.size() != _niter) _fwhm.resize(_niter,_fwhm[0]);
  // Reasonable FWHM
  for (unsigned int i=0; i<_fwhm.size(); i++) {
    if (_fwhm[i] < 0.0 || _fwhm[i] > 20.0) throw EddyInputError("--fwhm value outside valid range 0-20mm");
  }
  // Location ref value is reasonable.
  if (_ref_scan_no.value() < 0 || _ref_scan_no.value() > imahdr.tsize()-1) throw EddyInputError("--ref_scan_no out of range");
  // Shape reference reasonable if set. Must have exactly one ref per shell. 
  if (_shape_ref_scan_nos.set()) { 
    if (IsfMRI()) { // fMRI counts as a single shell
      if (_shape_ref_scan_nos.value().size() != 1) throw EddyInputError("--shape_ref_scan_nos has wrong number of entries");
      if (_shape_ref_scan_nos.value()[0] < 0 || _shape_ref_scan_nos.value()[0] > imahdr.tsize()-1) throw EddyInputError("--shape_ref_scan_nos out of range");
      _fmri._shape_reference = _shape_ref_scan_nos.value()[0];
      _shape_references_set = true;
    }
    else {  // It is diffusion, so we must check it against shells
      NEWMAT::Matrix bvalsM;
      try { // First read bvals file so we know what the shells are
	bvalsM = MISCMATHS::read_ascii_matrix(_diff._bvals.value());
	if (std::min(bvalsM.Nrows(),bvalsM.Ncols()) != 1 || std::max(bvalsM.Nrows(),bvalsM.Ncols()) != imahdr.tsize()) {
	  throw EddyInputError("--bvals should contain a 1xN or Nx1 matrix where N is the number of volumes in --imain");
	}
      }
      catch (...) { throw EddyInputError("Error when attempting to read --bvals file " + _diff._bvals.value()); }
      // Now create a vector of DiffPara objects so we can use existing functions
      if (bvalsM.Nrows() < bvalsM.Ncols()) bvalsM = bvalsM.t();
      std::vector<EDDY::DiffPara> dpvec;
      for (int i=0; i<bvalsM.Nrows(); i++) dpvec.push_back(EDDY::DiffPara(bvalsM(i+1,1)));
      std::vector<double> shell_b_values;
      std::vector<unsigned int> shell_indicies;
      double original_b_range = EddyUtils::GetbRange(); 
      EddyUtils::SetbRange(_diff._brange.value());      // Temporarily set b-range to value we will check later
      if (!EddyUtils::GetGroups(dpvec,shell_indicies,shell_b_values) && !_diff._data_is_shelled.value()) throw EddyInputError("Data not shelled");
      EddyUtils::SetbRange(original_b_range);
      if (shell_b_values.size() != _shape_ref_scan_nos.value().size()) throw EddyInputError("There must be one index per shell in --shape_ref_scan_nos");
      std::vector<int> shells(_shape_ref_scan_nos.value().size()); // The shell index for each of the shape-refs
      for (unsigned int i=0; i<_shape_ref_scan_nos.value().size(); i++) shells[i] = shell_indicies[_shape_ref_scan_nos.value()[i]];
      for (unsigned int i=0; i<shells.size(); i++) {
	for (unsigned int j=i+1; j<shells.size(); j++) {
	  if (shells[i]==shells[j]) throw EddyInputError("There must be one index per shell in --shape_ref_scan_nos");
	}
      }
      // If we are here, it means all is fine with --shape_ref_scan_nos
      this->set_shape_refs(_shape_ref_scan_nos.value(),shells,shell_b_values);
    }
  }
  // Valid interpolation method?
  if (_interp.value() != string("spline") && _interp.value() != string("trilinear")) throw EddyInputError("Invalid --interp parameter");
  // Valid extrapolation method?
  if (_extrap.value() != string("mirror") && _extrap.value() != string("periodic")) throw EddyInputError("Invalid --extrap parameter");
  // Correct extrapolation method for valid extrapolation
  if (_extrap.value() == string("mirror") && _epvalid.value()) throw EddyInputError("--extrap=mirror cannot be combined with --epvalid");

  // Parameters related to outlier detection/replacement
  // Reasonable values for outlier detection
  if (_ol_nstd.value() < 1.96) throw EddyInputError("--ol_nstd value too low (below 1.96)");
  if (_ol_nvox.value() < 250) throw EddyInputError("--ol_nvox value below 250");
  if (IsDiffusion()) _oldef = EDDY::OutlierDefinition(_ol_nstd.value(),_ol_nvox.value(),_diff._ol_pos.value(),_diff._ol_sqr.value());
  else throw EddyInputError("Must look into outlier definition for eddy-fmri");
  if (_ol_ec.value() != 1 && _ol_ec.value() != 2) throw EddyInputError("--ol_ec value must be 1 or 2");
  if (_ol_type.value() != string("sw") && _ol_type.value() != string("gw") && _ol_type.value() != string("both")) throw EddyInputError("--ol_type must be \"sw\", \"gw\" or \"both\"");
  if ((_ol_type.value() == string("gw") || _ol_type.value() == string("both")) && MultiBand().MBFactor() < 2) throw EddyInputError("--ol_type indicating mb-groups without providing mb structure");
  if (_ol_jacut.value() < 0.0) throw EddyInputError("--ol_jacut must be a positive number");

  // Parameters related to slice-to-volume motion correction
  // Make sure that slices only specified in one way (i.e. that only one of --mb, --slspec and --json has been used)
  if ((_mb.set() + _slspec.set() + _json.set()) > 1) throw EddyInputError("--mb, --slspec and --json mutually exclusive");
  // Reasonable values for mb factor and offset
  if (_mb.value() < 1 || _mb.value() > 10) throw EddyInputError("--mb value outside valid range 1-10");
  if (std::abs(_diff._mb_offs.value()) > 1) throw EddyInputError("--mb_offs must be -1, 0 or 1");
  // MB offset only makes sense in context of MB
  if (!_mb.set() && _diff._mb_offs.set()) throw EddyInputError("--mb_offs makes no sense without --mb");
  // Check for valid multi-band structure
  EDDY::MultiBandGroups mbg(imahdr.zsize());
  if (_mb.set() || _slspec.set() || _json.set()) {
    try {
      mbg = this->MultiBand();
    }
    catch (const std::exception& e) { throw EddyInputError("Failed to decode valid multi-band structure. Exception thrown with message: " + std::string(e.what())); }
    if (mbg.NSlices() != imahdr.zsize()) throw EddyInputError("Mismatch between no of slices in multi-band structure and in image file.");
  }
  // Check that movement model order not greater than number of groups
  if (*std::max_element(_mp_order.value().begin(),_mp_order.value().end()) > mbg.NGroups()) throw EddyInputError("--mporder can not be greater than number of slices/mb-groups");
  // Valid slice-to-vol interpolation method?
  if (_s2v_interp.value() != string("spline") && _s2v_interp.value() != string("trilinear")) throw EddyInputError("Invalid --s2v_interp parameter");
  // Valid and reasonable parameters for slice-to-vol (throws with error message from within constructor)
  _s2vparam = EDDY::S2VParam(_mp_order.value(),_s2v_lambda.value(),_s2v_fwhm.value(),_s2v_niter.value());

  // Parameters related to susceptibility-by-movement interaction
  // Check for reasonable values pertaining to MBS (Movement By Susceptibility).
  if (_init_mbs.set() && _init_mbs.value() != std::string("")) {
    if (_mbs_niter.value() < 0 || _mbs_niter.value() > 50) throw EddyInputError("--mbs_niter value outside valid range 0-50");
  }
  else if (_mbs_niter.value() < 1 || _mbs_niter.value() > 50) throw EddyInputError("--mbs_niter value outside valid range 1-50");
  if (_mbs_lambda.value() < 1.0 || _mbs_lambda.value() > 10000.0) throw EddyInputError("--mbs_lambda value outside valid range 1.0-10000.0");
  if (_mbs_ksp.value() < 2.0 || _mbs_ksp.value() > 30.0) throw EddyInputError("--mbs_ksp value outside valid range 2.0-30.0");

  // Parameters related to the Gaussian process
  // Valid hypar cost-function?
  if (_hyparcostfunc.value() != string("MML") && _hyparcostfunc.value() != string("CV") && _hyparcostfunc.value() != string("GPP") && _hyparcostfunc.value() != string("CC")) throw EddyInputError("Invalid --hpcf parameter");
  // Reasonable # of voxels for estimation of hyperparameters
  if (_nvoxhp.value() < 100  || _nvoxhp.value() > 50000) throw EddyInputError("--nvoxhp value outside valid range 100-50000");
  _nvoxhp_internal = static_cast<unsigned int>(_nvoxhp.value());
  // Mimic old behaviour by setting srand to one if --initrand is set without an explicit value
  if (_initrand.set()) {
    if (_initrand.value() == 0) _init_rand = 1;
    else _init_rand = _initrand.value();
  }
  else _init_rand = 0;
  // Reasonable fudge factor for hyperparameters
  if (_hyparfudgefactor.value() < 1.0 || _hyparfudgefactor.value() > 10.0) throw EddyInputError("--ff value outside valid range 1.0-10.0");
  else _hypar_ff_internal = static_cast<double>(_hyparfudgefactor.value());
  // User specified hyperparameters?
  if (_hypar.value().size()) { _fixed_hpar = true; _hypar_internal = _hypar.value(); } else { _fixed_hpar = false; }

  // Parameters related to initialisation of movement/distortion parameters
  // Initialisation file exists if one has been specified
  if (_init.set() && _init.value() != string("")) {
    NEWMAT::Matrix  tmp = MISCMATHS::read_ascii_matrix(_init.value());
    if (tmp.Nrows() == 0) throw EddyInputError("Error when attempting to read --init file " + _init.value());
  }
  // Initialisation for s2v not set if EC/volumetric initialisation not set
  if (_init_s2v.set() && _init_s2v.value() != string("") && (!_init.set() || _init.value() == string(""))) {
    throw EddyInputError("--init_s2v can not be used without --init");
  }
  // Check that s2v initialisation file exists if one has been specified
  if (_init_s2v.set() && _init_s2v.value() != string("")) {
    NEWMAT::Matrix  tmp = MISCMATHS::read_ascii_matrix(_init_s2v.value());
    if (tmp.Nrows() == 0) throw EddyInputError("Error when attempting to read --init_s2v file " + _init_s2v.value());
  }
  // Check that we are not to run any vol-to-vol iterations if s2v initialisation is set
  if (_init_s2v.set() && _init_s2v.value() != string("") && _niter != 0) {
    throw EddyInputError("--init_s2v does not make sense unless --niter=0");
  }
  // Movement-by-susceptibility file exists and has correct dimensions
  if (_init_mbs.set() && _init_mbs.value() != string("")) {
    NEWIMAGE::volume4D<float> mbshdr;
    try { NEWIMAGE::read_volume_hdr_only(mbshdr,_init_mbs.value()); }  // Throws if there is a problem
    catch(...) { throw EddyInputError("Error when attempting to read --init_mbs file " + _init_mbs.value()); }
    if (!samesize(imahdr[0],mbshdr,3,true)) throw EddyInputError("--imain and --init_mbs images must have the same dimensions");
    if (mbshdr.tsize() != int(this->MoveBySuscParam().size())) throw EddyInputError("--init_mbs has wrong number of volumes");
  }

  // Parameters related to the writing of various "extra" output
  // Check that --with_outliers hasn't been set unless --repol is also set
  if (_with_outliers.value() && !_rep_ol.value()) throw EddyInputError("--with_outliers does not make sense without also specifying --repol");

  // Parameters related to things that can be useful for debugging/development
  // Make sure that debug requests are reasonable
  if (_debug_tmp.value() < 0 || _debug_tmp.value() > 4) { throw EddyInputError("--debug must be a value between 0 and 4"); }
  _debug = static_cast<unsigned int>(_debug_tmp.value());
  if (_debug > 0 && (!_dbg_indx_str.set() || _dbg_indx_str.value() == string(""))) {
    _dbg_indx = DebugIndexClass(0,imahdr.tsize()-1);
  }
  else if (_dbg_indx_str.set() && _dbg_indx_str.value() != string("")) {
    _dbg_indx = DebugIndexClass(_dbg_indx_str.value());
  }
  if (_dbg_indx.Max() > static_cast<unsigned int>(imahdr.tsize()-1)) { throw EddyInputError("--dbg_indx out of range"); }
  // If debug is set we should also set initialisation of rand
  if (_debug > 0) _init_rand = true;

  // Make sure we don't attempt to use multiple CPU threads for GPU version
  #ifdef COMPILE_GPU
  if (_nthr.value() > 1) {
    throw EddyInputError("The version compiled for GPU can only use 1 CPU thread (i.e. --nthr=1)");
  }
  #endif

  // Make sure we don't try to use more CPU threads than the hardware supports
  if (_nthr.value() > std::thread::hardware_concurrency()) {
    if (std::thread::hardware_concurrency()) { // This means that if hardware_concurrency isn't implemented we trust the user.
      _nthr.set_T(std::thread::hardware_concurrency());
    }
  }
  else if (_nthr.value() < 1) _nthr.set_T(1);

  // Finally some generally helpful parameters
  // If very_verbose is set we should also set verbose
  if (_very_verbose.value() && !_verbose.value()) _verbose.set_T(true);

  // The session parameter is mandatory for fMRI, but not for diffusion.
  // That is why there are two different fields for it.
  NEWMAT::Matrix sessM;
  // Check for valid session file
  if (IsfMRI() || _diff._session.set()) { // There is/should be a session file
    try {
      if (IsDiffusion()) sessM = MISCMATHS::read_ascii_matrix(_diff._session.value());
      else sessM = MISCMATHS::read_ascii_matrix(_fmri._session.value());
      if (std::min(sessM.Nrows(),sessM.Ncols()) != 1 || std::max(sessM.Nrows(),sessM.Ncols()) != imahdr.tsize()) {
	throw EddyInputError("--session must be an 1xN or Nx1 matrix where N is the number of volumes in --imain");
      }
    }
    catch (...) { 
      if (IsDiffusion()) throw EddyInputError("Error when attempting to read --session file " + _diff._session.value()); 
      else throw EddyInputError("Error when attempting to read --session file " + _fmri._session.value()); 
    }
    if (!session_indicies_kosher(sessM)) throw EddyInputError("Values in --session file seems dodgy");
  }
  else { // There is no (and should not be) a session file.
    _sessvec.resize(imahdr.tsize(),1); // Set single session
  }
  
  // This ends the checking of common parameters
  // Next, check on diffusion parameters.

  if (IsDiffusion()) {
    // First mandatory parameters
    // bvals and bvecs exists and are consistent with image file
    try {
      NEWMAT::Matrix bvecsM = MISCMATHS::read_ascii_matrix(_diff._bvecs.value());
      if (std::min(bvecsM.Nrows(),bvecsM.Ncols()) != 3 || std::max(bvecsM.Nrows(),bvecsM.Ncols()) != imahdr.tsize()) {
	throw EddyInputError("--bvecs should contain a 3xN or Nx3 matrix where N is the number of volumes in --imain");
      }
    }
    catch (...) { throw EddyInputError("Error when attempting to read --bvecs file " + _diff._bvecs.value()); }
    try {
      NEWMAT::Matrix bvalsM = MISCMATHS::read_ascii_matrix(_diff._bvals.value());
      if (std::min(bvalsM.Nrows(),bvalsM.Ncols()) != 1 || std::max(bvalsM.Nrows(),bvalsM.Ncols()) != imahdr.tsize()) {
	throw EddyInputError("--bvals should contain a 1xN or Nx1 matrix where N is the number of volumes in --imain");
      }
    }
    catch (...) { throw EddyInputError("Error when attempting to read --bvals file " + _diff._bvals.value()); }
    if (_diff._bdeltas.value() != std::string("")) {
      try {
	NEWMAT::Matrix bdeltasM = MISCMATHS::read_ascii_matrix(_diff._bdeltas.value());
	if (bdeltasM.Nrows()>bdeltasM.Ncols()) bdeltasM = bdeltasM.t();
	if (bdeltasM.Nrows()!=1 || bdeltasM.Ncols()!=imahdr.tsize()) {
	  throw EddyInputError("--bdeltas should contain a 1xN or Nx1 matrix where N is the number of volumes in --imain");
	}
	for (int i=0; i<bdeltasM.Ncols(); i++) {
	  if (std::fabs(bdeltasM(1,i+1)-1.0)>1e-4 && std::fabs(bdeltasM(1,i+1)+0.5)>1e-4 && std::fabs(bdeltasM(1,i+1))>1e-4) {
	    throw EddyInputError("--bdeltas should only contain values 1.0, -0.5 or 0.0");
	  }
	}
      }
      catch (...) { throw EddyInputError("Error when attempting to read --bdeltas file " + _diff._bdeltas.value()); }
    }

    // Paramaters related to the general running of vanilla eddy
    // Valid range of b-values within a single shell
    if (_diff._brange.value() < 0.0 || _diff._brange.value() > 200.0) throw EddyInputError("Invalid --brange parameter");
    // Valid final resampling method?
    if (_diff._resamp.value() != string("jac") && _diff._resamp.value() != string("jac_nn") && _diff._resamp.value() != string("lsr")) {
      throw EddyInputError("Invalid --resamp parameter");
    }
    // Check that output masking is not off for LSR resampling
    if (_diff._resamp.value() == string("lsr") && _dont_mask_output.value()) throw EddyInputError("You cannot combine --resamp=lsr with --dont_mask_output");
    // Reasonable LSR regularisation weighting
    if (_diff._lsr_lambda.value() < 0.0 || _diff._lsr_lambda.value() > 1.0) throw EddyInputError("--lsr_lambda value outside valid range");

    // Parameters related to eddy current modelling
    // Valid first level model?
    if (_diff._flm.value() != string("movement") && _diff._flm.value() != string("linear") && _diff._flm.value() != string("quadratic") && _diff._flm.value() != string("cubic")) throw EddyInputError("Invalid --flm parameter");
    // Valid second level model?
    if (_diff._slm_str.value() == string("none")) _diff._slm = EDDY::SecondLevelECModelType::None;
    else if (_diff._slm_str.value() == string("linear")) _diff._slm = EDDY::SecondLevelECModelType::Linear;
    else if (_diff._slm_str.value() == string("quadratic")) _diff._slm = EDDY::SecondLevelECModelType::Quadratic;
    else throw EddyInputError("Invalid --slm parameter");
    // Valid first level b0-model?
    if (_diff._b0_flm.value() != string("movement") && _diff._b0_flm.value() != string("linear") && _diff._b0_flm.value() != string("quadratic")) throw EddyInputError("Invalid --b0_flm parameter");
    // Valid second level b0 model?
    if (_diff._b0_slm_str.value() == string("none")) _diff._b0_slm = EDDY::SecondLevelECModelType::None;
    else if (_diff._b0_slm_str.value() == string("linear")) _diff._b0_slm = EDDY::SecondLevelECModelType::Linear;
    else if (_diff._b0_slm_str.value() == string("quadratic")) _diff._b0_slm = EDDY::SecondLevelECModelType::Quadratic;
    else throw EddyInputError("Invalid --b0_slm parameter");
    // Valid model for long (time varying) EC?
    if (_diff._long_ec_str.value() == string("none")) _diff._long_ec_model = EDDY::LongECModelType::None;
    else if (_diff._long_ec_str.value() == string("weights")) _diff._long_ec_model = EDDY::LongECModelType::Individual;
    else if (_diff._long_ec_str.value() == string("tc")) _diff._long_ec_model = EDDY::LongECModelType::IndividualTimeConstant;
    else if (_diff._long_ec_str.value() == string("jointweights")) _diff._long_ec_model = EDDY::LongECModelType::Joint;
    else if (_diff._long_ec_str.value() == string("jointtc")) _diff._long_ec_model = EDDY::LongECModelType::JointTimeConstant;
    else throw EddyInputError("Invalid --lecm parameter");
    if (_diff._long_ec_model != EDDY::LongECModelType::None && mbg.MBFactor() < 2) {
      throw EddyInputError("--lecm can only be used for multi-band data");
    }
    if (_diff._long_ec_str.value() == string("none") && _diff._long_ec_dont_reest.value() == true) {
      throw EddyInputError("--lecm_dont_reest only makes sense together with --lecm");
    }
    if (_diff._long_ec_str.value() == string("none") && _diff._long_ec_sep_offs_move.value() == true) {
      throw EddyInputError("--lecm_sep_offs_move only makes sense together with --lecm");
    }
    
    // Parameters related to outlier detection/replacement
    if (_diff._ol_ss.value() != string("sw") && _diff._ol_ss.value() != string("pooled")) throw EddyInputError("--ol_ss must be \"sw\" (shell-wise) or \"pooled\"");

    // Parameters related to the Gaussian process
    // Valid covariance function?
    if (_diff._covfunc.value() != string("spheri") && _diff._covfunc.value() != string("expo") && _diff._covfunc.value() != string("old")) throw EddyInputError("Invalid --covfunc parameter");

    // Miscallenous parameters related to missing planes/output masking/bvec-rotation
    // Set "internal" copy of _rbvde to allow us to set/reset internally
    _diff._rbvde_internal = _diff._rbvde.value();

    // Parameters related to "tricky" movements for diffusion
    // Valid offset model?
    if (_diff._offset_model.value() != string("linear") && _diff._offset_model.value() != string("linear_lag") &&
	_diff._offset_model.value() != string("quadratic") && _diff._offset_model.value() != string("linear_lag_quadratic") &&
	_diff._offset_model.value() != string("linear_lag_quadratic_lag")) throw EddyInputError("Invalid --offset_model parameter");
    // Set internal parameter for "dont_sep_offs_move"
    _diff._dont_sep_offs_move_internal = _diff._dont_sep_offs_move.value();
    // Check that peas options are compatible
    if (_diff._dont_peas.value() && _diff._use_b0s_for_peas.value()) throw EddyInputError("--dont_peas and --b0_peas cannot both be set");

    // Parameters related to writing of extra output
    // Check that --write_scatter_brain_predictions hasn't been specified for volume-to-volume registration
    if (_diff._write_scatter_brain_predictions.value() && !IsSliceToVol()) throw EddyInputError("--write_scatter_brain_predictions does not make sense without also specifying --mporder > 0");
    // Check that --write_scatter_brain_predictions hasn't been specified without --write_predictions
    if (_diff._write_scatter_brain_predictions.value() && !_write_predictions.value()) throw EddyInputError("--write_scatter_brain_predictions can only be used together with --write_predictions");

    // Parameters related to things that can be useful for debugging/development
    // Make sure the "only" flags are mutually exclusive
    _diff._rb0 = true; _diff._rdwi = true;
    if (_diff._dwi_only.value() && _diff._b0_only.value()) throw EddyInputError("--dwi_only and --b0_only cannot both be set");
    if (_diff._dwi_only.value()) _diff._rb0 = false;
    else if (_diff._b0_only.value()) _diff._rdwi = false;
    // Make sure that test_rot vector is valid if set
    if (_diff._test_rot.set()) {
      if (_diff._test_rot.value().size() < 1 || _diff._test_rot.value().size() > 3) throw EddyInputError("--test_rot must be one, two or three angles");
    }
  } // Here ends the checking of diffusion parameters
  else if (IsfMRI()) { // Here starts checking of parameters specific for fMRI
    // Check for reasonable repetition time
    if (_fmri._rep_time.value() < 0.5 || _fmri._rep_time.value() > 10.0) throw EddyInputError("--rep_time must be between 0.5 and 10 seconds");
    // Check for reasonable echo time
    if (_fmri._echo_time.set() && (_fmri._echo_time.value() < 10.0 || _fmri._echo_time.value() > 100.0)) throw EddyInputError("--echo_time must be between 10 and 100 milliseconds");
  }
}
catch (const EddyInputError& e)
{
  cout << e.what() << endl;
  cout << "Terminating program" << endl;
  exit(EXIT_FAILURE);
} EddyCatch

EDDY::FinalResamplingType EddyCommandLineOptions::ResamplingMethod() const EddyTry
{
  if (IsfMRI() || _diff._resamp.value() == std::string("jac")) return(FinalResamplingType::Jac);
  else if (_diff._resamp.value() == std::string("jac_nn")) return(FinalResamplingType::Jac_NN);
  else if (_diff._resamp.value() == std::string("lsr")) return(FinalResamplingType::LSR);
  else return(FinalResamplingType::Unknown);
} EddyCatch

EDDY::CovarianceFunctionType EddyCommandLineOptions::CovarianceFunction() const EddyTry
{
  if (IsDiffusion()) {
    if (_diff._covfunc.value() == std::string("spheri")) return(EDDY::CovarianceFunctionType::NewSpherical);
    else if (_diff._covfunc.value() == std::string("expo")) return(EDDY::CovarianceFunctionType::Exponential);
    else if (_diff._covfunc.value() == std::string("old")) return(EDDY::CovarianceFunctionType::Spherical);
    else return(EDDY::CovarianceFunctionType::Unknown);
  }
  else if (IsfMRI()) {
    if (_fmri._covfunc.value() == std::string("sqexp")) return(EDDY::CovarianceFunctionType::SquaredExponential);
    else return(EDDY::CovarianceFunctionType::Unknown);
  }
  else return(EDDY::CovarianceFunctionType::Unknown);
} EddyCatch

EDDY::HyParCostFunctionType EddyCommandLineOptions::HyParCostFunction() const EddyTry
{
  if (_hyparcostfunc.value() == std::string("MML")) return(EDDY::HyParCostFunctionType::MML);
  else if (_hyparcostfunc.value() == std::string("CV")) return(EDDY::HyParCostFunctionType::CV);
  else if (_hyparcostfunc.value() == std::string("GPP")) return(EDDY::HyParCostFunctionType::GPP);
  else if (_hyparcostfunc.value() == std::string("CC")) return(EDDY::HyParCostFunctionType::CC);
  else return(EDDY::HyParCostFunctionType::Unknown);
} EddyCatch

float EddyCommandLineOptions::FWHM(unsigned int iter) const EddyTry
{
  if (_fwhm.size()==1) return(_fwhm[0]);
  else if (iter >= _fwhm.size()) throw EddyException("EddyCommandLineOptions::FWHM: iter out of range");
  return(_fwhm[iter]);
} EddyCatch

void EddyCommandLineOptions::SetNIterAndFWHM(unsigned int niter, const std::vector<float>& fwhm) EddyTry
{
  if (fwhm.size() != 1 && fwhm.size() != niter) throw EddyException("EddyCommandLineOptions::SetNIterAndFWHM: mismatch between niter and fwhm");
  _niter = niter;
  if (fwhm.size() == niter) _fwhm = fwhm;
  else _fwhm = std::vector<float>(niter,fwhm[0]);
} EddyCatch

void EddyCommandLineOptions::SetS2VParam(unsigned int order, float lambda, float fwhm, unsigned int niter) EddyTry
{
  std::vector<int>   lorder(1,static_cast<int>(order));
  std::vector<float> llambda(1,lambda);
  std::vector<float> lfwhm(1,fwhm);
  std::vector<int>   lniter(1,static_cast<int>(niter));

  _s2vparam = EDDY::S2VParam(lorder,llambda,lfwhm,lniter);
} EddyCatch

unsigned int S2VParam::NIter(unsigned int  oi) const EddyTry
{
  if (oi >= _niter.size()) throw EddyException("S2VParam::NIter: oi out of range");
  return(static_cast<unsigned int>(_niter[oi]));
} EddyCatch

unsigned int S2VParam::NOrder() const EddyTry
{
  if (_order.size() > 1) return(_order.size());
  else return(static_cast<unsigned int>(_order[0] > 0 ? 1 : 0));
} EddyCatch

unsigned int S2VParam::Order(unsigned int oi) const EddyTry
{
  if (oi >= _order.size()) throw EddyException("S2VParam::Order: oi out of range");
  return(static_cast<unsigned int>(_order[oi]));
} EddyCatch

std::vector<unsigned int> S2VParam::Order() const EddyTry
{
  std::vector<unsigned int> rval(_order.size());
  for (unsigned int i=0; i<_order.size(); i++) rval[i] = static_cast<unsigned int>(_order[i]);
  return(rval);
} EddyCatch

double S2VParam::Lambda(unsigned int oi) const EddyTry
{
  if (oi >= _lambda.size()) throw EddyException("S2VParam::Lambda: oi out of range");
  return(static_cast<double>(_lambda[oi]));
} EddyCatch

std::vector<double> S2VParam::Lambda() const EddyTry
{
  std::vector<double> rval(_lambda.size());
  for (unsigned int i=0; i<_lambda.size(); i++) rval[i] = static_cast<double>(_lambda[i]);
  return(rval);
} EddyCatch

float S2VParam::FWHM(unsigned int oi, unsigned int iter) const EddyTry
{
  if (oi >= this->NOrder()) throw EddyException("S2VParam::FWHM: oi out of range");
  if (iter >= this->NIter(oi)) throw EddyException("S2VParam::FWHM: iter out of range");
  if (_fwhm.size()==1) return(_fwhm[0]);
  else {
    unsigned int indx=0;
    for (unsigned int i=0; i<oi; i++) indx += this->NIter(i);
    return(_fwhm[indx+iter]);
  }
} EddyCatch

std::vector<float> S2VParam::FWHM(unsigned int oi) const EddyTry
{
  if (oi >= NOrder()) throw EddyException("S2VParam::FWHM: oi out of range");
  if (_fwhm.size()==1) {
    std::vector<float> rval(this->NIter(oi),_fwhm[0]);
    return(rval);
  }
  else {
    unsigned int indx=0;
    for (unsigned int i=0; i<oi; i++) indx += this->NIter(i);
    std::vector<float> rval(this->NIter(oi),0.0);
    for (unsigned int i=0; i<this->NIter(oi); i++) rval[i] = _fwhm[indx+i];
    return(rval);
  }
} EddyCatch

EDDY::MultiBandGroups EddyCommandLineOptions::MultiBand() const EddyTry
{
  NEWIMAGE::volume<float> mhdr;
  NEWIMAGE::read_volume_hdr_only(mhdr,_mask.value());
  if (_slspec.set()) {
    EDDY::MultiBandGroups mbg(_slspec.value());
    return(mbg);
  }
  else if (_json.set()) {
    EDDY::JsonReader jr(_json.value());
    EDDY:: MultiBandGroups mbg(jr.SliceOrdering());
    return(mbg);
  }
  else {
    EDDY::MultiBandGroups mbg(mhdr.zsize(),static_cast<unsigned int>(_mb.value()),static_cast<unsigned int>(_diff._mb_offs.value()));
    if (IsDiffusion()) { if (_diff._slorder.set() && _diff._slorder.value() != string("")) mbg.SetTemporalOrder(get_slorder(_diff._slorder.value(),mbg.NGroups())); }
    return(mbg);
  }
} EddyCatch

EDDY::ECModelType EddyCommandLineOptions::FirstLevelModel() const EddyTry
{
  if (IsDiffusion()) {
    if (_diff._flm.value() == string("movement")) return(EDDY::ECModelType::NoEC);
    else if (_diff._flm.value() == string("linear")) return(EDDY::ECModelType::Linear);
    else if (_diff._flm.value() == string("quadratic")) return(EDDY::ECModelType::Quadratic);
    else if (_diff._flm.value() == string("cubic")) return(EDDY::ECModelType::Cubic);
    return(EDDY::ECModelType::Unknown);
  }
  return(EDDY::ECModelType::NoEC); // No EC for fMRI
} EddyCatch

EDDY::ECModelType EddyCommandLineOptions::b0_FirstLevelModel() const EddyTry
{
  if (IsDiffusion()) {
    if (_diff._b0_flm.value() == string("movement")) return(EDDY::ECModelType::NoEC);
    else if (_diff._b0_flm.value() == string("linear")) return(EDDY::ECModelType::Linear);
    else if (_diff._b0_flm.value() == string("quadratic")) return(EDDY::ECModelType::Quadratic);
    return(EDDY::ECModelType::Unknown);
  }
  return(EDDY::ECModelType::NoEC); // No EC for fMRI
} EddyCatch

std::string EddyCommandLineOptions::LongECModelString() const EddyTry
{
  LongECModelType lmt = this->LongECModel();
  switch (lmt) {
  case LongECModelType::None:
    return("No modelling of long EC");
  case LongECModelType::Individual:
    return("Modelling weight for each MB-group of each volume separately.");
  case LongECModelType::Joint:
    return("Modelling weight for each MB-group pooled across all volumes.");
  case LongECModelType::IndividualTimeConstant:
    return("Modelling a separate time constant for each volume.");
  case LongECModelType::JointTimeConstant:
    return("Modelling a single time constant across all volumes.");
  default:
    throw EddyException("EddyCommandLineOptions::LongECModelString: Unknown model type.");
  }
} EddyCatch

enum class OffsetModelType { Linear, LinearPlusLag, Quadratic, QuadraticPlusLinearLag, QuadraticPlusLag, Unknown };

EDDY::OffsetModelType EddyCommandLineOptions::OffsetModel() const EddyTry
{
  if (IsDiffusion()) {
    if (_diff._offset_model.value() == string("linear")) return(EDDY::OffsetModelType::Linear);
    else if (_diff._offset_model.value() == string("linear_lag")) return(EDDY::OffsetModelType::LinearPlusLag);
    else if (_diff._offset_model.value() == string("quadratic")) return(EDDY::OffsetModelType::Quadratic);
    else if (_diff._offset_model.value() == string("linear_lag_quadratic")) return(EDDY::OffsetModelType::QuadraticPlusLinearLag);
    else if (_diff._offset_model.value() == string("linear_lag_quadratic_lag")) return(EDDY::OffsetModelType::QuadraticPlusLag);
    return(EDDY::OffsetModelType::Unknown);
  }
  return(EDDY::OffsetModelType::Unknown);
} EddyCatch

EDDY::OLSumStats EddyCommandLineOptions::OLSummaryStatistics() const EddyTry 
{ if (IsDiffusion()) { 
    if (_diff._ol_ss.value()==std::string("sw")) return(OLSumStats::ShellWise);
    else if (_ol_type.value()==std::string("pooled")) return(OLSumStats::Pooled);
    else return(OLSumStats::Unknown);
  }
  else return(OLSumStats::Unknown); 
} EddyCatch

NEWMAT::ColumnVector EddyCommandLineOptions::HyperParValues() const EddyTry
{
  NEWMAT::ColumnVector rval(_hypar_internal.size());
  for (unsigned int i=0; i<_hypar_internal.size(); i++) rval(i+1) = _hypar_internal[i];
  return(rval);
} EddyCatch

void EddyCommandLineOptions::SetHyperParValues(const NEWMAT::ColumnVector& p) EddyTry
{
  _hypar_internal.resize(p.Nrows());
  for (unsigned int i=0; i<_hypar_internal.size(); i++) _hypar_internal[i] = p(i+1);
  this->SetHyperParFixed(true);
  return;
} EddyCatch

std::vector<float> EddyCommandLineOptions::TestRotAngles() const EddyTry
{
  std::vector<float> rot = _diff._test_rot.value();
  if (rot.size() < 3) rot.resize(3,0.0);
  return(rot);
} EddyCatch

std::pair<int,double> EddyCommandLineOptions::ShellShapeReference(unsigned int i) const EddyTry
{
  if (IsDiffusion()) {
    if (i >= _diff._shell_shape_references.size()) throw EddyException("EddyCommandLineOptions::ShellShapeReference: Index out of bounds");
    return(std::pair<int,double>(_diff._shell_shape_references[i]._indx,_diff._shell_shape_references[i]._b_val));
  }
  return(std::pair<int,double>(-1,-999.0)); // Not difusion
} EddyCatch

NEWIMAGE::interpolation EddyCommandLineOptions::InterpolationMethod() const EddyTry
{
  if (_interp.value() == string("spline")) return(NEWIMAGE::spline);
  else if (_interp.value() == string("trilinear")) return(NEWIMAGE::trilinear);
  else throw EddyException("EddyCommandLineOptions::InterpolationMethod: Invalid interpolation option.");
} EddyCatch

NEWIMAGE::extrapolation EddyCommandLineOptions::ExtrapolationMethod() const EddyTry
{
  if (_extrap.value() == string("periodic")) return(NEWIMAGE::periodic);
  else if (_extrap.value() == string("mirror")) return(NEWIMAGE::mirror);
  else throw EddyException("EddyCommandLineOptions::ExtrapolationMethod: Invalid extrapolation option.");
} EddyCatch

NEWIMAGE::interpolation EddyCommandLineOptions::S2VInterpolationMethod() const EddyTry
{
  if (_s2v_interp.value() == string("spline")) return(NEWIMAGE::spline);
  else if (_s2v_interp.value() == string("trilinear")) return(NEWIMAGE::trilinear);
  else throw EddyException("EddyCommandLineOptions::S2VInterpolationMethod: Invalid interpolation option.");
} EddyCatch

DebugIndexClass::DebugIndexClass(const std::string& in) EddyTry
{
  _indx = parse_commaseparated_numbers(in);
} EddyCatch

std::vector<unsigned int> DebugIndexClass::parse_commaseparated_numbers(const std::string& list) const EddyTry
{
  std::vector<std::string> str_list = parse_commaseparated_list(list);
  std::vector<unsigned int> number_list(str_list.size(),0);
  for (unsigned int i=0; i<str_list.size(); i++) {
    number_list[i] = atoi(str_list[i].c_str());
  }

  return(number_list);
} EddyCatch

std::vector<std::string> DebugIndexClass::parse_commaseparated_list(const std::string&  list) const EddyTry
{
  std::vector<std::string> ostr;

  size_t cur_pos = 0;
  size_t new_pos = 0;
  unsigned int n=0;
  while ((new_pos = list.find_first_of(',',cur_pos)) != string::npos) {
    ostr.resize(++n);
    ostr[n-1] = list.substr(cur_pos,new_pos-cur_pos);
    cur_pos = new_pos+1;
  }
  ostr.resize(++n);
  ostr[n-1] = list.substr(cur_pos,string::npos);

  return(ostr);
} EddyCatch

S2VParam::S2VParam(const std::vector<int>& order, const std::vector<float>& lambda, const std::vector<float>& fwhm,
		   const std::vector<int>& niter) EddyTry : _order(order), _lambda(lambda), _fwhm(fwhm), _niter(niter)
{
  if (_order.size() != _lambda.size()) throw EddyException("Size of --s2v_lambda must match size of --mporder");
  if (_order.size() != _niter.size()) throw EddyException("Size of --s2v_niter must match size of --mporder");
  if (_fwhm.size() != 1 && _fwhm.size() != this->total_niter()) throw EddyException("--s2v_fwhm value must be given once or once per iteration");
  if (this->total_niter() > 100) throw EddyException("You have asked for more than 100 slice-to-volume iterations. Seriously?");
  // Reasonable slice-to-vol FWHM
  for (unsigned int i=0; i<_fwhm.size(); i++) {
    if (_fwhm[i] < 0.0 || _fwhm[i] > 20.0) throw EddyException("--s2v_fwhm value outside valid range 0-20mm");
  }
} EddyCatch

void EddyCommandLineOptions::do_initial_parsing(int argc, char *argv[]) EddyTry
{
  Utilities::OptionParser options(_title,_examples);

  try {
    // First mandatory parameters
    options.add(_imain);
    options.add(_mask);
    options.add(_acqp);
    options.add(_index);
    options.add(_out);
    if (IsDiffusion()) {
      options.add(_diff._bvecs);
      options.add(_diff._bvals);
      options.add(_diff._bdeltas);         // Not mandatory, but groups naturally with bvals and bvecs
      options.add(_diff._data_is_shelled); // Not mandatory, but groups naturally with bvals and bvecs
    }
    if (IsfMRI()) {
      options.add(_fmri._session);
      options.add(_fmri._rep_time);
      options.add(_fmri._echo_time);       // Not mandatory, but groups naturally with repetition time
    }
    // Parameters related to input susceptibility field
    options.add(_topup);
    options.add(_field);
    options.add(_field_mat);
    // Paramaters related to the general running of vanilla eddy
    options.add(_niter_tmp);
    options.add(_fwhm_tmp);
    options.add(_ref_scan_no);
    options.add(_interp);
    options.add(_extrap);
    options.add(_epvalid);
    if (IsDiffusion()) {
      options.add(_diff._resamp);
      options.add(_diff._lsr_lambda);
      options.add(_diff._brange);
    }
    // Parameters related to eddy current modelling
    if (IsDiffusion()) {
      options.add(_diff._flm);
      options.add(_diff._slm_str);
      options.add(_diff._b0_flm);
      options.add(_diff._b0_slm_str);
      options.add(_diff._rep_time);
      options.add(_diff._long_ec_str);
      options.add(_diff._long_ec_dont_reest);
      options.add(_diff._long_ec_sep_offs_move);	
    }
    // Parameters related to outlier detection/replacement
    options.add(_rep_ol);
    if (IsDiffusion()) options.add(_diff._rep_noise);
    options.add(_ol_nstd);
    options.add(_ol_nvox);
    options.add(_ol_ec);
    options.add(_ol_type);
    options.add(_ol_jacut);
    if (IsDiffusion()) {
      options.add(_diff._ol_ss);
      options.add(_diff._ol_pos);
      options.add(_diff._ol_sqr);
    }
    // Parameters related to slice-to-volume motion correction
    options.add(_mb);
    options.add(_slspec);
    options.add(_json);
    options.add(_mp_order);
    options.add(_s2v_lambda);
    options.add(_s2v_niter);
    options.add(_s2v_fwhm);
    options.add(_s2v_interp);
    options.add(_shape_ref_scan_nos);
    // Parameters related to susceptibility-by-movement interaction estimation
    options.add(_estimate_mbs);
    options.add(_mbs_niter);
    options.add(_mbs_lambda);
    options.add(_mbs_ksp);
    if (IsfMRI()) options.add(_fmri._fast_map);
    // Parameters related to the Gaussian process
    if (IsDiffusion()) options.add(_diff._covfunc);
    else if (IsfMRI()) options.add(_fmri._covfunc);
    options.add(_hyparcostfunc);
    options.add(_nvoxhp);
    options.add(_initrand);
    options.add(_hyparfudgefactor);
    options.add(_hypar);
    // Miscallenous parameters related to missing planes/output masking/bvec-rotation
    options.add(_fep);
    options.add(_dont_mask_output);
    options.add(_diff._rbvde);
    // Parameters related to initialisation of movement/distortion parameters
    options.add(_init);
    options.add(_init_s2v);
    options.add(_init_mbs);
    // Parameters related to "tricky" movements for diffusion
    if (IsDiffusion()) {
      options.add(_diff._dont_sep_offs_move);
      options.add(_diff._offset_model);
      options.add(_diff._dont_peas);
      options.add(_diff._use_b0s_for_peas);
      options.add(_diff._session);
    }
    // Parameters related to the writing of various "extra" output
    if (IsDiffusion()) {
      options.add(_diff._fields);
      options.add(_diff._write_cnr_maps);
      options.add(_diff._write_range_cnr_maps);
      options.add(_diff._write_scatter_brain_predictions);
    }
    if (IsfMRI()) options.add(_fmri._write_snr_maps);
    options.add(_dfields);
    options.add(_residuals);
    options.add(_with_outliers);
    options.add(_history);
    options.add(_write_slice_stats);
    options.add(_write_predictions);
    options.add(_no_text_files);
    // Parameters related to things that can be useful for debugging/development
    if (IsDiffusion()) {
      options.add(_diff._dwi_only);
      options.add(_diff._b0_only);
      options.add(_diff._test_rot);
      options.add(_diff._print_mi_values);
      options.add(_diff._print_mi_planes);
    }
    options.add(_debug_tmp);
    options.add(_dbg_indx_str);
    options.add(_log_timings);
    options.add(_nthr);
    // Finally some generally helpful parameters
    options.add(_verbose);
    options.add(_very_verbose);
    options.add(_help);
    if (IsDiffusion()) { // Here we add all the deprecated parameters from the "old" diffusion eddy
      options.add(_diff._sep_offs_move);
      options.add(_diff._peas);
      options.add(_diff._rms);
      options.add(_diff._mb_offs);
      options.add(_diff._slorder);
    }

    // Copies of argc and argv
    char **cargv = argv;
    int cargc = argc;
    // Remove fmri/diffusion tags if they are there
    if (argc > 1) {
      std::string tag(argv[1]);
      std::for_each(tag.begin(),tag.end(),[](char& c){c = std::tolower(c);});
      if (tag==std::string("fmri") || tag==std::string("diffusion")) {
	cargc = argc-1;
	cargv = new char*[cargc];
	cargv[0] = new char[strlen(argv[0])+1];
	std::strcpy(cargv[0],argv[0]);
	for (int i=1; i<cargc; i++) {
	  cargv[i] = new char[strlen(argv[i+1])+1];
	  std::strcpy(cargv[i],argv[i+1]);
	}
      }
    }
    // Parse command line
    // cout << "cargv[0] = " << cargv[0] << ", cargv[1] = " << cargv[1] << ", cargv[2] = " << cargv[2] << endl;
    int i=options.parse_command_line(cargc, cargv);
    if (i < cargc) {
      for (; i<cargc; i++) {
        cerr << "Unknown input: " << cargv[i] << endl;
      }
      exit(EXIT_FAILURE);
    }
    // Clean up if we modified argv
    if (cargv != argv) {
      for (int i=0; i<cargc; i++) delete[] cargv[i];
      delete[] cargv;
    }
    if (_help.value() || !options.check_compulsory_arguments(true)) {
      options.usage();
      exit(EXIT_FAILURE);
    }
  }
  catch(const Utilities::X_OptionError& e) {
    options.usage();
    cerr << endl << e.what() << endl;
    exit(EXIT_FAILURE);
  }

  // Write text files with
  // a. The command that the user used with any expansions performed by the shell
  // b. A file with values for all the parameters.
  ofstream cf;  // Command file
  cf.exceptions(ofstream::failbit | ofstream::badbit);
  try {
    cf.open(_out.value()+std::string(".eddy_command_txt"));
    cf << argv[0];
    for (int ii=1; ii<argc; ii++) cf << " " << argv[ii];
    cf << endl;
    cf.close();
  }
  catch(const ofstream::failure& e) {
    cout << "Writing file " << _out.value()+std::string(".eddy_command_txt") << " failed with message " << e.what() << endl;
    throw EddyException("EddyCommandLineOptions::do_initial_parsing: Writing file failed. Likely misspecification of --out");
  }
  ofstream dpf; // Detailed Parameter File
  dpf.exceptions(ofstream::failbit | ofstream::badbit);
  try {
    Utilities::detailed_output(dpf);
    dpf.open(_out.value()+std::string(".eddy_values_of_all_input_parameters"));
    dpf << "This file contains values (including defaults) for all input parameters to eddy." << endl;
    dpf << options << endl;
    dpf.close();
  }
  catch(const ofstream::failure& e) {
    cout << "Writing file " << _out.value()+std::string(".eddy_values_of_all_input_parameters") << " failed with message " << e.what() << endl;
    throw EddyException("EddyCommandLineOptions::do_initial_parsing: Writing file failed. Likely misspecification of --out");
  }
} EddyCatch

// Makes sure that the indicies makes sense w.r.t. pointing in to
// the acquistition parameter file, and stores them in _indvec.
// N.B. It is _not_ a mistake that the matrix is passed by value.
bool EddyCommandLineOptions::indicies_kosher(NEWMAT::Matrix indx, NEWMAT::Matrix acqp) EddyTry
{
  if (indx.Ncols() > indx.Nrows()) indx = indx.t();
  _indvec.resize(indx.Nrows());
  unsigned int max_indx = static_cast<unsigned int>(indx(1,1));
  for (int i=0; i<indx.Nrows(); i++) {
    _indvec[i] = static_cast<unsigned int>(indx(i+1,1));
    if (fabs(static_cast<double>(_indvec[i])-indx(i+1,1)) > 1e-6) return(false);
    max_indx = std::max(max_indx,_indvec[i]);
  }
  if (max_indx > static_cast<unsigned int>(acqp.Nrows())) return(false);
  else return(true);
} EddyCatch

// Makes sure that the session indicies makes sense and stores them 
// in _sessvec. N.B. It is _not_ a mistake that the matrix is passed 
// by value. 
bool EddyCommandLineOptions::session_indicies_kosher(NEWMAT::Matrix indx) EddyTry
{
  if (indx.Ncols() > indx.Nrows()) indx = indx.t();
  _sessvec.resize(indx.Nrows());
  std::vector<unsigned int> bins;
  for (int i=0; i<indx.Nrows(); i++) {
    unsigned int ii = static_cast<unsigned int>(indx(i+1,1));
    _sessvec[i] = ii;
    if (bins.size() < ii) bins.resize(ii,0);
    bins[ii-1]++;
  }
  // Checks that all sessions have a "reasonable" number of volumes and
  // that they are not too unevenly distributed.
  unsigned int maxbin = 0; unsigned int minbin = 1e6;
  double meanbin = 0.0;
  for (int i=0; i<bins.size(); i++) {
    maxbin = std::max(maxbin,bins[i]);
    minbin = std::min(minbin,bins[i]);
    meanbin += static_cast<double>(bins[i]);
  }
  meanbin /= static_cast<double>(bins.size());
  if (minbin==0 || maxbin > 3*minbin || maxbin > 6*meanbin) return(false);
  return(true);
} EddyCatch

/****************************************************************//**
*
* \brief Routine for sorting shape-reference indicies
*
* Sorts the indices in 'shape_refs' using the shell indicies
* in 'shells' as a key.
*
********************************************************************/
void EddyCommandLineOptions::set_shape_refs(const std::vector<int>&   shape_refs,
					    const std::vector<int>&   shells,
					    const std::vector<double> shell_b_vals) EddyTry
{
  struct shape_ref_shell_pair {
    shape_ref_shell_pair(int ri, int si) : _ri(ri), _si(si) {}
    int _ri; // Index (into vector of all scans) of refs
    int _si; // Index of shells
  };
  std::vector<shape_ref_shell_pair> ri_si_v;
  for (unsigned int i=0; i<shape_refs.size(); i++) ri_si_v.push_back(shape_ref_shell_pair(shape_refs[i],shells[i]));
  std::sort(ri_si_v.begin(),ri_si_v.end(),[](shape_ref_shell_pair& p1,shape_ref_shell_pair& p2){ return(p1._si < p2._si); });
  unsigned int offset = 0;
  if (EddyUtils::Isb0(EDDY::DiffPara(shell_b_vals[0]))) {
    _diff._b0_shape_reference = ri_si_v[0]._ri;
    offset = 1;
  }
  _diff._shell_shape_references.clear();
  for (unsigned int i=offset; i<ri_si_v.size(); i++) {
    _diff._shell_shape_references.push_back(diff::shell_shape_ref_struct(ri_si_v[i]._ri,shell_b_vals[i]));
  }
  _shape_references_set = true;
} EddyCatch

/****************************************************************//**
*
* \brief Routine for reading slice-orders-spec file
*
* The file should be a text-file with as many entries as there are
* slice-groups in the data. In the simplest case that is simply the
* # of slices, and for example in the case of multi-band 3 it would
* be #of_slices/3.
* If the content of the file is for example
* 1 3 5 ... 2 4 6 ...
* It means that the first group (that containing the most basal slice)
* was acquired first followed by the third group (that containing the
* third most basal slice) etc.
*
********************************************************************/
std::vector<unsigned int> EddyCommandLineOptions::get_slorder(const std::string& fname,
							      unsigned int       ngrp) const EddyTry
{
  std::vector<unsigned int> rval;
  try {
    NEWMAT::Matrix tmp = MISCMATHS::read_ascii_matrix(fname);
    if (tmp.Ncols() > tmp.Nrows()) tmp = tmp.t();
    int n = tmp.Nrows();
    int one = tmp.Ncols();
    rval.resize(n);
    if (n != int(ngrp) || one != 1) throw EddyException("Size mismatch between --imain and --slorder file");
    for (int i=0; i<n; i++) {
      if (tmp(i+1,1)<0 || tmp(i+1,1)>ngrp-1) throw EddyException("--slorder file contains invalid entry");
      rval[i] = tmp(i+1,1);
    }
  }
  catch(...) { throw EddyException("Error when attempting to read --slorder file " + fname); }

  return(rval);
} EddyCatch