diff --git a/doc/building.rst b/doc/building.rst index 6e60a27..de97700 100644 --- a/doc/building.rst +++ b/doc/building.rst @@ -5,16 +5,17 @@ I have an official FSL distribution installed --------------------------------------------- You can build Fabber using the FSL build system. First you need to set -up your development environment: - -:: +up your development environment:: source $FSLDIR/etc/fslconf/fsl-devel.sh export FSLDEVDIR= -Then building Fabber should be a case of: +``FSLDEVDIR`` is an anternate prefix to ``FSLDIR`` which is used to +store updated code separately from the official FSL release. Most +FSL-based scripts should use code installed in ``FSLDEVDIR`` in preference +to the main FSL release code. -:: +Building Fabber should be a case of:: cd fabber_core make install @@ -23,75 +24,59 @@ This approach uses the same build tools as the rest of FSL which is important on some platforms, notably OSX. It will install the updated code into whatever prefix you selected as ``FSLDEVDIR``. -I don't have an FSL distribution which supports development ------------------------------------------------------------ - -Fabber has an alternative build system using ``cmake`` as its build -tool. CMake is cross-platform and is designed for out-of-source builds, -so you create a separate build directory and all the compiled files end -up there. This helps to keep your source tree free of build artifacts. - -Some convenience scripts are provided to build and install Fabber. You -need to ensure that ``FSLDIR`` is set to point to wherever the FSL -dependencies are installed, then for example to do a build which -includes debug symbols on Linux or OSX you run: - -:: +Building new or updated model libraries +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - scripts/build.sh debug +Model libraries are distributed separately from the Fabber core. +If you need an updated version of a model library, for example +the ASL model library, you first need to get the source code +for the models library. A number of model libraries are +available in our `Github repositories `_ +all named ``fabber_models_``. -The executables and libraries will end up in ``build_debug``. Other -options are: +Then to build and install the updated model libraries you would then +run, for example:: -:: + cd fabber_models_asl + make install - scripts/build.sh release - - scripts/build.sh relwithdebinfo - -The last example is quite useful - it produces a release build with full -optimization (typically about 2x faster than a debug build) but keeps -the debug symbols in the executable so debugging is possible. +Adding your own models +---------------------- -You can also look at the scripts - they are very simple - to see the -actual CMake commands being run. +If you want to create your own model to use with the Fabber core +model fitting engine, see `Building a new model`_. Once you've +designed and coded your model there are two ways to incorporate +it into the Fabber system: -Pitfalls -~~~~~~~~ +Adding models directly to the core +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -gcc vs Clang on OSX -^^^^^^^^^^^^^^^^^^^ +If you wish, you can add your own models directly into the Fabber source +tree and build the executable as above. This is not generally +recommended because your model will be built into the core executable, however +it can be the quickest way to get an existing model built in. You will +need to follow these steps: -OSX has Apple’s ``Clang`` compiler and the LLVM C++ standard library -``libc++`` as its default for compiling C++. However FSL uses ``gcc`` -and the GNU ``libstdc++`` library for OSX builds. These are not -compatible. If you are building against a standard FSL installation you -need to force CMake to use ``gcc`` and ``libstdc++``. To do this: +1. Add your model source code into the fabber_core directory, e.g.:: -1. Uncomment the following line in CMakeLists.txt + fabber_core/fwdmodel_mine.cc + fabber_core/fwdmodel_mine.h - set(CMAKE_CXX_FLAGS “${CMAKE_CXX_FLAGS} -std=c++11 - -stdlib=libstdc++”) +2. Edit ``Makefile`` to add your model to the list of core objects, e.g.:: -2. Add the following directives to the ``cmake`` command itself (edit - the scripts/build.sh script if you are using it): + COREOBJS = fwdmodel_mine.o noisemodel.o fwdmodel.o inference.o fwdmodel_linear.o fwdmodel_poly.o convergence.o motioncorr.o priors.o transforms.o - -DCMAKE_C_COMPILER=/usr/bin/gcc -DCMAKE_CXX_COMPILER=/usr/bin/g++ +3. Run ``make install`` again to build and install a new executable -Adding your own models -~~~~~~~~~~~~~~~~~~~~~~ +Creating a new models library +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -If you wish, you can add your own models directly into the Fabber source -tree and build the executable as above. This is not generally -recommended because your model will be built into libfabbercore, however -it can be the quickest way to get an existing model built in. You will -need to follow these steps: +This is the preferred approach if you want to distribute your new models. A template +for a new model library including a simple sine-function implementation is +included with the Fabber source code in ``fabber_core/examples``. See +`Building a new model`_ for a full tutorial on this example which includes +how to set up the build scripts. -1. Add your model source code into the fabber_core directory, - e.g. \ ``fabber_core/fwdmodel_mine.cc`` and - ``fabber_core/fwdmodel_mine.h`` +.. _Building a new model: models -2. Edit CMakeLists.txt to add your model sources as follows:: - # Core objects - things that implement the framework for inference - set(CORE_SRC noisemodel.cc fwdmodel.cc inference.cc) diff --git a/doc/exp_test_biexp_good.png b/doc/exp_test_biexp_good.png new file mode 100644 index 0000000..c074bcf Binary files /dev/null and b/doc/exp_test_biexp_good.png differ diff --git a/doc/exp_test_biexp_improved.png b/doc/exp_test_biexp_improved.png new file mode 100644 index 0000000..0b3e7cd Binary files /dev/null and b/doc/exp_test_biexp_improved.png differ diff --git a/doc/exp_test_single.png b/doc/exp_test_single.png new file mode 100644 index 0000000..5f976a0 Binary files /dev/null and b/doc/exp_test_single.png differ diff --git a/doc/getting.rst b/doc/getting.rst index 6ee5e3a..daf16b2 100644 --- a/doc/getting.rst +++ b/doc/getting.rst @@ -1,16 +1,37 @@ Getting Fabber ============== +From FSL +-------- + Fabber is distributed as part of `FSL `_, and this is the easiest way to get Fabber if you want to use existing models for fMRI data . This documentation describes the -version of Fabber included with FSL v6.0.1 and above. +version of Fabber included with *FSL v6.0.1 and above*. + +Addition tools that can use Fabber will work with a correctly installed +FSL distribution although currently not all models we have developed are +available in FSL. + +Standalone Fabber distribution +------------------------------ + +Standalone versions of Fabber including a selection of model libraries are available +for a number of platforms. These may be useful if you don't want the rest of FSL +or if you need a more up to date version of Fabber than the one included with FSL. + +The current standalone release can be found at `https://github.com/ibme-qubic/fabber_core/releases`_. + +This distribution can be used with tools requiring a Fabber installation such as +the `Python API `_, or Fabber-based plugins +for `Quantiphyse `_. You should set the environment +variable ``FABBERDIR`` to the unpacked distribution directory to ensure these tools +can find Fabber. -If you want to build your own models, or if you need a more up to date -version of Fabber than the one included with FSL, you can build Fabber from -the source code available in the `Github repository `_. -Pre-built binaries may be available for some platforms. +Building from source code +------------------------- -To build the source code, see `Building Fabber`_. +You can build Fabber from the source code available in the `Github repository `_. +You will need an FSL installation for this. For instructions see `Building Fabber`_. .. _Building Fabber: building diff --git a/doc/index.rst b/doc/index.rst index 14ae4dd..d61f6b5 100644 --- a/doc/index.rst +++ b/doc/index.rst @@ -46,8 +46,8 @@ for ASL, CEST, DSC, DCE and dual echo fMRI data. Fabber is distributed as part of `FSL `_, however you should ensure that you are using **FSL v6.0.1 as a minimum version**. -A GUI tool for processing ASL, CEST and (soon) DSC and DCE data using Fabber -is `Quantiphyse `_. +The `Quantiphyse `_ visual analysis +tool contains plugins to analyse ASL, CEST, DSC and DCE fMRI data using Fabber. .. toctree:: :maxdepth: 1 diff --git a/doc/models.rst b/doc/models.rst index c0c89dc..99c29c0 100644 --- a/doc/models.rst +++ b/doc/models.rst @@ -1,35 +1,34 @@ -Building new models -=================== +Building a new model +==================== -For most applications, a model will need to be constructed. This will +For most new applications, a model will need to be constructed. This will include adjustable parameters which Fabber will then fit. A complete example model is provided in the ``examples`` subdirectory of the Fabber source code. This provides an easy template to implement a new model. In the next section we will go through this example. +We will assume only some basic knowledge of ``C++`` for this example. + A simple example ---------------- To create a new Fabber model it is necessary to create an instance of the class ``FwdModel``. As an example, we will create a model which fits the -data to a sine function with an amplitude, frequency and phase shift, -plus an optional constant offset: +data to sum of exponential functions, each with an amplitude and decay rate. .. math:: - A\sin(B(t+C)) [+D] + \sum{A_n\exp(-R_nt) .. note:: - The source code and ``CMakeLists.txt`` file for this example are in the - fabber source code, in the examples subdirectory. We will assume you + The source code and ``Makefile`` file for this example are in the + Fabber source code, in the ``examples`` subdirectory. We will assume you have this to hand as we go through the process! -First we will create the interface .h file which shows the methods we -will need to implement: - -.. code:: +First we will create the interface ``fwdmodel_exp.h`` file which shows the methods we +will need to implement:: - // fwdmodel_sine.h - A simple sine curve fitting model + // fwdmodel_exp.h - A simple exponential sum model #pragma once #include "fabber_core/fwdmodel.h" @@ -39,12 +38,12 @@ will need to implement: #include #include - class SineFwdModel : public FwdModel { + class ExpFwdModel : public FwdModel { public: static FwdModel* NewInstance(); - SineFwdModel() - : m_include_offset(false) + ExpFwdModel() + : m_num(1), m_dt(1.0) { } @@ -56,125 +55,171 @@ will need to implement: void EvaluateModel(const NEWMAT::ColumnVector ¶ms, NEWMAT::ColumnVector &result, const std::string &key="") const; - + + protected: + void GetParameterDefaults(std::vector ¶ms) const; + private: - bool m_include_offset; - static FactoryRegistration registration; + int m_num; + double m_dt; + static FactoryRegistration registration; }; We have not made our methods virtual, so nobody will be able to create a -subclass of our model. If we wanted this to be the case all the +subclass of our model. If we wanted this to be the case all themake non-static methods would need to be virtual, and we would need to add a -virtual destructor. +virtual destructor. This is sometimes useful when you want to create +variations on a basic model. -We have also declared a private variable to hold whether we are going to -use the optional offset or not. +Most of the code above is completely generic to any model. The only parts which +are specific to our exp-function model are: -We will now implement these methods one by one. Many of them are -straightforward. We start our implementation file as follows + - The name ``ExpFwdModel`` + - The private variables ``m_num`` (the number of exponentials in our sum) and + ``m_dt`` (the time between data points). -:: +We will now implement these methods one by one. Many of them are +straightforward. We start our implementation file ``fwdmodel_exp.cc`` as follows:: - // fwdmodel_sine.cc - Implements a simple sine curve fitting model - #include "fwdmodel_sine.h" + // fwdmodel_exp.cc - Implements a simple exp curve fitting model + #include "fwdmodel_exp.h" - #include "fabber_core/fwdmodel.h" + #include #include using namespace std; using namespace NEWMAT; -Here are the first few methods: +This just declares some standard headers we will use. If you prefer to fully qualify your +namespaces you can leave out the ``using namespace`` lines. -:: +We need to implement a couple of methods to ensure that our model is visible to the +Fabber system:: - FactoryRegistration SineFwdModel::registration("sine"); + FactoryRegistration ExpFwdModel::registration("exp"); - FwdModel* SineFwdModel::NewInstance() + FwdModel* ExpFwdModel::NewInstance() { - return new SineFwdModel(); + return new ExpFwdModel(); } The first line here registers our model so that it is known to Fabber by -the name *sine*. The second line is a *Factory method* used so that +the name ``exp`` The second line is a *Factory method* used so that Fabber can create a new instance of our model when its name appears on -the command line. +the command line:: -:: - - string SineFwdModel::ModelVersion() const - { - return "1.0"; - } + string ExpFwdModel::ModelVersion() const + { + return "1.0"; + } - string SineFwdModel::GetDescription() const + string ExpFwdModel::GetDescription() const { - return "Example model which uses a sine function"; + return "Example model of a sum of exponentials"; } We’ve given our model a version number, if we update it at some later stage we should change the number returned so anybody using the model will know it has changed and what version they have. There's also -a brief description which fabber will return when the user requests help on the model. - -:: +a brief description which fabber will return when the user requests help on the model:: - static int NUM_OPTIONS = 1; - static OptionSpec OPTIONS[] = - { - {"use-offset", OPT_BOOL, "If True, allow an additional constant offset parameter", OPT_NONREQ, "false"}, - {""} - }; + static OptionSpec OPTIONS[] = { + { "dt", OPT_FLOAT, "Time separation between samples", OPT_REQ, "" }, + { "num-exps", OPT_INT, "Number of independent exponentials in sum", OPT_NONREQ, "1" }, + { "" } + }; - void SineFwdModel::GetOptions(vector &opts) const - { - for (int i = 0; OPTIONS[i].name != ""; i++) - { - opts.push_back(OPTIONS[i]); - } - } + void ExpFwdModel::GetOptions(vector &opts) const + { + for (int i = 0; OPTIONS[i].name != ""; i++) + { + opts.push_back(OPTIONS[i]); + } + } This is the suggested way to declare the options that your model can -take. It is a little cumbersome when there is only one option, but if -you have many options it will make it clear to see what they are. - -The OptionSpec definition contains, in order, the option name, a type -indicator (OPT_STR, OPT_INT, OPT_BOOL, OPT_FILE) which indicates what -kind of data is expected, a human readable description, whether the -option must be specified (OPT_REQ) or not (OPT_NONREQ), and finally the -default value if any. - -We have a single non-mandatory Boolean option - whether to allow the -extra constant offset. If not specified, it defaults to false, so not -including an offset. +take - in this case the user can choose how many exponentials to include in the sum +and what the time resolution in the data is. Each option is listed in the ``OPTIONS`` array which +**ends with an empty option** (important!). + +An option is described by: + + - It's name which generally should *not* include underscores (hyphen is OK as in this + case). This translates into a command line option ``--use-offset``. + - An option type. Possibilities are: + - ``OPT_BOOL`` for a Yes/No boolean + - ``OPT_FLOAT`` for a decimal number + - ``OPT_INT`` for a whole number + - ``OPT_STR`` for text + - ``OPT_MATRIX`` for a small matrix (specified by giving the filename of + a text file which contains the matrix data in tab-separated form) + - ``OPT_IMAGE`` for a 3D image specified as a Nifti file + - ``OPT_TIMESERIES`` for a 4D image specified as a Nifti file + - ``OPT_FILE`` for a generic filename + - A brief description of the option. This will be displayed when ``--help`` is + requested for the model + - ``OPT_NONREQ`` if the option is not mandatory (does not need to be specified) + or ``OPT_REQ`` if the option must be provided by the user. + - An indication of the default value. This value is not actually used to initialize + anything but is shown in ``--help`` to explain to the user what the default is + if the option is not given. So it can contain any text (e.g. ``"0.7 for PASL, 1.3 for pCASL"``. + You should not specify a default for a mandatory option (``OPT_REQ``) + +In this case we have made the time resolution option mandatory, but the number +of exponentials defaults to 1 if not specified. + +This option system is a little cumbersome when there is only one option, but if +you have many it will make it clear to see what they are. Most +real models will have many configuration options, for example an ASL +model will need to know details of the sequence such as the TIs/PLDs, +the bolus duration, the labelling method, number of repeats, etc... + +Options specified by the user are captured in the ``FabberRunData`` +object which we use to set the variables in our model class +in the ``Initialize`` method. ``Initialize`` is called before the model +will be used. Its purpose is to allow the model to set up any internal +variables based on the user-supplied options. Here we capture the +time resolution option and the number of exponentials - note that +the latter has a default value. + + void ExpFwdModel::Initialize(FabberRunData& rundata) + { + m_dt = rundata.GetDouble("dt"); + m_num = rundata.GetIntDefault("num-exps", 1); + } -:: +We use the term *Options* to distinguish user-specified or default model +configuration from *Parameters* which are the parts of the model inferred by +the Fabber process. Next we need to specify what parameters our model +includes:: - void SineFwdModel::GetParameterDefaults(std::vector ¶ms) const + void ExpFwdModel::GetParameterDefaults(std::vector ¶ms) const { params.clear(); int p=0; - params.push_back(Parameter(p++, "a", DistParams(1, 1e6), DistParams(1, 1e6))); - params.push_back(Parameter(p++, "b", DistParams(1, 1e6), DistParams(1, 1e6))); - params.push_back(Parameter(p++, "c", DistParams(0, 1e6), DistParams(0, 1e6))); - if (m_include_offset) { - params.push_back(Parameter(p++, "d", DistParams(0, 1e6), DistParams(0, 1e6))); + for (int i=0; i`` and ``r`` for each exponential +in the sum, where ```` is 1, 2, ... As well as a name, each parameter has two ``DistParams`` instances defining the *prior* and *initial posterior* distribution for the parameter. -The ``DistParams`` instances take two parameters - a mean and a variance. +``DistParams`` take two parameters - a mean and a variance. At this +point we will diverge slightly to explain what these mean. Priors and Posteriors ~~~~~~~~~~~~~~~~~~~~~ + *Priors* are central to Bayesian inference, and describe the extent of our belief about a parameter's -value before we have seen any data. +value *before we have seen any data*. For example if a parameter represents the T_1 value of grey matter in the brain there is a well known range of plausible values. By declaring a @@ -192,137 +237,333 @@ The second ``DistParams`` instance represents the initial *posterior*. This is t point for the optimisation as it tries to find the best values for each parameter. Usually this does not matter too much and can often be set to be identical to the prior. -Sometimes it -is helpful to set this to a more restrictive value (lower variance) to avoid numerical -instability. it is also possible to adjust this on a per-voxel basis - i.e. when the data -being fitted is available. We will not do that here, but it can be useful when fitting, for +Sometimes, however, it may be helpful to give the initial posterior a more restrictive (lower) +variance to avoid numerical instability. + +It is also possible to adjust the initial posterior on a per-voxel basis using the actual +voxel data. We will not do that here, but it can be useful when fitting, for example, a constant offset, where we can tell the optimisation to start with a value that is the mean of the data. This may help avoid instability and local minima. In general it is against the spirit of the Bayesian approach to modify the priors on the -basis of the data, and no means are provided to do this. It is possible to modify the priors -on a global basis but this is not encouraged and in general a model should try to provide +basis of the data, and no means are provided to do this. It is possible for the user to modify +the priors on a global basis but this is not encouraged and in general a model should try to provide good priors that will not need modification. -:: +We now go back to our model code where we finally reach the point where we write +the code to calculate our model:: - void SineFwdModel::Initialize(FabberRunData& rundata) - { - m_include_offset = rundata.GetBool("use-offset"); - } + void ExpFwdModel::EvaluateModel(const NEWMAT::ColumnVector ¶ms, + NEWMAT::ColumnVector &result, + const std::string &key) const + { + result.ReSize(data.Nrows()); + result = 0; + + for (int i=0; i/libfabber_models_sine.so + #!/bin/env python + import sys + import traceback -Both of these should contain the new model, as you can show by using the -``--listmodels`` option + from fabber import self_test, FabberException + + save = "--save" in sys.argv + try: + rundata= { + "model" : "exp", # Exponential model + "num-exps" : 1, # Single exponential function + "dt" : 0.02, # With 100 time points time values will range from 0 to 2 + } + params = { + "amp1" : [1, 0.5], # Amplitude + "r1" : [1.0, 0.8], # Decay rate + } + test_config = { + "nt" : 100, # Number of time points + "noise" : 0.1, # Amplitude of Gaussian noise to add to simulated data + "patchsize" : 20, # Each patch is 20 voxels along each dimension + } + result, log = self_test("exp", rundata, params, save_input=save, save_output=save, invert=True, **test_config) + except FabberException, e: + print e.log + traceback.print_exc() + except: + traceback.print_exc() + +The test script generates a test Nifti image containing 'patches' of +data chequerboard style, each of which corresponds to a combination +of true parameter values. As Fabber is designed to work on 3D timeseries +data you can only vary three model parameters in each test - others +must have fixed values. + +The test data is generated both 'clean' and with added Gaussian +noise of specified amplitude. The model is then run on the noisy +data to determine how closely the true parameter values can +be recovered. In this case we get the following output:: + + python test_single.py --save + + Running self test for model exp + Saving test data to Nifti file: test_data_exp + Saving clean data to Nifti file: test_data_exp_clean + Inverting test data - running Fabber: 100% + + Parameter: amp1 + Input 1.000000 -> 0.999701 Output + Input 0.500000 -> 0.500674 Output + Parameter: r1 + Input 1.000000 -> 1.000728 Output + Input 0.800000 -> 0.801230 Output + Noise: Input 0.100000 -> 0.099521 Output + +For each parameter, the input (`ground truth`) value is given and +also the mean inferred value across the patch. In this case +it has recovered the parameters pretty well on average. An +example plot of a single voxel might look like this: + +.. image:: exp_test_single.png + +The orange line is the noisy data it's trying to fit while the +two smooth lines represent the 'true' data and the model fit. +In fact for this example typically the model fit is much closer +to the true data - we have chosen this voxel as an example +so it is possible to see them separately! + +Testing the model - bi-exponential +---------------------------------- -Running our simple example --------------------------- +Fitting to a single exponential is not too challenging - here +we will test fitting to a bi-exponential where there are two +different decay rates. We will find that we need to improve +the model to get a better fit. -We will show how to run the model using the GUI tool as it makes the -results more visually obvious. However you can run from the command line -if you prefer. +First we can modify the test script to test a bi-exponential +(``test_biexp.py`` in examples):: -On creating a new Fabber configuration file we need to ensure that we -are using the correct Fabber executable with our model built in to it. -Once this is set, you will see that the ‘sine’ model appears in the list -of forward models. Clicking on ‘Options’ then shows our single option. -We will leave it unset initially. + #!/bin/env python -.. figure:: /uploads/e50fd8068e61fa5d778e582df8be0ec3/sine_options.PNG - :alt: sine_options + import sys + import traceback - sine_options + from fabber import self_test, FabberException -We now click on ‘Run’ using the nlls inference method and a sample fMRI -data set. After some time, Fabber completes and we can see a list of our -output files. If we make the model fit visible and click on points in -the image data we can see how well our model fits the data. + save = "--save" in sys.argv + try: + rundata= { + "model" : "exp", + "num-exps" : 2, + "dt" : 0.02, + "max-iterations" : 50, + } + params = { + "amp1" : [1, 0.5], # Amplitude first exponential + "amp2" : 0.5, # Amplitude second exponential + "r1" : [1.0, 0.8], # Decay rate of first exponential + "r2" : 6.0, # Decay rate of second exponential + } + test_config = { + "nt" : 100, # Number of time points + "noise" : 0.1, # Amplitude of Gaussian noise to add to simulated data + "patchsize" : 20, # Each patch is 20 voxels along each dimension + } + result, log = self_test("exp", rundata, params, save_input=save, save_output=save, invert=True, **test_config) + except FabberException, e: + print e.log + traceback.print_exc() + except: + traceback.print_exc() + +This is similar to the last test but we have set ``num-exps`` to 2 and added +parameters for a fixed second exponential curve with a faster decay rate. +If we run this we get output something like this:: + + python test_biexp.py --save + Running self test for model exp + Saving test data to Nifti file: test_data_exp + Saving clean data to Nifti file: test_data_exp_clean + Inverting test data - running Fabber: 100% + + Parameter: amp1 + Input 1.000000 -> 0.633822 Output + Input 0.500000 -> 0.309912 Output + Parameter: r1 + Input 1.000000 -> 19693700210770313216.000000 Output + Input 0.800000 -> -324689116576874496.000000 Output + Noise: Input 0.100000 -> 0.150277 Output + +This isn't looking too encouraging. If we examine the model fit +against the data we find that actually most voxels have fitted +quite well: + +.. image:: exp_test_biexp_good.png + +However a few voxels have ended up with very unrealistic +parameter values. This kind of behaviour is a risk with model fitting - +in trying to find the best solution the inference can end up +finding a local minimum which is a long way from the true +minimum. + +We will show two additions we can make to our model to improve this +behaviour. + +Initialising the posterior +~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The initial posterior is a 'first guess' at the parameter values +and can be based on the data. Fabber models can use their knowledge +of the model to make a better guess by overriding the ``InitVoxelPosterior`` +method. We firstly add this method to ``fwdmodel_exp.h`` + + void InitVoxelPosterior(MVNDist &posterior) const; + +Now we implement it in ``fwdmodel_exp.cc`` + + void ExpFwdModel::InitVoxelPosterior(MVNDist &posterior) const + { + double data_max = data.Maximum(); + + for (int i=0; i ¶ms) const + { + params.clear(); + + int p=0; + for (int i=0; i`` at the top of ``fwdmodel_exp.cc``. - modelfit_sine_nlls_no_offset +With these changes we still retain some bad fitting voxels but +fewer than previously. The output of the test script is now:: -Not very well - because we allowed no offset, the best fit is nearly -linear. Let’s switch that option on, and try again. + python test_biexp.py --save + Running self test for model exp + Saving test data to Nifti file: test_data_exp + Saving clean data to Nifti file: test_data_exp_clean + Inverting test data - running Fabber: 100% -.. figure:: /uploads/03e57bbd242d9d8de1ef3c47627d4dac/modelfit_sine_nlls_offset.PNG - :alt: modelfit_sine_nlls_offset + Parameter: amp1 + Input 1.000000 -> 9.651313 Output + Input 0.500000 -> 0.499837 Output + Parameter: r1 + Input 1.000000 -> 6.453277 Output + Input 0.800000 -> 124170.593750 Output + Noise: Input 0.100000 -> 0.099531 Output - modelfit_sine_nlls_offset +So we have a reduction in the number of extreme values. In this case we can't actually trust the +self-test output because sometimes the inference 'swaps' the exponentials around making +``amp1`` = ``amp2`` and ``r1`` = ``r2``. But viewing the model fit +visually shows sensible fitting in the overwhelming majority of voxels:: -That’s a bit better - it’s able to pick up the general variation in the -signal. Although it’s clear that this example is not giving us any -physical information we can see how the model has tried to reproduce the -general shape of the data by varying the allowed parameters. +.. image:: exp_test_biexp_improved.png Changing the example to your own model -------------------------------------- -To implement a single new model, it should be as simple as: +To summaries, these are the main steps you'll need to take to +change this example into your own new model: -- Edit the source files, ``Makefile`` and ``CMakeLists.txt`` to change - references to ``sine`` to the name of your model -- Rename source files, e.g. \ ``fwdmodel_sine.cc`` -> - ``fwdmodel_.cc`` +- Edit the ``Makefile`` to change references to ``exp`` to the name of your model +- Rename source files, e.g. ``fwdmodel_exp.cc`` -> ``fwdmodel_.cc`` - Add your model options to the options list in the ``.cc`` file -- Implement the ``Evaluate`` and ``HardcodedInitialDists`` methods for - your model +- Add any model-specific private variables in the ``.h`` file +- Implement the ``Initialize``, ``GetParameterDefaults``, ``Evaluate`` methods for + your model. +- If required, implement ``InitVoxelPosterior``