diff --git a/src/basis.cpp b/src/basis.cpp
index 1b8d1aeb..b7a674dc 100644
--- a/src/basis.cpp
+++ b/src/basis.cpp
@@ -4,14 +4,14 @@
 
 namespace rl {
 
-Basis::Basis(
-  Eigen::ArrayXXf const &dyn,
+void Basis(
+  Eigen::ArrayXXf const &dynamics,
   float const thresh,
   Index const nBasis,
   bool const demean,
   bool const rotate,
-  bool const normalize)
-  : dynamics{dyn}
+  bool const normalize,
+  HD5::Writer &writer)
 {
   // Calculate SVD - observations are in cols
   Eigen::ArrayXXf d = normalize ? dynamics.colwise().normalized() : dynamics;
@@ -31,7 +31,7 @@ Basis::Basis(
   }
   Log::Print("Retaining {} basis vectors, cumulative energy: {}", nRetain, cumsum.head(nRetain).transpose());
 
-  basis = svd.V.leftCols(nRetain);  
+  Eigen::MatrixXf basis = svd.V.leftCols(nRetain);
   if (rotate) {
     Log::Print("Bastardised Gram-Schmidt");
     basis.col(0) = basis.rowwise().mean().normalized();
@@ -45,13 +45,10 @@ Basis::Basis(
 
   Log::Print("Computing dictionary");
   basis *= std::sqrt(basis.rows());
-  dict = basis.transpose() * dynamics.matrix();
-  norm = dict.colwise().norm();
+  Eigen::MatrixXf dict = basis.transpose() * dynamics.matrix();
+  Eigen::ArrayXf norm = dict.colwise().norm();
   dict = dict.array().rowwise() / norm.transpose();
-}
 
-void Basis::write(HD5::Writer &writer)
-{
   writer.writeMatrix(Eigen::MatrixXf(basis.transpose()), HD5::Keys::Basis);
   writer.writeMatrix(dict, HD5::Keys::Dictionary);
   writer.writeMatrix(norm, HD5::Keys::Norm);
diff --git a/src/basis.hpp b/src/basis.hpp
index 59e2c9df..b8ff9405 100644
--- a/src/basis.hpp
+++ b/src/basis.hpp
@@ -6,20 +6,13 @@
 
 namespace rl {
 
-struct Basis
-{
-  Basis(
-    Eigen::ArrayXXf const &dynamics,
-    float const thresh,
-    Index const nB,
-    bool const demean,
-    bool const rotate,
-    bool const normalize);
-  void write(HD5::Writer &writer);
-
-  Eigen::ArrayXXf dynamics;
-  Eigen::MatrixXf basis, dict;
-  Eigen::ArrayXf norm;
-};
+void Basis(
+  Eigen::ArrayXXf const &dynamics,
+  float const thresh,
+  Index const nB,
+  bool const demean,
+  bool const rotate,
+  bool const normalize,
+  HD5::Writer &writer);
 
 } // namespace rl
diff --git a/src/cmd/basis-img.cpp b/src/cmd/basis-img.cpp
index 8f2a1704..795f26d6 100644
--- a/src/cmd/basis-img.cpp
+++ b/src/cmd/basis-img.cpp
@@ -14,6 +14,8 @@ int main_basis_img(args::Subparser &parser)
   args::Positional<std::string> iname(parser, "INPUT", "Input image file");
   args::Positional<std::string> oname(parser, "OUTPUT", "Name for the basis file");
 
+  args::ValueFlag<Index> ss(parser, "S", "Subsample image", {"subsamp"}, 1);
+  args::ValueFlag<Index> spf(parser, "S", "Spokes per frame", {"spf"}, 1);
   args::ValueFlag<Index> nBasis(parser, "N", "Number of basis vectors to retain (overrides threshold)", {"nbasis"}, 5);
   args::Flag demean(parser, "C", "Mean-center dynamics", {"demean"});
   args::Flag rotate(parser, "V", "Rotate basis", {"rotate"});
@@ -24,7 +26,10 @@ int main_basis_img(args::Subparser &parser)
   }
 
   HD5::Reader reader(iname.Get());
-  Cx4 const img = reader.readSlab<Cx4>(HD5::Keys::Image, 0);
+  Cx4 img = reader.readSlab<Cx4>(HD5::Keys::Image, 0);
+  if (ss) {
+    img = Cx4(img.stride(Sz4{1, ss.Get(), ss.Get(), ss.Get()}));
+  }
   Sz4 const shape = img.dimensions();
   Cx4 const ref = img.slice(Sz4{shape[0] - 1, 0, 0, 0}, AddFront(LastN<3>(shape), 1));
   Re4 const real = (img / ref.broadcast(Sz4{shape[0], 1, 1, 1})).real();
@@ -32,16 +37,17 @@ int main_basis_img(args::Subparser &parser)
   Re3 const toMask = ref.chip<0>(0).abs();
   auto const toMaskMat = CollapseToArray(toMask);
   auto const [thresh, count] = Otsu(toMaskMat);
-  Eigen::MatrixXf dynamics(shape[0], count);
+  Index const fullN = shape[0] * spf.Get();
+  Eigen::MatrixXf dynamics(fullN, count);
   Index col = 0;
+  auto inds = Eigen::ArrayXi::LinSpaced(fullN, 0, shape[0] - 1);
   for (Index ii = 0; ii < realMat.cols(); ii++) {
     if (toMaskMat(ii) >= thresh) {
-      dynamics.col(col) = realMat.col(ii);
+      dynamics.col(col) = realMat.col(ii)(inds).transpose();
       col++;
     }
   }
-  Basis basis(dynamics, 99.f, nBasis.Get(), demean, rotate, normalize);
   HD5::Writer writer(oname.Get());
-  basis.write(writer);
+  Basis(dynamics, 99.f, nBasis.Get(), demean, rotate, normalize, writer);
   return EXIT_SUCCESS;
 }
diff --git a/src/cmd/basis-sim.cpp b/src/cmd/basis-sim.cpp
index 962fe5a5..da787eb0 100644
--- a/src/cmd/basis-sim.cpp
+++ b/src/cmd/basis-sim.cpp
@@ -31,7 +31,7 @@ auto Run(rl::Settings const &s, Index const nsamp, std::vector<std::vector<float
   parameters.setZero();
   Index totalP = 0;
   for (size_t ii = 0; ii < los.size(); ii++) {
-    Log::Print("Parameter set {}/{}. Low {} High {}", ii+1, los.size(), fmt::join(los[ii], "/"), fmt::join(his[ii], "/"));
+    Log::Print("Parameter set {}/{}. Low {} High {}", ii + 1, los.size(), fmt::join(los[ii], "/"), fmt::join(his[ii], "/"));
     auto p = simulator.parameters(nsamp, los[ii], his[ii]);
     parameters.middleCols(totalP, p.cols()) = p;
     totalP += p.cols();
@@ -135,9 +135,8 @@ int main_basis_sim(args::Subparser &parser)
   case Sequences::DWI: std::tie(pars, dyns) = Run<rl::DWI>(settings, nsamp.Get(), pLo.Get(), pHi.Get()); break;
   }
 
-  Basis basis(dyns, thresh.Get(), nBasis.Get(), demean, rotate, normalize);
   HD5::Writer writer(oname.Get());
-  basis.write(writer);
+  Basis(dyns, thresh.Get(), nBasis.Get(), demean, rotate, normalize, writer);
   writer.writeMatrix(pars, HD5::Keys::Parameters);
   return EXIT_SUCCESS;
 }