diff --git a/.gitattributes b/.gitattributes
new file mode 100644
index 0000000000000000000000000000000000000000..07f04632fc22dee40c5f68bde9406197619cb128
--- /dev/null
+++ b/.gitattributes
@@ -0,0 +1 @@
+*.nii.gz filter=lfs diff=lfs merge=lfs -text
diff --git a/fslstats.cc b/fslstats.cc
index 83a987301fa76c02be5faef947a71d6acf180485..8e90fd22080a2e5a0917f1cccbf5eba447e2b46e 100644
--- a/fslstats.cc
+++ b/fslstats.cc
@@ -8,21 +8,27 @@
 
 #include <limits>
 #include <string>
+#include <typeinfo>
 
 #include "utils/fsl_isfinite.h"
 #include "armawrap/newmat.h"
 #include "miscmaths/miscmaths.h"
 #include "newimage/newimageall.h"
 #include "newimage/costfns.h"
+#include "nlohmann/json.hpp"
 
 using namespace std;
 using namespace NEWMAT;
 using namespace NEWIMAGE;
+using json = nlohmann::json;
 
 int print_usage(const string &progname)
 {
   cout << "Usage: fslstats [preoptions] <input> [options]" << endl
        << endl;
+  // add --json output format option
+  // keys will correspond to the options above
+  cout << "preoption -json: output in JSON format to standard out" << endl;
   cout << "preoption -t will give a separate output line for each 3D volume of a 4D timeseries" << endl;
   cout << "preoption -K < indexMask > will generate seperate n submasks from indexMask, for indexvalues 1..n where n is the maximum index value in indexMask, and generate statistics for each submask" << endl;
   cout << "Note - options are applied in order, e.g. -M -l 10 -M will report the non-zero mean, apply a threshold and then report the new nonzero mean" << endl
@@ -51,11 +57,7 @@ int print_usage(const string &progname)
   cout << "-k <mask>    : use the specified image (filename) for masking - overrides lower and upper thresholds" << endl;
   cout << "-d <image>   : take the difference between the base image and the image specified here" << endl;
   cout << "-h <nbins>   : output a histogram (for the thresholded/masked voxels only) with nbins" << endl;
-  cout << "-H <nbins> <min> <max>   : output a histogram (for the thresholded/masked voxels only) with nbins and histogram limits of min and max" << endl;
-  // add --json output format option
-  // keys will correspond to the options above
-  cout << "--json       : output in JSON format" << endl
-       << endl;
+  cout << "-H <nbins> <min> <max>   : output a histogram (for the thresholded/masked voxels only) with nbins and histogram limits of min and max" << endl; 
   cout << "Note - thresholds are not inclusive ie lthresh<allowed<uthresh" << endl;
   return 1;
 }
@@ -120,38 +122,35 @@ int generate_masks(volume4D<float> &mask, volume4D<float> &masknz, const volume4
   return generateNonZeroMask(mask, masknz, input);
 }
 
-// print function that can handle outputting in JSON format or not
-void print_value(const string &key, const string  &value, string &json, const bool jsonOutput)
+// print function that can handle outputting in JSON format or not. Use a template to allow for different types
+template <typename T>
+void update_json(string key, T &value, json &j)
 {
-  // append the key and value to the json string
-  if (jsonOutput)
-  {
-    json.append("\"" + key + "\": \"" + value + "\", ");
-  }
-  else
-  {
-    cout << value << " ";
-  }
+  j[key] = value;
 }
 
+template <typename T>
+void update_json(string key, initializer_list<T> &value, json &j)
+{
+  j[key] = value;
+}
 
-int fmrib_main_float(int argc, char *argv[], const bool timeseriesMode, const string &indexMaskName)
+template <typename T>
+void update_json(string key, vector<T> &value, json &j)
 {
-  // determine if the user wants to output in JSON format
-  bool jsonOutput = false;
-  for (int i = 1; i < argc; i++)
-  {
-    if (strcmp(argv[i], "--json") == 0)
-    {
-      jsonOutput = true;
-      break;
-    }
-  }
-  // if JSON output is requested, then we need to store the output as a string
-  // and then output it at the end via cout
-  // don't forget to close the json at the end with a "}"
-  string json = "";
-  
+  j[key] = value;
+}
+
+void update_json(string key, ColumnVector &value, json &j)
+{
+  j[key] = value;
+}
+
+int fmrib_main_float(int argc, char *argv[], const bool timeseriesMode, const string &indexMaskName, const bool jsonMode = false)
+{
+
+  json js;
+
   cout.setf(ios::dec);
   cout.setf(ios::fixed, ios::floatfield);
   cout.setf(ios::left, ios::adjustfield);
@@ -167,27 +166,16 @@ int fmrib_main_float(int argc, char *argv[], const bool timeseriesMode, const st
     read_volume4D(indexMask, indexMaskName);
   int nTimepoints((timeseriesMode) ? inputMaster.tsize() : 1), nIndices((indexMaskName != "") ? indexMask.max() : 1);
 
+  // initialise the volumes array in the json object. This will get appended to later
+  js["volumes"] = json::array();
+
   for (int timepoint = 0; timepoint < nTimepoints; timepoint++)
   {
-    if (jsonOutput)
-    {
-      if (timepoint == 0){
-        json.append("{\"volumes\": [");
-      } else {
-        json.append(", ");
-      }
-    }
+    json tp = json::object();
+    tp["indices"] = json::array();
     for (int index = 1; index <= nIndices; index++)
     {
-      // if json mode is enabled
-      if (jsonOutput)
-      {
-        if (index == 1){
-          json.append("{\"indices\": [{");
-        } else {
-          json.append(", {");
-        }
-      }
+      json idx = json::object();
       if (timeseriesMode)
         vol = inputMaster[timepoint];
       volume<float> mask, masknz;
@@ -229,23 +217,41 @@ int fmrib_main_float(int argc, char *argv[], const bool timeseriesMode, const st
         else if (sarg == "-m")
         {
           if (mask.nvoxels() > 0){
-            string val = to_string(vol.mean(mask));
-            print_value(sarg, val, json, jsonOutput);
+            auto val = vol.mean(mask);
+            if (jsonMode){
+              update_json(sarg, val, idx);
+            } else {
+              cout << val << " ";
+            }
           } else {
-            string val = to_string(vol.mean());
-            print_value(sarg, val, json, jsonOutput);
+            // print_value(sarg, val, json, jsonMode);
+            auto val = vol.mean();
+            if (jsonMode){
+              update_json(sarg, val, idx);
+            } else {
+              cout << val << " ";
+            }
           }
         }
         else if (sarg == "-M")
         {
-          if (masknz.nvoxels() > 0)
-            print_value(sarg, to_string(vol.mean(masknz)), json, jsonOutput);
+          if (masknz.nvoxels() > 0){
+            auto val = vol.mean(masknz);
+            if (jsonMode){
+              update_json(sarg, val, idx);
+            } else {
+              cout << val << " ";
+            }
+          }
           else
           {
             double nzmean = 0;
             nzmean = nonzeromean(vol);
-            string val = to_string(nzmean);
-            print_value(sarg, val, json, jsonOutput);
+            if (jsonMode){
+              update_json(sarg, nzmean, idx);
+            } else {
+              cout << nzmean << " ";
+            }
           }
         }
         else if (sarg == "-X")
@@ -259,17 +265,12 @@ int fmrib_main_float(int argc, char *argv[], const bool timeseriesMode, const st
           auto coord1 = MISCMATHS::round(coord(1));
           auto coord2 = MISCMATHS::round(coord(2));
           auto coord3 = MISCMATHS::round(coord(3));
-          string coordString = "";
-          coordString.append(to_string(coord1));
-          coordString.append(" ");
-          coordString.append(to_string(coord2));
-          coordString.append(" ");
-          coordString.append(to_string(coord3));
-          print_value(
-            sarg,
-            coordString,
-            json,
-            jsonOutput);
+          auto coordList = {coord1, coord2, coord3};
+          if (jsonMode){
+            update_json(sarg, coordList, idx);
+          } else {
+            cout << coord1 << " " << coord2 << " " << coord3 << " ";
+          }
           
         }
         else if (sarg == "-x")
@@ -283,17 +284,12 @@ int fmrib_main_float(int argc, char *argv[], const bool timeseriesMode, const st
           auto coord1 = MISCMATHS::round(coord(1));
           auto coord2 = MISCMATHS::round(coord(2));
           auto coord3 = MISCMATHS::round(coord(3));
-          string coordString = "";
-          coordString.append(to_string(coord1));
-          coordString.append(" ");
-          coordString.append(to_string(coord2));
-          coordString.append(" ");
-          coordString.append(to_string(coord3));
-          print_value(
-            sarg,
-            coordString,
-            json,
-            jsonOutput);
+          auto coordList = {coord1, coord2, coord3};
+          if (jsonMode){
+            update_json(sarg, coordList, idx);
+          } else {
+            cout << coord1 << " " << coord2 << " " << coord3 << " ";
+          }
         }
         else if (sarg == "-w")
         {
@@ -357,28 +353,12 @@ int fmrib_main_float(int argc, char *argv[], const bool timeseriesMode, const st
             zmax = zmin;
             zmin = tmp;
           }
-          // now output nifti coords
-          string coords = "";
-          coords.append(to_string(xmin));
-          coords.append(" ");
-          coords.append(to_string(1 + xmax - xmin));
-          coords.append(" ");
-          coords.append(to_string(ymin));
-          coords.append(" ");
-          coords.append(to_string(1 + ymax - ymin));
-          coords.append(" ");
-          coords.append(to_string(zmin));
-          coords.append(" ");
-          coords.append(to_string(1 + zmax - zmin));
-          coords.append(" ");
-          coords.append(to_string(tmin));
-          coords.append(" ");
-          coords.append(to_string(1 + tmax - tmin));
-          print_value(
-            sarg,
-            coords,
-            json,
-            jsonOutput);
+          vector<int> coords = {xmin, 1 + xmax - xmin, ymin, 1 + ymax - ymin, zmin, 1 + zmax - zmin, tmin, 1 + tmax - tmin};
+          if (jsonMode){
+            update_json(sarg, coords, idx);
+          } else {
+            cout << coords[0] << " " << coords[1] << " " << coords[2] << " " << coords[3] << " " << coords[4] << " " << coords[5] << " " << coords[6] << " " << coords[7] << " ";
+          }
         }
         else if (sarg == "-e")
         {
@@ -400,12 +380,11 @@ int fmrib_main_float(int argc, char *argv[], const bool timeseriesMode, const st
             }
           }
           entropy /= log((double)nbins);
-          print_value(
-            sarg,
-            to_string(entropy),
-            json,
-            jsonOutput
-          );
+          if (jsonMode){
+            update_json(sarg, entropy, idx);
+          } else {
+            cout << entropy << " ";
+          }
         }
         else if (sarg == "-E")
         {
@@ -427,12 +406,11 @@ int fmrib_main_float(int argc, char *argv[], const bool timeseriesMode, const st
             }
           }
           entropy /= log((double)nbins);
-          print_value(
-            sarg,
-            to_string(entropy),
-            json,
-            jsonOutput
-          );
+          if (jsonMode){
+            update_json(sarg, entropy, idx);
+          } else {
+            cout << entropy << " ";
+          }
         }
         else if (sarg == "-k")
         {
@@ -547,125 +525,125 @@ int fmrib_main_float(int argc, char *argv[], const bool timeseriesMode, const st
             if (mask.tsize() == 1)
               nvox = nvox * vol.tsize();
             
-            string nvox_str = "";
-            nvox_str.append(to_string(nvox));
-            nvox_str.append(" ");
-            nvox_str.append(to_string(nvox * vol.xdim() * vol.ydim() * vol.zdim()));
-            print_value(
-              sarg,
-              nvox_str,
-              json,
-              jsonOutput
-            );
+            auto volume_v = nvox * vol.xdim() * vol.ydim() * vol.zdim();
+            // promote nvox to float to maintain types in list
+            auto nvox_volume = {static_cast<float>(nvox), volume_v};
+            if (jsonMode){
+              update_json(sarg, nvox_volume, idx);
+            } else {
+              cout << nvox << " " << volume_v << " ";
+            }
           }
           else
           {
-            string nvox_str = "";
-            nvox_str.append(to_string((long int)vol.nvoxels() * vol.tsize()));
-            nvox_str.append(" ");
-            nvox_str.append(to_string(vol.nvoxels() * vol.tsize() * vol.xdim() * vol.ydim() * vol.zdim()));
-            print_value(
-              sarg,
-              nvox_str,
-              json,
-              jsonOutput);
+            long int nvox = vol.nvoxels() * vol.tsize();
+            auto volume_v = nvox * vol.xdim() * vol.ydim() * vol.zdim();
+            // promote nvox to float to maintain types in list
+            auto nvox_volume = {static_cast<float>(nvox), volume_v};
+            if (jsonMode){
+              update_json(sarg, nvox_volume, idx);
+            } else {
+              cout << nvox << " " << volume_v << " ";
+            }
           }
         }
         else if (sarg == "-V")
         {
           if (masknz.nvoxels() > 0)
           {
-            string nvox_str = "";
-            nvox_str.append(to_string((long int)masknz.sum()));
-            nvox_str.append(" ");
-            nvox_str.append(to_string(masknz.sum() * vol.xdim() * vol.ydim() * vol.zdim()));
-            print_value(
-              sarg,
-              nvox_str,
-              json,
-              jsonOutput);
+            long int nvox = masknz.sum();
+            auto volume_V = nvox * vol.xdim() * vol.ydim() * vol.zdim();
+            // promote nvox to float to maintain types in list
+            auto nvox_volume = {static_cast<float>(nvox), volume_V};
+            if (jsonMode){
+              update_json(sarg, nvox_volume, idx);
+            } else {
+              cout << nvox << " " << volume_V << " ";
+            }
           }
           else
           {
             long int nzvox;
             nzvox = nonzerocount(vol);
-            string nzvox_str = "";
-            nzvox_str.append(to_string(nzvox));
-            nzvox_str.append(" ");
-            nzvox_str.append(to_string(nzvox * vol.xdim() * vol.ydim() * vol.zdim()));
-            print_value(
-              sarg,
-              nzvox_str,
-              json,
-              jsonOutput);
+            auto volume_V = nzvox * vol.xdim() * vol.ydim() * vol.zdim();
+            // promate nzvox to float to maintain types in list
+            auto nvox_volume = {static_cast<float>(nzvox), volume_V};
+            if (jsonMode){
+              update_json(sarg, nvox_volume, idx);
+            } else {
+              cout << nzvox << " " << volume_V << " ";
+            }
           }
         }
         else if (sarg == "-D")
         {
           // hidden debug option!
-          print_value(
-            sarg,
-            to_string(vol.sum()),
-            json,
-            jsonOutput
-          );
+          auto sum = vol.sum();
+          if (jsonMode){
+            update_json(sarg, sum, idx);
+          } else {
+            cout << sum << " ";
+          }
         }
         else if (sarg == "-s")
         {
-          if (mask.nvoxels() > 0)
-            print_value(
-              sarg,
-              to_string(vol.stddev(mask)),
-              json,
-              jsonOutput
-            );
-          else
-            print_value(
-              sarg,
-              to_string(vol.stddev()),
-              json,
-              jsonOutput
-            );
+          if (mask.nvoxels() > 0){
+            auto stddev = vol.stddev(mask);
+            if (jsonMode){
+              update_json(sarg, stddev, idx);
+            } else {
+              cout << stddev << " ";
+            }
+          }
+          else {
+            auto stddev = vol.stddev();
+            if (jsonMode){
+              update_json(sarg, stddev, idx);
+            } else {
+              cout << stddev << " ";
+            }
+          }
         }
         else if (sarg == "-S")
         {
           if (masknz.nvoxels() > 0)
           {
-            print_value(
-              sarg,
-              to_string(vol.stddev(masknz)),
-              json,
-              jsonOutput
-            );
+            auto stddev = vol.stddev(masknz);
+            if (jsonMode){
+              update_json(sarg, stddev, idx);
+            } else {
+              cout << stddev << " ";
+            }
           }
           else
           {
-            print_value(
-              sarg,
-              to_string(nonzerostddev(vol)),
-              json,
-              jsonOutput
-            );
+            auto stddev = nonzerostddev(vol);
+            if (jsonMode){
+              update_json(sarg, stddev, idx);
+            } else {
+              cout << stddev << " ";
+            }
           }
         }
         else if (sarg == "-r")
         {
           vector<float> limits(vol.robustlimits(mask));
-          print_value(
-            sarg,
-            to_string(limits[0]).append(" ").append(to_string(limits[1])),
-            json,
-            jsonOutput
-          );
+          if (jsonMode){
+            update_json(sarg, limits, idx);
+          } else {
+            cout << limits[0] << " " << limits[1] << " ";
+          }
         }
         else if (sarg == "-R")
         {
-          print_value(
-            sarg,
-            to_string(vol.min(mask)).append(" ").append(to_string(vol.max(mask))),
-            json,
-            jsonOutput
-          );
+          auto min_R = vol.min(mask);
+          auto max_R = vol.max(mask);
+          auto minmax_R = {min_R, max_R};
+          if (jsonMode){
+            update_json(sarg, minmax_R, idx);
+          } else {
+            cout << min_R << " " << max_R << " ";
+          }
         }
         else if (sarg == "-c")
         {
@@ -673,12 +651,15 @@ int fmrib_main_float(int argc, char *argv[], const bool timeseriesMode, const st
           cog.SubMatrix(1, 3, 1, 1) = vol.cog();
           cog(4) = 1.0;
           cog = vol.newimagevox2mm_mat() * cog;
-          print_value(
-            sarg,
-            to_string(cog(1)).append(" ").append(to_string(cog(2))).append(" ").append(to_string(cog(3))),
-            json,
-            jsonOutput
-          );
+          // remove the last element
+          cog = cog.Rows(1, 3);
+          // reshape to a row vector
+          cog = cog.t();
+          if (jsonMode){
+            update_json(sarg, cog, idx);
+          } else {
+            cout << cog(1) << " " << cog(2) << " " << cog(3) << " ";
+          }
         }
         else if (sarg == "-C")
         {
@@ -687,12 +668,15 @@ int fmrib_main_float(int argc, char *argv[], const bool timeseriesMode, const st
           cog.SubMatrix(1, 3, 1, 1) = vol.cog();
           cog(4) = 1.0;
           cog = vol.niftivox2newimagevox_mat().i() * cog;
-          print_value(
-            sarg,
-            to_string(cog(1)).append(" ").append(to_string(cog(2))).append(" ").append(to_string(cog(3))),
-            json,
-            jsonOutput
-          );
+          // remove the last element
+          cog = cog.Rows(1, 3);
+          // reshape to a row vector
+          cog = cog.t();
+          if (jsonMode){
+            update_json(sarg, cog, idx);
+          } else {
+            cout << cog(1) << " " << cog(2) << " " << cog(3) << " ";
+          }
         }
         else if (sarg == "-p")
         {
@@ -712,20 +696,22 @@ int fmrib_main_float(int argc, char *argv[], const bool timeseriesMode, const st
             cerr << "Percentile must be between 0 and 100" << endl;
             exit(1);
           }
-          if (mask.nvoxels() > 0)
-            print_value(
-              sarg,
-              to_string(vol.percentile((float)n / 100.0, mask)),
-              json,
-              jsonOutput
-            );
-          else
-            print_value(
-              sarg,
-              to_string(vol.percentile((float)n / 100.0)),
-              json,
-              jsonOutput
-            );
+          if (mask.nvoxels() > 0){
+            auto percentile_p = vol.percentile((float)n / 100.0, mask);
+            if (jsonMode){
+              update_json(sarg, percentile_p, idx);
+            } else {
+              cout << percentile_p << " ";
+            }
+          }
+          else {
+            auto percentile_p = vol.percentile((float)n / 100.0);
+            if (jsonMode){
+              update_json(sarg, percentile_p, idx);
+            } else {
+              cout << percentile_p << " ";
+            }
+          }
         }
         else if (sarg == "-P")
         {
@@ -750,12 +736,12 @@ int fmrib_main_float(int argc, char *argv[], const bool timeseriesMode, const st
             generate_masks(mask, masknz, vol, lthr, uthr);
             vol *= mask;
           }
-          print_value(
-            sarg,
-            to_string(vol.percentile((float)n / 100.0, masknz)),
-            json,
-            jsonOutput
-          );
+          auto percentile_P = vol.percentile((float)n / 100.0, masknz);
+          if (jsonMode){
+            update_json(sarg, percentile_P, idx);
+          } else {
+            cout << percentile_P << " ";
+          }
         }
         else if (sarg == "-h")
         {
@@ -779,41 +765,20 @@ int fmrib_main_float(int argc, char *argv[], const bool timeseriesMode, const st
           if (mask.nvoxels() > 0)
           {
             auto hist = vol.histogram(nbins, vol.min(), vol.max(), mask);
-            // loop over histogram values and and build a string to print
-            string hist_str = "";
-            for (int i = 0; i < hist.Nrows(); i++)
-            {
-              hist_str += to_string(hist(i+1));
-              if (i < hist.Nrows() - 1)
-              {
-                hist_str += " ";
-              }
+            if (jsonMode){
+              update_json(sarg, hist, idx);
+            } else {
+              cout << hist << " ";
             }
-            print_value(
-              sarg,
-              hist_str,
-              json,
-              jsonOutput
-            );
           }
           else
           {
             auto hist = vol.histogram(nbins, vol.min(), vol.max());
-            string hist_str = "";
-            for (int i = 0; i < hist.Nrows(); i++)
-            {
-              hist_str += to_string(hist(i+1));
-              if (i < hist.Nrows() - 1)
-              {
-                hist_str += " ";
-              }
+            if (jsonMode){
+              update_json(sarg, hist, idx);
+            } else {
+              cout << hist << " ";
             }
-            print_value(
-              sarg,
-              hist_str,
-              json,
-              jsonOutput
-            );
           }
         }
         else if (sarg == "-H")
@@ -859,91 +824,41 @@ int fmrib_main_float(int argc, char *argv[], const bool timeseriesMode, const st
           }
           if (mask.nvoxels() > 0)
           {
-            auto hist = vol.histogram(nbins, min, max, mask);
-            string hist_str = "";
-            for (int i = 0; i < hist.Nrows(); i++)
-            {
-              hist_str += to_string(hist(i+1));
-              if (i < hist.Nrows() - 1)
-              {
-                hist_str += " ";
-              }
+            auto hist = vol.histogram(nbins,min,max,mask);
+            if (jsonMode){
+              update_json(sarg, hist, idx);
+            } else {
+              cout << hist << " ";
             }
-            print_value(
-              sarg,
-              hist_str,
-              json,
-              jsonOutput
-            );
           }
           else
           {
             auto hist = vol.histogram(nbins, min, max);
-            string hist_str = "";
-            for (int i = 0; i < hist.Nrows(); i++)
-            {
-              hist_str += to_string(hist(i+1));
-              if (i < hist.Nrows() - 1)
-              {
-                hist_str += " ";
-              }
+            if (jsonMode){
+              update_json(sarg, hist, idx);
+            } else {
+              cout << hist << " ";
             }
-            print_value(
-              sarg,
-              hist_str,
-              json,
-              jsonOutput
-            );
           }
         }
         else
         {
-          // only print cerr if not the --json option
-          if (sarg != "--json")
-          {
-            cerr << "Unrecognised option: " << sarg << endl;
-            exit(3);
-          } else if (sarg == "--json" && last) {
-              // remove the trailing comma and space
-              json.pop_back();
-              json.pop_back();
-          }
+          cerr << "Unrecognised option: " << sarg << endl;
+          exit(3);
         }
-
         narg++;
       }
-      // if json mode is enabled
-      if (jsonOutput)
-      {
-        // close the json object
-        json.append("}");
-      }
-      if (index < nIndices && !jsonOutput)
-      {
+      tp["indices"] += idx;
+      if (!jsonMode)
         cout << endl;
-      }
-      // cout << endl;
     } // end nIndices
-    // if json mode is enabled
-    if (jsonOutput)
-    {
-      // close the json object
-      json.append("]}");
-    }
-    if (timepoint < nTimepoints - 1 && !jsonOutput)
-    {
-      cout << endl;
-    }
-    // cout << endl;
+    // append the timepoint to the the volumes array
+    js["volumes"] += tp;
   } // end timepoints
-  if (jsonOutput)
-  {
-    // close the json object
-    json.append("]}");
-    cout << json;
-  } else {
-    cout << endl;
-  }
+  if (jsonMode){
+    // print the json with 4 spaces of indentation
+    cout << js.dump(4) << endl;
+  }   
   return 0;
 }
 
@@ -953,11 +868,14 @@ int main(int argc, char *argv[])
   Tracer tr("main");
   string progname(argv[0]);
   bool timeseriesMode(false);
+  bool jsonMode(false);
   string indexMask("");
-  while (argc > 2 && (string(argv[1]) == "-t" || string(argv[1]) == "-K"))
+  while (argc > 2 && (string(argv[1]) == "-t" || string(argv[1]) == "-K" || string(argv[1]) == "-json"))
   {
     if (string(argv[1]) == "-t")
       timeseriesMode = true;
+    if (string(argv[1]) == "-json")
+      jsonMode = true;
     if (string(argv[1]) == "-K")
     {
       indexMask = string(argv[2]);
@@ -971,7 +889,7 @@ int main(int argc, char *argv[])
     return print_usage(progname);
   try
   {
-    return fmrib_main_float(argc, argv, timeseriesMode, indexMask);
+    return fmrib_main_float(argc, argv, timeseriesMode, indexMask, jsonMode);
   }
   catch (std::exception &e)
   {
diff --git a/tests/fslstats/feedsRun b/tests/fslstats/feedsRun
index c6d130521f54644eec25efea7243e2b4c08204b4..a007f1dd74c02c6dcb61231246aa462a4af03f61 100755
--- a/tests/fslstats/feedsRun
+++ b/tests/fslstats/feedsRun
@@ -53,6 +53,8 @@ options = [
     {"option": "-v", "expected": "1000 1000.000000"},
     {"option": "-V", "expected": "1000 1000.000000"},
     {"option": "-m", "expected": "0.495922"},
+    {"option": "-m -a", "expected": "0.495922"},
+    {"option": "-m -a -n", "expected": "0.495922"},
     {"option": "-M", "expected": "0.495922"},
     {"option": "-s", "expected": "0.290744"},
     {"option": "-S", "expected": "0.290744"},
@@ -63,6 +65,12 @@ options = [
     {"option": "-C", "expected": "4.498706 4.545873 4.613188"},
     {"option": "-p 50", "expected": "0.482584"},
     {"option": "-P 50", "expected": "0.482584"},
+    {"option": "-m -h 10 -c", "expected": "0.495922 101.000000 \n99.000000 \n108.000000 \n107.000000 \n102.000000 \n93.000000 \n99.000000 \n88.000000 \n95.000000 \n108.000000 \n 4.498706 4.545873 4.613188"},
+    {"option": "-m -h 10", "expected": "0.495922 101.000000 \n99.000000 \n108.000000 \n107.000000 \n102.000000 \n93.000000 \n99.000000 \n88.000000 \n95.000000 \n108.000000"},
+    {"option": "-h 10", "expected": "101.000000 \n99.000000 \n108.000000 \n107.000000 \n102.000000 \n93.000000 \n99.000000 \n88.000000 \n95.000000 \n108.000000"},
+    {"option": "-H 10 0 1", "expected": "98.000000 \n102.000000 \n107.000000 \n108.000000 \n102.000000 \n93.000000 \n99.000000 \n88.000000 \n95.000000 \n108.000000"},
+    {"option": "-l 0.25 -u 0.75 -m", "expected": "0.495150"},
+    {"option": "-d imageForDiff.nii.gz -m", "expected": "0.000000"},
 ]
 
 # create a list of tests and expected results
@@ -82,55 +90,24 @@ tests_with_preoptions = [
         "preoptions": "-K",
         "mask": "mask.nii.gz",
         "options": "-m",
-        "expected": "0.948930 \n0.947671 \n1.003258 \n1.010696",
+        "expected": "1.078044 \n1.028827 \n0.986428 \n1.020869",
     },
     {
         "preoptions": "-t -K",
         "mask": "mask.nii.gz",
         "options": "-m",
-        "expected": "0.459736 \n0.476035 \n0.504080 \n0.549485 \n0.489194 \n0.471636 \n0.499178 \n0.461211",
+        "expected": "0.526682 \n0.515652 \n0.492337 \n0.511661 \n0.551362 \n0.513176 \n0.494091 \n0.509208",
     },
     {
         "preoptions": "-t",
         "mask": "",
         "options": "-m",
-        "expected": "0.496675 \n0.487950",
+        "expected": "0.503236 \n0.504035",
     },
 ]
 
-# taken from fslchpixdim test
-def create_image(shape, pixdim):
-    pixdim = list(pixdim)
-    data   = np.random.random(shape).astype(np.float32)
-    hdr    = nib.Nifti1Header()
-    hdr.set_data_dtype(np.float32)
-    hdr.set_data_shape(shape)
-    hdr.set_zooms(pixdim[:len(shape)])
-
-    return nib.Nifti1Image(data, np.eye(4), hdr)
-
-def create_mask(input_img, n_rois=4):
-    '''
-    Create a mask image from the input image.
-    The mask image will have n_rois different values with random sizes randomly placed in the image.
-    '''
-    data = input_img.get_fdata()
-    mask = np.zeros_like(data)
-    for i in range(n_rois + 1):
-        roi_size = np.random.randint(1, data.size)
-        roi_value = i
-        roi = np.random.choice(data.size, roi_size, replace=False)
-        mask.flat[roi] = roi_value
-    return nib.Nifti1Image(mask, input_img.affine, input_img.header)
-
-
 def test_fslstats():
     imgfile = 'test.nii.gz'
-    shape   = (10,10,10)
-    img     = create_image(shape, [1] * len(shape))
-    img.to_filename(imgfile)
-    mask = create_mask(img)
-    mask.to_filename('mask.nii.gz')
 
     # run the tests witoout preoptions
     for test in tests:
@@ -151,9 +128,6 @@ def test_fslstats():
 
     # run the tests with preoptions
     imgfile = 'test_4d.nii.gz'
-    shape   = (10,10,10, 2)
-    img     = create_image(shape, [1] * len(shape))
-    img.to_filename(imgfile)
     for test in tests_with_preoptions:
         cmd = f"fslstats {test['preoptions']} {test['mask']} {imgfile} {test['options']}"
         # remove any double spaces from empty test options
@@ -175,6 +149,4 @@ def test_fslstats():
 
 
 if __name__ == "__main__":
-    if len(sys.argv) > 1:
-        os.chdir(sys.argv[1])
     test_fslstats()
diff --git a/tests/fslstats/imageForDiff.nii.gz b/tests/fslstats/imageForDiff.nii.gz
new file mode 100644
index 0000000000000000000000000000000000000000..f4d00fe3675aa6b9dd763c59d6da23ec402c2b17
--- /dev/null
+++ b/tests/fslstats/imageForDiff.nii.gz
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:14105cd775bf8e192299695a90848e285690c2dafb086bdd8830c5665fb0394b
+size 3716
diff --git a/tests/fslstats/mask.nii.gz b/tests/fslstats/mask.nii.gz
new file mode 100644
index 0000000000000000000000000000000000000000..3a7c1553ca342e57efe8b8ffd49bd880b6832c7f
--- /dev/null
+++ b/tests/fslstats/mask.nii.gz
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:d9a0b1719ca95cf36b44ddcde54121fa28411d060ab0d288303c80935eb429b1
+size 786
diff --git a/tests/fslstats/test.nii.gz b/tests/fslstats/test.nii.gz
new file mode 100644
index 0000000000000000000000000000000000000000..f4d00fe3675aa6b9dd763c59d6da23ec402c2b17
--- /dev/null
+++ b/tests/fslstats/test.nii.gz
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:14105cd775bf8e192299695a90848e285690c2dafb086bdd8830c5665fb0394b
+size 3716
diff --git a/tests/fslstats/test_4d.nii.gz b/tests/fslstats/test_4d.nii.gz
new file mode 100644
index 0000000000000000000000000000000000000000..247023d54a81d5c7994b0561ed2fc17768cbf514
--- /dev/null
+++ b/tests/fslstats/test_4d.nii.gz
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:a48fc73f54d9aba889968a99d58339d495950301970f166166bede44dd38e623
+size 7304