make a header-only reader

Author: mattcieslak

cpp/cmd/connectome2tck.cpp

@@ -220,7 +220,7 @@ void usage() {
 //     1-based indices, matching trx2connectome's ordering convention.
 static std::vector<std::vector<node_t>>
 assignments_from_trx_groups(const std::string &trx_path, const std::string &prefix_filter, node_t &max_node_index) {
-  auto trx = load_trx(trx_path);
+  auto trx = load_trx_header_only(trx_path);
   if (!trx || !trx->streamlines)
     throw Exception("Failed to load TRX file: " + trx_path);
   if (trx->groups.empty())

cpp/cmd/tckconvert.cpp

@@ -139,7 +139,13 @@ void usage() {
       .allow_multiple()
       + Argument ("name").type_text()
       + Argument ("datatype").type_text()
-      + Argument ("path").type_file_in();
+      + Argument ("path").type_file_in()
+
+    + Option ("positions_datatype",
+        "datatype of the positions array in the output TRX file (float16, float32, or float64). "
+        "Only applies to TRX output. Default: preserves the source datatype for TRX→TRX conversion; "
+        "float32 for all other input formats.")
+      + Argument ("spec").type_text();
 
 }
 // clang-format on
@@ -832,7 +838,23 @@ void run() {
       trx->add_dpv_from_tsf(spec.name, spec.dtype, spec.path);
     }
     try {
-      trx->save(trx_save_path, ZIP_CM_STORE);
+      trx::TrxSaveOptions save_opts;
+      save_opts.compression_standard = ZIP_CM_STORE;
+      {
+        auto opt = get_options("positions_datatype");
+        if (!opt.empty()) {
+          const std::string spec = opt[0][0];
+          if (spec == "float16")
+            save_opts.positions_dtype = trx::TrxScalarType::Float16;
+          else if (spec == "float32")
+            save_opts.positions_dtype = trx::TrxScalarType::Float32;
+          else if (spec == "float64")
+            save_opts.positions_dtype = trx::TrxScalarType::Float64;
+          else
+            throw Exception("Unknown -positions_datatype '" + spec + "'; expected float16, float32, or float64");
+        }
+      }
+      trx->save(trx_save_path, save_opts);
       if (rename_trx_directory) {
         std::error_code ec;
         std::filesystem::remove_all(trx_output, ec);

cpp/cmd/tckinfo.cpp

@@ -56,7 +56,7 @@ void run() {
     std::cout << "  Tracks file: \"" << argument[i] << "\"\n";
 
     if (Tractography::TRX::is_trx(argument[i])) {
-      auto trx = Tractography::TRX::load_trx(argument[i]);
+      auto trx = Tractography::TRX::load_trx_header_only(argument[i]);
       if (!trx)
         throw Exception("Failed to load TRX file: " + std::string(argument[i]));
       Tractography::TRX::print_info(std::cout, *trx, prefix_depth, !prefix_depth_specified);

cpp/cmd/trx2connectome.cpp

@@ -193,7 +193,7 @@ void run() {
   if (!is_trx(tracks_path))
     throw Exception("Input must be a TRX file; use trxlabel first to add group assignments");
 
-  auto trx = load_trx(tracks_path);
+  auto trx = load_trx_header_only(tracks_path);
   if (!trx)
     throw Exception("Failed to load TRX file: " + tracks_path);
 

cpp/cmd/tsfinfo.cpp

@@ -67,7 +67,7 @@ void run() {
 
     if (DWI::Tractography::TRX::is_trx(path)) {
       // TRX mode: list dpv fields
-      auto trx = DWI::Tractography::TRX::load_trx(path);
+      auto trx = DWI::Tractography::TRX::load_trx_header_only(path);
       if (!trx)
         throw Exception("Failed to load TRX file: " + path);
 

cpp/cmd/tsfvalidate.cpp

@@ -64,7 +64,7 @@ void run() {
       throw Exception("Use -field to specify the dpv field to validate when the tracks argument is a TRX file");
 
     const std::string field_name(field_opt[0][0]);
-    auto trx = TRX::load_trx(tck_path);
+    auto trx = TRX::load_trx_header_only(tck_path);
     if (!trx)
       throw Exception("Failed to load TRX file: " + tck_path);
 

cpp/core/dwi/tractography/trx_utils.h

@@ -59,48 +59,33 @@ inline bool is_trx_field_name(std::string_view input_path, std::string_view outp
          output_arg.find('\\') == std::string_view::npos;
 }
 
-// Load a TRX file, auto-detecting zip vs directory format, and expose positions
-// as float32.
-//
-// For non-float32 positions, conversion is performed in bounded chunks directly
-// into the loaded float32 matrix to avoid allocating a second full-size copy.
+// Load a TRX file with positions guaranteed as float32.
+// For float16/float64 files, delegates to trx::load_float32_positions() which
+// performs chunked in-place conversion into owned storage — no segfault, no
+// double-size allocation.
+// Use this whenever streamline coordinates will be read or written.
 inline std::unique_ptr<trx::TrxFile<float>> load_trx(std::string_view path) {
   const std::string filename(path);
   const auto dtype = trx::detect_positions_scalar_type(filename, trx::TrxScalarType::Float32);
-  if (dtype == trx::TrxScalarType::Float32)
-    return trx::load<float>(filename);
-  if (dtype == trx::TrxScalarType::Float16) {
+  if (dtype == trx::TrxScalarType::Float16)
     WARN("TRX file '" + filename + "' has float16 positions; converting to float32 on load");
-  } else {
+  else if (dtype == trx::TrxScalarType::Float64)
     WARN("TRX file '" + filename + "' has float64 positions; converting to float32 (precision loss)");
-  }
-  auto trx_f32 = trx::load<float>(filename);
-  if (!trx_f32 || !trx_f32->streamlines)
-    return trx_f32;
-  const Eigen::Index n_vertices = static_cast<Eigen::Index>(trx_f32->num_vertices());
-  if (n_vertices <= 0)
-    return trx_f32;
-  const Eigen::Index chunk_rows = static_cast<Eigen::Index>(1) << 20;
-  trx::with_trx_reader(filename, [&](auto &reader, trx::TrxScalarType) -> int {
-    const auto &src = reader->streamlines->_data;
-    if (src.rows() != n_vertices || src.cols() != 3)
-      throw Exception("unexpected TRX positions shape while converting to float32");
-    for (Eigen::Index row0 = 0; row0 < n_vertices; row0 += chunk_rows) {
-      const Eigen::Index row1 = std::min<Eigen::Index>(n_vertices, row0 + chunk_rows);
-      for (Eigen::Index r = row0; r < row1; ++r) {
-        trx_f32->streamlines->_data(r, 0) = static_cast<float>(src(r, 0));
-        trx_f32->streamlines->_data(r, 1) = static_cast<float>(src(r, 1));
-        trx_f32->streamlines->_data(r, 2) = static_cast<float>(src(r, 2));
-      }
-    }
-    return 0;
-  });
-  return trx_f32;
+  return trx::load_float32_positions(filename);
+}
+
+// Load a TRX file for metadata inspection only (groups, dps, dpv field
+// names/shapes, streamline/vertex counts). Positions are mmapped with their
+// native dtype but never read, so this is safe for float16/float64 files and
+// avoids any conversion overhead.
+// Do NOT read streamline coordinates from the returned TrxFile; use load_trx().
+inline std::unique_ptr<trx::TrxFile<float>> load_trx_header_only(std::string_view path) {
+  return trx::load<float>(std::string(path));
 }
 
 // Count streamlines and total vertices in a TRX file
 inline std::pair<size_t, size_t> count_trx(std::string_view path) {
-  auto trx = load_trx(path);
+  auto trx = load_trx_header_only(path);
   if (!trx)
     throw Exception("Failed to load TRX file: " + std::string(path));
   return {trx->num_streamlines(), trx->num_vertices()};
@@ -338,7 +323,7 @@ append_group(std::string_view trx_path, const std::string &name, const std::vect
 // because TrxFile<float> already presents all values as float.
 // Group indices are copied as uint32_t.
 inline void copy_trx_sidecar_data(std::string_view src_path, std::string_view dst_path, bool include_dpv = false) {
-  auto src = load_trx(src_path);
+  auto src = load_trx_header_only(src_path);
   if (!src)
     return;
   const std::string dst(dst_path);
@@ -758,7 +743,7 @@ inline std::vector<float> resolve_dps_weights(std::string_view tractogram_path,
   }
   if (!is_trx(tractogram_path))
     throw Exception("cannot resolve \"" + field_name_or_path + "\": not an existing file and input is not a TRX file");
-  auto trx = load_trx(tractogram_path);
+  auto trx = load_trx_header_only(tractogram_path);
   auto it = trx->data_per_streamline.find(field_name_or_path);
   if (it == trx->data_per_streamline.end() || !it->second)
     throw Exception("TRX file has no dps field named \"" + field_name_or_path + "\"");
@@ -818,7 +803,7 @@ inline std::vector<float> resolve_dpv_scalars(std::string_view tractogram_path,
   }
   if (!is_trx(tractogram_path))
     throw Exception("cannot resolve \"" + field_name_or_path + "\": not an existing file and input is not a TRX file");
-  auto trx = load_trx(tractogram_path);
+  auto trx = load_trx_header_only(tractogram_path);
   auto it = trx->data_per_vertex.find(field_name_or_path);
   if (it == trx->data_per_vertex.end() || !it->second)
     throw Exception("TRX file has no dpv field named \"" + field_name_or_path + "\"");
@@ -838,7 +823,7 @@ inline std::vector<float> resolve_dpv_scalars(std::string_view tractogram_path,
 class TRXScalarReader {
 public:
   TRXScalarReader(std::string_view trx_path, const std::string &field_name) : current_(0) {
-    auto trx = load_trx(trx_path);
+    auto trx = load_trx_header_only(trx_path);
     if (!trx || !trx->streamlines)
       throw Exception("Failed to load TRX file: " + std::string(trx_path));