diff --git a/src/cmd/basis-img.cpp b/src/cmd/basis-img.cpp index 48e1d023..01801b26 100644 --- 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 iname(parser, "INPUT", "Input image file"); args::Positional oname(parser, "OUTPUT", "Name for the basis file"); - args::ValueFlag ss(parser, "S", "Subsample image", {"subsamp"}, 1); + args::Flag otsu(parser, "O", "Otsu mask", {"otsu"}); + args::ValueFlag st(parser, "S", "ROI Start", {"roi-start"}); + args::ValueFlag sz(parser, "S", "ROI size", {"roi-size"}); args::ValueFlag spf(parser, "S", "Spokes per frame", {"spf"}, 1); args::ValueFlag 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(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;