Add ROI selection for basis-img

spinicist · May 16, 2023 · ab4ebc7 · ab4ebc7
1 parent 2541693
commit ab4ebc7
Showing 1 changed file with 17 additions and 6 deletions.
diff --git a/src/cmd/basis-img.cpp b/src/cmd/basis-img.cpp
@@ -17,7 +17,9 @@ 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::Flag otsu(parser, "O", "Otsu mask", {"otsu"});
+  args::ValueFlag<Sz3, Sz3Reader> st(parser, "S", "ROI Start", {"roi-start"});
+  args::ValueFlag<Sz3, Sz3Reader> sz(parser, "S", "ROI size", {"roi-size"});
   args::ValueFlag<Index> spf(parser, "S", "Spokes per frame", {"spf"}, 1);
   args::ValueFlag<Index> order(parser, "O", "Interpolation order", {"interp-order"}, 3);
   args::Flag clamp(parser, "C", "Clamp interpolation", {"interp-clamp"});
@@ -33,22 +35,31 @@ int main_basis_img(args::Subparser &parser)
 
   HD5::Reader reader(iname.Get());
   Cx4 img = reader.readSlab<Cx4>(HD5::Keys::Image, 0);
-  if (ss) {
-    img = Cx4(img.stride(Sz4{1, ss.Get(), ss.Get(), ss.Get()}));
+  if (st && sz) {
+    img = Cx4(img.slice(AddFront(st.Get(), 0), AddFront(sz.Get(), img.dimension(0))));
   }
   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();
+  Re4 const real = (img / (ref / ref.abs()).broadcast(Sz4{shape[0], 1, 1, 1})).real();
   auto const realMat = CollapseToMatrix(real);
   Re3 const toMask = ref.chip<0>(0).abs();
   auto const toMaskMat = CollapseToArray(toMask);
-  auto const [thresh, count] = Otsu(toMaskMat);
+  float thresh;
+  Index count;
+  if (otsu) {
+    auto o = Otsu(toMaskMat);
+    thresh = o.thresh;
+    count = o.countAbove;
+  } else {
+    thresh = 0.f;
+    count = Product(LastN<3>(shape));
+  }
   Index const f = shape[0];
   Index const s = shape[0] * spf.Get();
   Eigen::MatrixXf dynamics(s, count);
   Index col = 0;
   Eigen::ArrayXi x1, x2, z1, z2;
-  Index f1, f2, s1, s2;
+  Index f1 = 0, f2 = 0, s1 = 0, s2 = 0;
   if (start2) {
     f1 = start2.Get();
     f2 = f - f1;