From 516205dbcf26bbea40b6f8bca06387badf876b79 Mon Sep 17 00:00:00 2001 From: frheault Date: Wed, 17 Jun 2026 14:22:06 -0400 Subject: [PATCH 01/11] Casting routine from files/header --- src/trx.cpp | 50 ++++++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 46 insertions(+), 4 deletions(-) diff --git a/src/trx.cpp b/src/trx.cpp index 869a734..380ab89 100644 --- a/src/trx.cpp +++ b/src/trx.cpp @@ -512,12 +512,54 @@ AnyTrxFile::_create_from_pointer(json header, if (dim != 1) { throw TrxFormatError("Wrong group dimensionality"); } - if (ext != "uint32") { + if (ext == "uint32") { + auto arr = make_typed_array(elem_filename, static_cast(size), 1, ext); + arr.materialize_to_owned(); + trx.groups.emplace(base, std::move(arr)); + } else if (ext == "int64" || ext == "uint64" || ext == "int32" || ext == "uint8" || ext == "int8" || ext == "uint16" || ext == "int16") { + if (ext == "int32" || ext == "uint8" || ext == "int8" || ext == "uint16" || ext == "int16") { + std::cerr << "Warning: Upcasting group from " << ext << " to uint32\n"; + } + if (ext == "int64" || ext == "uint64") { + uint64_t num_strs = static_cast(header["NB_STREAMLINES"].number_value()); + if (num_strs > 4294967295ULL) { + throw TrxFormatError("downcasting is unsafe because the number of streamlines exceeds the 32-bit limit"); + } + } + auto tmp_arr = make_typed_array(elem_filename, static_cast(size), 1, ext); + tmp_arr.materialize_to_owned(); + TypedArray arr; + arr.dtype = "uint32"; + arr.rows = static_cast(size); + arr.cols = 1; + arr.owned.resize(static_cast(size) * sizeof(uint32_t)); + uint32_t* dst = reinterpret_cast(arr.owned.data()); + if (ext == "int64") { + const int64_t* src = reinterpret_cast(tmp_arr.owned.data()); + for (size_t i = 0; i < size; ++i) dst[i] = static_cast(src[i]); + } else if (ext == "uint64") { + const uint64_t* src = reinterpret_cast(tmp_arr.owned.data()); + for (size_t i = 0; i < size; ++i) dst[i] = static_cast(src[i]); + } else if (ext == "uint8") { + const uint8_t* src = reinterpret_cast(tmp_arr.owned.data()); + for (size_t i = 0; i < size; ++i) dst[i] = static_cast(src[i]); + } else if (ext == "int8") { + const int8_t* src = reinterpret_cast(tmp_arr.owned.data()); + for (size_t i = 0; i < size; ++i) dst[i] = static_cast(src[i]); + } else if (ext == "uint16") { + const uint16_t* src = reinterpret_cast(tmp_arr.owned.data()); + for (size_t i = 0; i < size; ++i) dst[i] = static_cast(src[i]); + } else if (ext == "int16") { + const int16_t* src = reinterpret_cast(tmp_arr.owned.data()); + for (size_t i = 0; i < size; ++i) dst[i] = static_cast(src[i]); + } else if (ext == "int32") { + const int32_t* src = reinterpret_cast(tmp_arr.owned.data()); + for (size_t i = 0; i < size; ++i) dst[i] = static_cast(src[i]); + } + trx.groups.emplace(base, std::move(arr)); + } else { throw TrxDTypeError("Unsupported group dtype: " + ext); } - auto arr = make_typed_array(elem_filename, static_cast(size), 1, ext); - arr.materialize_to_owned(); - trx.groups.emplace(base, std::move(arr)); } else { throw TrxFormatError("Entry is not part of a valid TRX structure: " + elem_filename); } From d5e182a9e036564b6e1462da1eacfc33d31d6be7 Mon Sep 17 00:00:00 2001 From: frheault Date: Mon, 22 Jun 2026 10:32:46 -0400 Subject: [PATCH 02/11] Introduce legacy IO namespace to parse and export TRK, TCK, and VTK files without metadata loss --- CMakeLists.txt | 1 + include/trx/legacy_io.h | 60 +++++ src/legacy_io.cpp | 550 ++++++++++++++++++++++++++++++++++++++++ 3 files changed, 611 insertions(+) create mode 100644 include/trx/legacy_io.h create mode 100644 src/legacy_io.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 0897a8b..9e7f30f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -174,6 +174,7 @@ endif() # ── Core library ──────────────────────────────────────────────────────────── add_library(trx src/trx.cpp + src/legacy_io.cpp src/detail/dtype_helpers.cpp include/trx/trx.h include/trx/trx.tpp diff --git a/include/trx/legacy_io.h b/include/trx/legacy_io.h new file mode 100644 index 0000000..53408cb --- /dev/null +++ b/include/trx/legacy_io.h @@ -0,0 +1,60 @@ +#ifndef TRX_LEGACY_IO_H +#define TRX_LEGACY_IO_H + +#include +#include +#include +#include + +namespace trx { +namespace legacy { + +struct Tractogram { + std::vector pts; + std::vector offsets; + json11::Json header; + std::shared_ptr original_trx; +}; + +#pragma pack(push, 1) +struct TrkHeader { + char magic_number[6]; + int16_t dimensions[3]; + float voxel_sizes[3]; + float origin[3]; + int16_t nb_scalars_per_point; + char scalar_name[10][20]; + int16_t nb_properties_per_streamline; + char property_name[10][20]; + float voxel_to_rasmm[4][4]; + char reserved[444]; + char voxel_order[4]; + char pad2[4]; + float image_orientation_patient[6]; + char pad1[2]; + char invert_x; + char invert_y; + char invert_z; + char swap_xy; + char swap_yz; + char swap_zx; + int32_t nb_streamlines; + int32_t version; + int32_t hdr_size; +}; +#pragma pack(pop) + +bool load_trx(const std::string &filename, Tractogram &tr); +bool load_trk(const std::string &filename, Tractogram &tr); +bool load_tck(const std::string &filename, Tractogram &tr); +bool load_vtk(const std::string &filename, Tractogram &tr); + +bool save_trx(const Tractogram &tr, const std::string &out_path); +bool save_trk(const Tractogram &tr, const std::string &out_path, const std::string &original_filename = ""); +bool save_tck(const Tractogram &tr, const std::string &out_path); +bool save_vtk(const Tractogram &tr, const std::string &out_path); + +} // namespace legacy +} // namespace trx + +#endif // TRX_LEGACY_IO_H diff --git a/src/legacy_io.cpp b/src/legacy_io.cpp new file mode 100644 index 0000000..205421d --- /dev/null +++ b/src/legacy_io.cpp @@ -0,0 +1,550 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace trx { +namespace legacy { + +inline float swap_float(float f) { + union { + float f; + uint32_t i; + } u; + u.f = f; + u.i = __builtin_bswap32(u.i); + return u.f; +} + +inline int32_t swap_int32(int32_t i) { + return __builtin_bswap32(i); +} + + + +bool load_trx(const std::string &filename, Tractogram &tr) { + try { + auto trx = trx::AnyTrxFile::load(filename); + size_t num_streamlines = trx.num_streamlines(); + size_t num_points = trx.num_vertices(); + + tr.pts.resize(num_points * 3); + tr.offsets.resize(num_streamlines + 1); + tr.header = trx.header; + + // Load offsets + if (!trx.offsets.empty()) { + if (trx.offsets.dtype == "uint32") { + auto mat = trx.offsets.as_matrix(); + for (size_t i = 0; i <= num_streamlines; ++i) tr.offsets[i] = mat.data()[i]; + } else if (trx.offsets.dtype == "uint64") { + auto mat = trx.offsets.as_matrix(); + for (size_t i = 0; i <= num_streamlines; ++i) tr.offsets[i] = mat.data()[i]; + } else if (trx.offsets.dtype == "int32") { + auto mat = trx.offsets.as_matrix(); + for (size_t i = 0; i <= num_streamlines; ++i) tr.offsets[i] = mat.data()[i]; + } else if (trx.offsets.dtype == "int64") { + auto mat = trx.offsets.as_matrix(); + for (size_t i = 0; i <= num_streamlines; ++i) tr.offsets[i] = mat.data()[i]; + } + } + + // Load positions quickly (bulk copy / fast casting) + if (!trx.positions.empty()) { + if (trx.positions.dtype == "float32") { + auto mat = trx.positions.as_matrix(); + std::memcpy(tr.pts.data(), mat.data(), num_points * 3 * sizeof(float)); + } else if (trx.positions.dtype == "float16") { + auto mat = trx.positions.as_matrix(); + for (size_t i = 0; i < num_points * 3; ++i) { + tr.pts[i] = static_cast(mat.data()[i]); + } + } else if (trx.positions.dtype == "float64") { + auto mat = trx.positions.as_matrix(); + for (size_t i = 0; i < num_points * 3; ++i) { + tr.pts[i] = static_cast(mat.data()[i]); + } + } + } + + tr.original_trx = std::make_shared(std::move(trx)); + return true; + } catch (const std::exception &e) { + std::cerr << "Error loading TRX file: " << e.what() << std::endl; + return false; + } +} + +bool load_trk(const std::string &filename, Tractogram &tr) { + std::ifstream f(filename, std::ios::binary | std::ios::ate); + if (!f.is_open()) return false; + + std::streamsize size = f.tellg(); + f.seekg(0, std::ios::beg); + + std::vector buffer(size); + if (!f.read(buffer.data(), size)) return false; + if (buffer.size() < 1000) return false; + + const TrkHeader* header = reinterpret_cast(buffer.data()); + if (std::string(header->magic_number, 5) != "TRACK") return false; + + int16_t n_scalars = header->nb_scalars_per_point; + int16_t n_properties = header->nb_properties_per_streamline; + + // Store metadata + tr.header = json11::Json::object { + { "DIMENSIONS", json11::Json::array { header->dimensions[0], header->dimensions[1], header->dimensions[2] } }, + { "VOXEL_TO_RASMM", json11::Json::array { + json11::Json::array { header->voxel_to_rasmm[0][0], header->voxel_to_rasmm[0][1], header->voxel_to_rasmm[0][2], header->voxel_to_rasmm[0][3] }, + json11::Json::array { header->voxel_to_rasmm[1][0], header->voxel_to_rasmm[1][1], header->voxel_to_rasmm[1][2], header->voxel_to_rasmm[1][3] }, + json11::Json::array { header->voxel_to_rasmm[2][0], header->voxel_to_rasmm[2][1], header->voxel_to_rasmm[2][2], header->voxel_to_rasmm[2][3] }, + json11::Json::array { header->voxel_to_rasmm[3][0], header->voxel_to_rasmm[3][1], header->voxel_to_rasmm[3][2], header->voxel_to_rasmm[3][3] } + } } + }; + + tr.offsets.clear(); + tr.offsets.push_back(0); + tr.pts.clear(); + + size_t offset = 1000; + while (offset + sizeof(int32_t) <= buffer.size()) { + int32_t n_points = *reinterpret_cast(buffer.data() + offset); + offset += sizeof(int32_t); + + tr.offsets.push_back(tr.offsets.back() + n_points); + + for (int32_t j = 0; j < n_points; ++j) { + float x = *reinterpret_cast(buffer.data() + offset); + float y = *reinterpret_cast(buffer.data() + offset + 4); + float z = *reinterpret_cast(buffer.data() + offset + 8); + tr.pts.push_back(x); + tr.pts.push_back(y); + tr.pts.push_back(z); + + offset += (3 + n_scalars) * sizeof(float); + } + offset += n_properties * sizeof(float); + } + + return true; +} + +bool load_tck(const std::string &filename, Tractogram &tr) { + std::ifstream f(filename, std::ios::binary | std::ios::ate); + if (!f.is_open()) return false; + + std::streamsize size = f.tellg(); + f.seekg(0, std::ios::beg); + + std::vector buffer(size); + if (!f.read(buffer.data(), size)) return false; + + std::string_view view(buffer.data(), buffer.size()); + size_t file_pos = view.find("file: . "); + if (file_pos == std::string_view::npos) return false; + size_t offset_pos = file_pos + 8; + size_t offset_end = view.find_first_not_of("0123456789", offset_pos); + if (offset_end == std::string_view::npos) return false; + size_t offset = std::stoull(std::string(view.substr(offset_pos, offset_end - offset_pos))); + + if (offset >= buffer.size()) return false; + + const float* data = reinterpret_cast(buffer.data() + offset); + size_t num_floats = (buffer.size() - offset) / sizeof(float); + size_t num_triplets = num_floats / 3; + + tr.offsets.clear(); + tr.offsets.push_back(0); + tr.pts.clear(); + + bool in_streamline = false; + size_t current_pts = 0; + + for (size_t i = 0; i < num_triplets; ++i) { + float x = data[i * 3]; + float y = data[i * 3 + 1]; + float z = data[i * 3 + 2]; + + if (std::isinf(x) && std::isinf(y) && std::isinf(z)) { + if (in_streamline) { + tr.offsets.push_back(tr.offsets.back() + current_pts); + current_pts = 0; + in_streamline = false; + } + break; + } else if (std::isnan(x) && std::isnan(y) && std::isnan(z)) { + if (in_streamline) { + tr.offsets.push_back(tr.offsets.back() + current_pts); + current_pts = 0; + in_streamline = false; + } + } else { + in_streamline = true; + tr.pts.push_back(x); + tr.pts.push_back(y); + tr.pts.push_back(z); + current_pts++; + } + } + + if (in_streamline) { + tr.offsets.push_back(tr.offsets.back() + current_pts); + } + + return true; +} + +bool load_vtk(const std::string &filename, Tractogram &tr) { + std::ifstream f(filename, std::ios::binary); + if (!f.is_open()) return false; + + std::string line; + size_t num_points = 0; + bool is_double = false; + while (std::getline(f, line)) { + if (line.rfind("POINTS ", 0) == 0) { + size_t space1 = line.find(" ", 7); + num_points = std::stoull(line.substr(7, space1 - 7)); + if (line.find("double", space1) != std::string::npos) { + is_double = true; + } + break; + } + } + if (num_points == 0) return false; + + tr.pts.resize(num_points * 3); + if (is_double) { + std::vector dpts(num_points * 3); + f.read(reinterpret_cast(dpts.data()), num_points * 3 * sizeof(double)); + for (size_t i = 0; i < num_points * 3; ++i) { + uint64_t val; + std::memcpy(&val, &dpts[i], 8); + val = ((val & 0xFF00000000000000ULL) >> 56) | ((val & 0x00FF000000000000ULL) >> 40) | + ((val & 0x0000FF0000000000ULL) >> 24) | ((val & 0x000000FF00000000ULL) >> 8) | + ((val & 0x00000000FF000000ULL) << 8) | ((val & 0x0000000000FF0000ULL) << 24) | + ((val & 0x000000000000FF00ULL) << 40) | ((val & 0x00000000000000FFULL) << 56); + double swapped; + std::memcpy(&swapped, &val, 8); + tr.pts[i] = static_cast(swapped); + } + } else { + f.read(reinterpret_cast(tr.pts.data()), num_points * 3 * sizeof(float)); + for (size_t i = 0; i < num_points * 3; ++i) { + tr.pts[i] = swap_float(tr.pts[i]); + } + } + + size_t num_streamlines = 0; + while (std::getline(f, line)) { + if (line.rfind("LINES ", 0) == 0) { + num_streamlines = std::stoull(line.substr(6, line.find(" ", 6) - 6)); + break; + } + } + if (num_streamlines == 0) return false; + + auto pos_before_offsets = f.tellg(); + std::getline(f, line); + if (!line.empty() && line.back() == '\r') line.pop_back(); + bool has_offsets = (line.rfind("OFFSETS", 0) == 0); + bool is_int64 = (line.find("int64") != std::string::npos); + + if (has_offsets) { + tr.offsets.resize(num_streamlines); + for (size_t i = 0; i < num_streamlines; ++i) { + if (is_int64) { + uint64_t val; + f.read(reinterpret_cast(&val), 8); + val = ((val & 0xFF00000000000000ULL) >> 56) | ((val & 0x00FF000000000000ULL) >> 40) | + ((val & 0x0000FF0000000000ULL) >> 24) | ((val & 0x000000FF00000000ULL) >> 8) | + ((val & 0x00000000FF000000ULL) << 8) | ((val & 0x0000000000FF0000ULL) << 24) | + ((val & 0x000000000000FF00ULL) << 40) | ((val & 0x00000000000000FFULL) << 56); + tr.offsets[i] = val; + } else { + uint32_t val; + f.read(reinterpret_cast(&val), 4); + val = swap_int32(val); + tr.offsets[i] = val; + } + } + return true; + } + f.seekg(pos_before_offsets); + + tr.offsets.clear(); + tr.offsets.push_back(0); + + for (size_t i = 0; i < num_streamlines; ++i) { + int32_t n_pts; + f.read(reinterpret_cast(&n_pts), sizeof(int32_t)); + if (!f) break; + n_pts = swap_int32(n_pts); + if (n_pts == 0) continue; + tr.offsets.push_back(tr.offsets.back() + n_pts); + + // Skip cell indices + f.seekg(n_pts * sizeof(int32_t), std::ios::cur); + } + + return true; +} + +bool save_trx(const Tractogram &tr, const std::string &out_path) { + try { + if (tr.original_trx) { + tr.original_trx->save(out_path, trx::TrxCompression::None); + return true; + } + size_t nb_vertices = tr.pts.size() / 3; + size_t nb_streamlines = tr.offsets.size() - 1; + + trx::TrxFile trx(nb_vertices, nb_streamlines); + + // Copy positions + std::memcpy(trx.streamlines->_data.data(), tr.pts.data(), tr.pts.size() * sizeof(float)); + + // Copy offsets + for (size_t i = 0; i <= nb_streamlines; ++i) { + trx.streamlines->_offsets(i, 0) = tr.offsets[i]; + } + + // Compute lengths + for (size_t i = 0; i < nb_streamlines; ++i) { + trx.streamlines->_lengths(i, 0) = tr.offsets[i+1] - tr.offsets[i]; + } + + // Copy header + trx.header = tr.header; + + trx.save(out_path, trx::TrxCompression::None); + trx.close(); + + return true; + } catch (const std::exception &e) { + std::cerr << "Error saving TRX file: " << e.what() << std::endl; + return false; + } +} + +bool save_trk(const Tractogram &tr, const std::string &out_path, const std::string &original_filename) { + std::ofstream f(out_path, std::ios::binary); + if (!f.is_open()) return false; + + TrkHeader header; + std::memset(&header, 0, sizeof(header)); + std::memcpy(header.magic_number, "TRACK", 5); + + // Default dimensions, voxel sizes and affine + header.dimensions[0] = 256; header.dimensions[1] = 256; header.dimensions[2] = 256; + header.voxel_sizes[0] = 1.0f; header.voxel_sizes[1] = 1.0f; header.voxel_sizes[2] = 1.0f; + for (int r = 0; r < 4; ++r) { + for (int c = 0; c < 4; ++c) { + header.voxel_to_rasmm[r][c] = (r == c) ? 1.0f : 0.0f; + } + } + + // Attempt to extract from JSON header + if (tr.header["DIMENSIONS"].is_array()) { + auto dims = tr.header["DIMENSIONS"].array_items(); + if (dims.size() >= 3) { + header.dimensions[0] = static_cast(dims[0].number_value()); + header.dimensions[1] = static_cast(dims[1].number_value()); + header.dimensions[2] = static_cast(dims[2].number_value()); + } + } + if (tr.header["VOXEL_TO_RASMM"].is_array()) { + auto rows = tr.header["VOXEL_TO_RASMM"].array_items(); + if (rows.size() >= 4) { + float vox_to_ras[4][4]; + for (int r = 0; r < 4; ++r) { + auto cols = rows[r].array_items(); + if (cols.size() >= 4) { + for (int c = 0; c < 4; ++c) { + vox_to_ras[r][c] = static_cast(cols[c].number_value()); + header.voxel_to_rasmm[r][c] = vox_to_ras[r][c]; + } + } + } + header.voxel_sizes[0] = std::sqrt(vox_to_ras[0][0]*vox_to_ras[0][0] + vox_to_ras[1][0]*vox_to_ras[1][0] + vox_to_ras[2][0]*vox_to_ras[2][0]); + header.voxel_sizes[1] = std::sqrt(vox_to_ras[0][1]*vox_to_ras[0][1] + vox_to_ras[1][1]*vox_to_ras[1][1] + vox_to_ras[2][1]*vox_to_ras[2][1]); + header.voxel_sizes[2] = std::sqrt(vox_to_ras[0][2]*vox_to_ras[0][2] + vox_to_ras[1][2]*vox_to_ras[1][2] + vox_to_ras[2][2]*vox_to_ras[2][2]); + } + } + + std::memcpy(header.voxel_order, "RAS", 3); + header.nb_streamlines = static_cast(tr.offsets.size() - 1); + header.version = 2; + header.hdr_size = 1000; + + f.write(reinterpret_cast(&header), 1000); + + size_t num_streamlines = tr.offsets.size() - 1; + std::vector chunk; + chunk.reserve(4 * 1024 * 1024); + + for (size_t i = 0; i < num_streamlines; ++i) { + size_t start = tr.offsets[i]; + size_t end = tr.offsets[i+1]; + int32_t n_pts = static_cast(end - start); + + // Push n_pts + const char* p_n_pts = reinterpret_cast(&n_pts); + chunk.insert(chunk.end(), p_n_pts, p_n_pts + 4); + + // Push points + for (size_t j = start; j < end; ++j) { + float x = tr.pts[j*3]; + float y = tr.pts[j*3 + 1]; + float z = tr.pts[j*3 + 2]; + const char* px = reinterpret_cast(&x); + const char* py = reinterpret_cast(&y); + const char* pz = reinterpret_cast(&z); + chunk.insert(chunk.end(), px, px + 4); + chunk.insert(chunk.end(), py, py + 4); + chunk.insert(chunk.end(), pz, pz + 4); + } + + if (chunk.size() >= 4000000) { + f.write(chunk.data(), chunk.size()); + chunk.clear(); + } + } + + if (!chunk.empty()) { + f.write(chunk.data(), chunk.size()); + } + + return true; +} + +bool save_tck(const Tractogram &tr, const std::string &out_path) { + std::ofstream f(out_path, std::ios::binary); + if (!f.is_open()) return false; + + size_t num_streamlines = tr.offsets.size() - 1; + + // Build TCK header + std::string header; + size_t offset = 80; + while (true) { + char buf[256]; + snprintf(buf, sizeof(buf), "mrtrix tracks\ncount: %010zu\ndatatype: Float32LE\nfile: . %zu\nEND\n", num_streamlines, offset); + std::string h(buf); + if (h.length() <= offset) { + h.append(offset - h.length(), ' '); + header = h; + break; + } + offset = h.length(); + } + f.write(header.data(), header.size()); + + // Payload writing + std::vector chunk; + chunk.reserve(1024 * 1024); + + for (size_t i = 0; i < num_streamlines; ++i) { + size_t start = tr.offsets[i]; + size_t end = tr.offsets[i+1]; + + for (size_t j = start; j < end; ++j) { + chunk.push_back(tr.pts[j*3]); + chunk.push_back(tr.pts[j*3 + 1]); + chunk.push_back(tr.pts[j*3 + 2]); + if (chunk.size() >= 1000000) { + f.write(reinterpret_cast(chunk.data()), chunk.size() * sizeof(float)); + chunk.clear(); + } + } + // Delimiter + chunk.push_back(std::numeric_limits::quiet_NaN()); + chunk.push_back(std::numeric_limits::quiet_NaN()); + chunk.push_back(std::numeric_limits::quiet_NaN()); + } + + // EOF Delimiter + chunk.push_back(std::numeric_limits::infinity()); + chunk.push_back(std::numeric_limits::infinity()); + chunk.push_back(std::numeric_limits::infinity()); + + if (!chunk.empty()) { + f.write(reinterpret_cast(chunk.data()), chunk.size() * sizeof(float)); + } + + return true; +} + +bool save_vtk(const Tractogram &tr, const std::string &out_path) { + std::ofstream f(out_path, std::ios::binary); + if (!f.is_open()) return false; + + size_t num_streamlines = tr.offsets.size() - 1; + size_t num_points = tr.pts.size() / 3; + + // Write ASCII header + char header[512]; + snprintf(header, sizeof(header), "# vtk DataFile Version 3.0\nvtk output\nBINARY\nDATASET POLYDATA\nPOINTS %zu float\n", num_points); + f.write(header, std::strlen(header)); + + // Write POINTS binary block (big-endian floats) + std::vector pts_buf; + pts_buf.reserve(1024 * 1024); + + for (size_t i = 0; i < num_points * 3; ++i) { + pts_buf.push_back(swap_float(tr.pts[i])); + if (pts_buf.size() >= 1000000) { + f.write(reinterpret_cast(pts_buf.data()), pts_buf.size() * sizeof(float)); + pts_buf.clear(); + } + } + if (!pts_buf.empty()) { + f.write(reinterpret_cast(pts_buf.data()), pts_buf.size() * sizeof(float)); + } + + // Write LINES header + size_t cell_array_size = num_streamlines + num_points; + char lines_hdr[128]; + snprintf(lines_hdr, sizeof(lines_hdr), "LINES %zu %zu\n", num_streamlines, cell_array_size); + f.write(lines_hdr, std::strlen(lines_hdr)); + + // Write LINES binary block (big-endian int32) + std::vector lines_buf; + lines_buf.reserve(1024 * 1024); + + int32_t current_point_idx = 0; + for (size_t i = 0; i < num_streamlines; ++i) { + size_t start = tr.offsets[i]; + size_t end = tr.offsets[i+1]; + int32_t n_pts = static_cast(end - start); + + lines_buf.push_back(swap_int32(n_pts)); + for (int32_t j = 0; j < n_pts; ++j) { + lines_buf.push_back(swap_int32(current_point_idx++)); + } + + if (lines_buf.size() >= 1000000) { + f.write(reinterpret_cast(lines_buf.data()), lines_buf.size() * sizeof(int32_t)); + lines_buf.clear(); + } + } + if (!lines_buf.empty()) { + f.write(reinterpret_cast(lines_buf.data()), lines_buf.size() * sizeof(int32_t)); + } + + return true; +} + +} // namespace legacy +} // namespace trx From 76f83e54ab9e1361f4bd28ab50e9642e8c7d80ac Mon Sep 17 00:00:00 2001 From: frheault Date: Mon, 22 Jun 2026 15:27:24 -0400 Subject: [PATCH 03/11] Added cross-endian NIfTI parsing and patched CLI argument handling edge cases --- include/trx/legacy_io.h | 6 +- main.cpp | 68 +++++++++ src/legacy_io.cpp | 326 +++++++++++++++++++++++++++++++++++++--- 3 files changed, 374 insertions(+), 26 deletions(-) create mode 100644 main.cpp diff --git a/include/trx/legacy_io.h b/include/trx/legacy_io.h index 53408cb..47e09cf 100644 --- a/include/trx/legacy_io.h +++ b/include/trx/legacy_io.h @@ -49,8 +49,10 @@ bool load_trk(const std::string &filename, Tractogram &tr); bool load_tck(const std::string &filename, Tractogram &tr); bool load_vtk(const std::string &filename, Tractogram &tr); -bool save_trx(const Tractogram &tr, const std::string &out_path); -bool save_trk(const Tractogram &tr, const std::string &out_path, const std::string &original_filename = ""); +bool load_nifti_header(const std::string &ref_path, json11::Json &out_header); + +bool save_trx(const Tractogram &tr, const std::string &out_path, const std::string &ref_nifti_path = ""); +bool save_trk(const Tractogram &tr, const std::string &out_path, const std::string &original_filename = "", const std::string &ref_nifti_path = ""); bool save_tck(const Tractogram &tr, const std::string &out_path); bool save_vtk(const Tractogram &tr, const std::string &out_path); diff --git a/main.cpp b/main.cpp new file mode 100644 index 0000000..3a01a1e --- /dev/null +++ b/main.cpp @@ -0,0 +1,68 @@ +#include +#include +#include +#include + +int main(int argc, char** argv) { + std::string input_file; + std::string output_file; + std::string ref_path; + + for (int i = 1; i < argc; ++i) { + std::string arg = argv[i]; + if (arg == "--ref") { + if (i + 1 < argc) { + ref_path = argv[++i]; + } else { + std::cerr << "Error: --ref requires an argument\n"; + return 1; + } + } else if (input_file.empty()) { + input_file = arg; + } else if (output_file.empty()) { + output_file = arg; + } + } + + if (input_file.empty() || output_file.empty()) { + std::cerr << "Usage: convert [--ref ]\n"; + return 1; + } + + trx::legacy::Tractogram tr; + bool success = false; + + auto ends_with = [](const std::string& str, const std::string& suffix) { + return str.size() >= suffix.size() && str.compare(str.size() - suffix.size(), suffix.size(), suffix) == 0; + }; + + if (ends_with(input_file, ".trx")) success = trx::legacy::load_trx(input_file, tr); + else if (ends_with(input_file, ".trk")) success = trx::legacy::load_trk(input_file, tr); + else if (ends_with(input_file, ".tck")) success = trx::legacy::load_tck(input_file, tr); + else if (ends_with(input_file, ".vtk")) success = trx::legacy::load_vtk(input_file, tr); + + if (!success) { + std::cerr << "Error loading input file\n"; + return 1; + } + + bool is_tck_vtk = ends_with(input_file, ".tck") || ends_with(input_file, ".vtk"); + bool is_trx_trk = ends_with(output_file, ".trx") || ends_with(output_file, ".trk"); + + if (is_tck_vtk && is_trx_trk && ref_path.empty()) { + std::cerr << "Error: TCK/VTK -> TRX/TRK conversion requires --ref \n"; + return 1; + } + + success = false; + if (ends_with(output_file, ".trx")) success = trx::legacy::save_trx(tr, output_file, ref_path); + else if (ends_with(output_file, ".trk")) success = trx::legacy::save_trk(tr, output_file, input_file, ref_path); + else if (ends_with(output_file, ".tck")) success = trx::legacy::save_tck(tr, output_file); + else if (ends_with(output_file, ".vtk")) success = trx::legacy::save_vtk(tr, output_file); + + if (!success) { + std::cerr << "Error saving output file\n"; + return 1; + } + return 0; +} diff --git a/src/legacy_io.cpp b/src/legacy_io.cpp index 205421d..d888c8c 100644 --- a/src/legacy_io.cpp +++ b/src/legacy_io.cpp @@ -1,5 +1,6 @@ #include #include +#include #include #include #include @@ -30,6 +31,26 @@ inline int32_t swap_int32(int32_t i) { return __builtin_bswap32(i); } +inline int16_t swap_int16(int16_t val) { + uint16_t uval = val; + uval = (uval << 8) | (uval >> 8); + return static_cast(uval); +} + +inline int64_t swap_int64(int64_t val) { + return __builtin_bswap64(val); +} + +inline double swap_double(double d) { + union { + double d; + uint64_t i; + } u; + u.d = d; + u.i = __builtin_bswap64(u.i); + return u.d; +} + bool load_trx(const std::string &filename, Tractogram &tr) { @@ -125,13 +146,26 @@ bool load_trk(const std::string &filename, Tractogram &tr) { tr.offsets.push_back(tr.offsets.back() + n_points); for (int32_t j = 0; j < n_points; ++j) { - float x = *reinterpret_cast(buffer.data() + offset); - float y = *reinterpret_cast(buffer.data() + offset + 4); - float z = *reinterpret_cast(buffer.data() + offset + 8); + float raw_x = *reinterpret_cast(buffer.data() + offset); + float raw_y = *reinterpret_cast(buffer.data() + offset + 4); + float raw_z = *reinterpret_cast(buffer.data() + offset + 8); + + float vx = header->voxel_sizes[0] > 0 ? header->voxel_sizes[0] : 1.0f; + float vy = header->voxel_sizes[1] > 0 ? header->voxel_sizes[1] : 1.0f; + float vz = header->voxel_sizes[2] > 0 ? header->voxel_sizes[2] : 1.0f; + + float cx = (raw_x / vx) - 0.5f; + float cy = (raw_y / vy) - 0.5f; + float cz = (raw_z / vz) - 0.5f; + + float x = cx * header->voxel_to_rasmm[0][0] + cy * header->voxel_to_rasmm[0][1] + cz * header->voxel_to_rasmm[0][2] + header->voxel_to_rasmm[0][3]; + float y = cx * header->voxel_to_rasmm[1][0] + cy * header->voxel_to_rasmm[1][1] + cz * header->voxel_to_rasmm[1][2] + header->voxel_to_rasmm[1][3]; + float z = cx * header->voxel_to_rasmm[2][0] + cy * header->voxel_to_rasmm[2][1] + cz * header->voxel_to_rasmm[2][2] + header->voxel_to_rasmm[2][3]; + tr.pts.push_back(x); tr.pts.push_back(y); tr.pts.push_back(z); - + offset += (3 + n_scalars) * sizeof(float); } offset += n_properties * sizeof(float); @@ -301,8 +335,198 @@ bool load_vtk(const std::string &filename, Tractogram &tr) { return true; } -bool save_trx(const Tractogram &tr, const std::string &out_path) { +bool load_nifti_header(const std::string &ref_path, json11::Json &out_header) { + std::ifstream f(ref_path, std::ios::binary); + if (!f.is_open()) { + std::cerr << "Error: Could not open reference NIfTI file: " << ref_path << "\n"; + return false; + } + char buf[540]; + f.read(buf, 540); + if (f.gcount() < 348) { + std::cerr << "Error: Invalid NIfTI file (too small)\n"; + return false; + } + + int32_t sizeof_hdr; + std::memcpy(&sizeof_hdr, buf, sizeof(int32_t)); + + bool swap_endian = false; + if (sizeof_hdr == 1543569408 || sizeof_hdr == 469893120) { + swap_endian = true; + sizeof_hdr = swap_int32(sizeof_hdr); + } + + std::vector dims(3); + float dx, dy, dz, qfac; + int sform_code, qform_code; + float srow_x[4], srow_y[4], srow_z[4]; + float qoffset_x, qoffset_y, qoffset_z, b, c, d; + + if (sizeof_hdr == 348) { // NIfTI-1 + int16_t dim[8]; + std::memcpy(dim, buf + 40, 8 * sizeof(int16_t)); + if (swap_endian) for(int i=0; i<8; i++) dim[i] = swap_int16(dim[i]); + dims[0] = dim[1]; dims[1] = dim[2]; dims[2] = dim[3]; + + float pixdim[8]; + std::memcpy(pixdim, buf + 76, 8 * sizeof(float)); + if (swap_endian) for(int i=0; i<8; i++) pixdim[i] = swap_float(pixdim[i]); + qfac = (pixdim[0] == 0.0f) ? 1.0f : pixdim[0]; + dx = pixdim[1]; dy = pixdim[2]; dz = pixdim[3]; + + int16_t sform16, qform16; + std::memcpy(&qform16, buf + 252, sizeof(int16_t)); + std::memcpy(&sform16, buf + 254, sizeof(int16_t)); + if (swap_endian) { qform16 = swap_int16(qform16); sform16 = swap_int16(sform16); } + qform_code = qform16; sform_code = sform16; + + std::memcpy(&b, buf + 256, sizeof(float)); + std::memcpy(&c, buf + 260, sizeof(float)); + std::memcpy(&d, buf + 264, sizeof(float)); + std::memcpy(&qoffset_x, buf + 268, sizeof(float)); + std::memcpy(&qoffset_y, buf + 272, sizeof(float)); + std::memcpy(&qoffset_z, buf + 276, sizeof(float)); + if (swap_endian) { + b = swap_float(b); c = swap_float(c); d = swap_float(d); + qoffset_x = swap_float(qoffset_x); qoffset_y = swap_float(qoffset_y); qoffset_z = swap_float(qoffset_z); + } + + std::memcpy(srow_x, buf + 280, 4 * sizeof(float)); + std::memcpy(srow_y, buf + 296, 4 * sizeof(float)); + std::memcpy(srow_z, buf + 312, 4 * sizeof(float)); + if (swap_endian) { + for(int i=0; i<4; i++) { + srow_x[i] = swap_float(srow_x[i]); + srow_y[i] = swap_float(srow_y[i]); + srow_z[i] = swap_float(srow_z[i]); + } + } + } else if (sizeof_hdr == 540) { // NIfTI-2 + if (f.gcount() < 540) { + std::cerr << "Error: Invalid NIfTI-2 file (too small)\n"; + return false; + } + int64_t dim[8]; + std::memcpy(dim, buf + 16, 8 * sizeof(int64_t)); + if (swap_endian) for(int i=0; i<8; i++) dim[i] = swap_int64(dim[i]); + dims[0] = static_cast(dim[1]); + dims[1] = static_cast(dim[2]); + dims[2] = static_cast(dim[3]); + + double pixdim[8]; + std::memcpy(pixdim, buf + 80, 8 * sizeof(double)); + if (swap_endian) for(int i=0; i<8; i++) pixdim[i] = swap_double(pixdim[i]); + qfac = (pixdim[0] == 0.0) ? 1.0f : static_cast(pixdim[0]); + dx = static_cast(pixdim[1]); + dy = static_cast(pixdim[2]); + dz = static_cast(pixdim[3]); + + int32_t sform32, qform32; + std::memcpy(&qform32, buf + 344, sizeof(int32_t)); + std::memcpy(&sform32, buf + 348, sizeof(int32_t)); + if (swap_endian) { qform32 = swap_int32(qform32); sform32 = swap_int32(sform32); } + qform_code = qform32; sform_code = sform32; + + double qb, qc, qd, qox, qoy, qoz; + std::memcpy(&qb, buf + 352, sizeof(double)); + std::memcpy(&qc, buf + 360, sizeof(double)); + std::memcpy(&qd, buf + 368, sizeof(double)); + std::memcpy(&qox, buf + 376, sizeof(double)); + std::memcpy(&qoy, buf + 384, sizeof(double)); + std::memcpy(&qoz, buf + 392, sizeof(double)); + if (swap_endian) { + qb = swap_double(qb); qc = swap_double(qc); qd = swap_double(qd); + qox = swap_double(qox); qoy = swap_double(qoy); qoz = swap_double(qoz); + } + b = static_cast(qb); c = static_cast(qc); d = static_cast(qd); + qoffset_x = static_cast(qox); qoffset_y = static_cast(qoy); qoffset_z = static_cast(qoz); + + double sx[4], sy[4], sz[4]; + std::memcpy(sx, buf + 400, 4 * sizeof(double)); + std::memcpy(sy, buf + 432, 4 * sizeof(double)); + std::memcpy(sz, buf + 464, 4 * sizeof(double)); + if (swap_endian) { + for(int i=0; i<4; i++) { + sx[i] = swap_double(sx[i]); + sy[i] = swap_double(sy[i]); + sz[i] = swap_double(sz[i]); + } + } + for(int i=0; i<4; i++) { + srow_x[i] = static_cast(sx[i]); + srow_y[i] = static_cast(sy[i]); + srow_z[i] = static_cast(sz[i]); + } + } else { + std::cerr << "Error: Unrecognized NIfTI file\n"; + return false; + } + + float v2r[4][4]; + if (sform_code > 0) { + for(int i=0; i<4; i++) { + v2r[0][i] = srow_x[i]; + v2r[1][i] = srow_y[i]; + v2r[2][i] = srow_z[i]; + } + v2r[3][0] = 0; v2r[3][1] = 0; v2r[3][2] = 0; v2r[3][3] = 1; + } else if (qform_code > 0) { + float b2 = b*b; + float c2 = c*c; + float d2 = d*d; + float a = std::sqrt(std::max(0.0f, 1.0f - b2 - c2 - d2)); + + float R[3][3]; + R[0][0] = a*a + b*b - c*c - d*d; + R[0][1] = 2.0f * (b*c - a*d); + R[0][2] = 2.0f * (b*d + a*c); + + R[1][0] = 2.0f * (b*c + a*d); + R[1][1] = a*a + c*c - b*b - d*d; + R[1][2] = 2.0f * (c*d - a*b); + + R[2][0] = 2.0f * (b*d - a*c); + R[2][1] = 2.0f * (c*d + a*b); + R[2][2] = a*a + d*d - c*c - b*b; + + v2r[0][0] = R[0][0] * dx; v2r[0][1] = R[0][1] * dy; v2r[0][2] = R[0][2] * qfac * dz; v2r[0][3] = qoffset_x; + v2r[1][0] = R[1][0] * dx; v2r[1][1] = R[1][1] * dy; v2r[1][2] = R[1][2] * qfac * dz; v2r[1][3] = qoffset_y; + v2r[2][0] = R[2][0] * dx; v2r[2][1] = R[2][1] * dy; v2r[2][2] = R[2][2] * qfac * dz; v2r[2][3] = qoffset_z; + v2r[3][0] = 0; v2r[3][1] = 0; v2r[3][2] = 0; v2r[3][3] = 1; + } else { + std::cerr << "Error: NIfTI file has no valid spatial transform\n"; + return false; + } + + out_header = json11::Json::object { + { "DIMENSIONS", json11::Json::array { dims[0], dims[1], dims[2] } }, + { "VOXEL_TO_RASMM", json11::Json::array { + json11::Json::array { v2r[0][0], v2r[0][1], v2r[0][2], v2r[0][3] }, + json11::Json::array { v2r[1][0], v2r[1][1], v2r[1][2], v2r[1][3] }, + json11::Json::array { v2r[2][0], v2r[2][1], v2r[2][2], v2r[2][3] }, + json11::Json::array { v2r[3][0], v2r[3][1], v2r[3][2], v2r[3][3] } + } } + }; + return true; +} + +bool save_trx(const Tractogram &tr, const std::string &out_path, const std::string &ref_nifti_path) { try { + json11::Json header_to_use = tr.header; + if (header_to_use.is_null() || !header_to_use["VOXEL_TO_RASMM"].is_array() || header_to_use["VOXEL_TO_RASMM"].array_items().empty()) { + if (ref_nifti_path.empty()) { + std::cerr << "Error: TCK/VTK -> TRX requires a reference NIfTI file\n"; + return false; + } + json11::Json ref_hdr; + if (!load_nifti_header(ref_nifti_path, ref_hdr)) return false; + auto obj = header_to_use.is_object() ? header_to_use.object_items() : std::map(); + obj["VOXEL_TO_RASMM"] = ref_hdr["VOXEL_TO_RASMM"]; + obj["DIMENSIONS"] = ref_hdr["DIMENSIONS"]; + header_to_use = obj; + } + if (tr.original_trx) { tr.original_trx->save(out_path, trx::TrxCompression::None); return true; @@ -326,7 +550,7 @@ bool save_trx(const Tractogram &tr, const std::string &out_path) { } // Copy header - trx.header = tr.header; + trx.header = header_to_use; trx.save(out_path, trx::TrxCompression::None); trx.close(); @@ -338,34 +562,76 @@ bool save_trx(const Tractogram &tr, const std::string &out_path) { } } -bool save_trk(const Tractogram &tr, const std::string &out_path, const std::string &original_filename) { +// Simple 4x4 matrix inversion helper for save_trk +bool invert_matrix4x4(const float m[4][4], float invOut[4][4]) { + float inv[16], det; + float m_1d[16]; + for(int i=0; i<4; i++) for(int j=0; j<4; j++) m_1d[i*4+j] = m[i][j]; + + inv[0] = m_1d[5] * m_1d[10] * m_1d[15] - m_1d[5] * m_1d[11] * m_1d[14] - m_1d[9] * m_1d[6] * m_1d[15] + m_1d[9] * m_1d[7] * m_1d[14] + m_1d[13] * m_1d[6] * m_1d[11] - m_1d[13] * m_1d[7] * m_1d[10]; + inv[4] = -m_1d[4] * m_1d[10] * m_1d[15] + m_1d[4] * m_1d[11] * m_1d[14] + m_1d[8] * m_1d[6] * m_1d[15] - m_1d[8] * m_1d[7] * m_1d[14] - m_1d[12] * m_1d[6] * m_1d[11] + m_1d[12] * m_1d[7] * m_1d[10]; + inv[8] = m_1d[4] * m_1d[9] * m_1d[15] - m_1d[4] * m_1d[11] * m_1d[13] - m_1d[8] * m_1d[5] * m_1d[15] + m_1d[8] * m_1d[7] * m_1d[13] + m_1d[12] * m_1d[5] * m_1d[11] - m_1d[12] * m_1d[7] * m_1d[9]; + inv[12] = -m_1d[4] * m_1d[9] * m_1d[14] + m_1d[4] * m_1d[10] * m_1d[13] + m_1d[8] * m_1d[5] * m_1d[14] - m_1d[8] * m_1d[6] * m_1d[13] - m_1d[12] * m_1d[5] * m_1d[10] + m_1d[12] * m_1d[6] * m_1d[9]; + inv[1] = -m_1d[1] * m_1d[10] * m_1d[15] + m_1d[1] * m_1d[11] * m_1d[14] + m_1d[9] * m_1d[2] * m_1d[15] - m_1d[9] * m_1d[3] * m_1d[14] - m_1d[13] * m_1d[2] * m_1d[11] + m_1d[13] * m_1d[3] * m_1d[10]; + inv[5] = m_1d[0] * m_1d[10] * m_1d[15] - m_1d[0] * m_1d[11] * m_1d[14] - m_1d[8] * m_1d[2] * m_1d[15] + m_1d[8] * m_1d[3] * m_1d[14] + m_1d[12] * m_1d[2] * m_1d[11] - m_1d[12] * m_1d[3] * m_1d[10]; + inv[9] = -m_1d[0] * m_1d[9] * m_1d[15] + m_1d[0] * m_1d[11] * m_1d[13] + m_1d[8] * m_1d[1] * m_1d[15] - m_1d[8] * m_1d[3] * m_1d[13] - m_1d[12] * m_1d[1] * m_1d[11] + m_1d[12] * m_1d[3] * m_1d[9]; + inv[13] = m_1d[0] * m_1d[9] * m_1d[14] - m_1d[0] * m_1d[10] * m_1d[13] - m_1d[8] * m_1d[1] * m_1d[14] + m_1d[8] * m_1d[2] * m_1d[13] + m_1d[12] * m_1d[1] * m_1d[10] - m_1d[12] * m_1d[2] * m_1d[9]; + inv[2] = m_1d[1] * m_1d[6] * m_1d[15] - m_1d[1] * m_1d[7] * m_1d[14] - m_1d[5] * m_1d[2] * m_1d[15] + m_1d[5] * m_1d[3] * m_1d[14] + m_1d[13] * m_1d[2] * m_1d[7] - m_1d[13] * m_1d[3] * m_1d[6]; + inv[6] = -m_1d[0] * m_1d[6] * m_1d[15] + m_1d[0] * m_1d[7] * m_1d[14] + m_1d[4] * m_1d[2] * m_1d[15] - m_1d[4] * m_1d[3] * m_1d[14] - m_1d[12] * m_1d[2] * m_1d[7] + m_1d[12] * m_1d[3] * m_1d[6]; + inv[10] = m_1d[0] * m_1d[5] * m_1d[15] - m_1d[0] * m_1d[7] * m_1d[13] - m_1d[4] * m_1d[1] * m_1d[15] + m_1d[4] * m_1d[3] * m_1d[13] + m_1d[12] * m_1d[1] * m_1d[7] - m_1d[12] * m_1d[3] * m_1d[5]; + inv[14] = -m_1d[0] * m_1d[5] * m_1d[14] + m_1d[0] * m_1d[6] * m_1d[13] + m_1d[4] * m_1d[1] * m_1d[14] - m_1d[4] * m_1d[2] * m_1d[13] - m_1d[12] * m_1d[1] * m_1d[6] + m_1d[12] * m_1d[2] * m_1d[5]; + inv[3] = -m_1d[1] * m_1d[6] * m_1d[11] + m_1d[1] * m_1d[7] * m_1d[10] + m_1d[5] * m_1d[2] * m_1d[11] - m_1d[5] * m_1d[3] * m_1d[10] - m_1d[9] * m_1d[2] * m_1d[7] + m_1d[9] * m_1d[3] * m_1d[6]; + inv[7] = m_1d[0] * m_1d[6] * m_1d[11] - m_1d[0] * m_1d[7] * m_1d[10] - m_1d[4] * m_1d[2] * m_1d[11] + m_1d[4] * m_1d[3] * m_1d[10] + m_1d[8] * m_1d[2] * m_1d[7] - m_1d[8] * m_1d[3] * m_1d[6]; + inv[11] = -m_1d[0] * m_1d[5] * m_1d[11] + m_1d[0] * m_1d[7] * m_1d[9] + m_1d[4] * m_1d[1] * m_1d[11] - m_1d[4] * m_1d[3] * m_1d[9] - m_1d[8] * m_1d[1] * m_1d[7] + m_1d[8] * m_1d[3] * m_1d[5]; + inv[15] = m_1d[0] * m_1d[5] * m_1d[10] - m_1d[0] * m_1d[6] * m_1d[9] - m_1d[4] * m_1d[1] * m_1d[10] + m_1d[4] * m_1d[2] * m_1d[9] + m_1d[8] * m_1d[1] * m_1d[6] - m_1d[8] * m_1d[2] * m_1d[5]; + + det = m_1d[0] * inv[0] + m_1d[1] * inv[4] + m_1d[2] * inv[8] + m_1d[3] * inv[12]; + if (det == 0) return false; + det = 1.0f / det; + for (int i = 0; i < 16; i++) { + invOut[i/4][i%4] = inv[i] * det; + } + return true; +} + +bool save_trk(const Tractogram &tr, const std::string &out_path, const std::string &original_filename, const std::string &ref_nifti_path) { std::ofstream f(out_path, std::ios::binary); if (!f.is_open()) return false; + json11::Json header_to_use = tr.header; + if (header_to_use.is_null() || !header_to_use["VOXEL_TO_RASMM"].is_array() || header_to_use["VOXEL_TO_RASMM"].array_items().empty()) { + if (ref_nifti_path.empty()) { + std::cerr << "Error: TCK/VTK -> TRK requires a reference NIfTI file\n"; + return false; + } + json11::Json ref_hdr; + if (!load_nifti_header(ref_nifti_path, ref_hdr)) return false; + auto obj = header_to_use.is_object() ? header_to_use.object_items() : std::map(); + obj["VOXEL_TO_RASMM"] = ref_hdr["VOXEL_TO_RASMM"]; + obj["DIMENSIONS"] = ref_hdr["DIMENSIONS"]; + header_to_use = obj; + } + TrkHeader header; std::memset(&header, 0, sizeof(header)); std::memcpy(header.magic_number, "TRACK", 5); - - // Default dimensions, voxel sizes and affine - header.dimensions[0] = 256; header.dimensions[1] = 256; header.dimensions[2] = 256; - header.voxel_sizes[0] = 1.0f; header.voxel_sizes[1] = 1.0f; header.voxel_sizes[2] = 1.0f; - for (int r = 0; r < 4; ++r) { - for (int c = 0; c < 4; ++c) { + + // Initialize vox_to_rasmm to identity + for (int r = 0; r < 4; ++r) + for (int c = 0; c < 4; ++c) header.voxel_to_rasmm[r][c] = (r == c) ? 1.0f : 0.0f; - } - } // Attempt to extract from JSON header - if (tr.header["DIMENSIONS"].is_array()) { - auto dims = tr.header["DIMENSIONS"].array_items(); + if (header_to_use["DIMENSIONS"].is_array()) { + auto dims = header_to_use["DIMENSIONS"].array_items(); if (dims.size() >= 3) { header.dimensions[0] = static_cast(dims[0].number_value()); header.dimensions[1] = static_cast(dims[1].number_value()); header.dimensions[2] = static_cast(dims[2].number_value()); } } - if (tr.header["VOXEL_TO_RASMM"].is_array()) { - auto rows = tr.header["VOXEL_TO_RASMM"].array_items(); + if (header_to_use["VOXEL_TO_RASMM"].is_array()) { + auto rows = header_to_use["VOXEL_TO_RASMM"].array_items(); if (rows.size() >= 4) { float vox_to_ras[4][4]; for (int r = 0; r < 4; ++r) { @@ -390,6 +656,16 @@ bool save_trk(const Tractogram &tr, const std::string &out_path, const std::stri f.write(reinterpret_cast(&header), 1000); + Eigen::Matrix4f mat = Eigen::Matrix4f::Identity(); + for (int r = 0; r < 4; ++r) + for (int c = 0; c < 4; ++c) + mat(r, c) = header.voxel_to_rasmm[r][c]; + Eigen::Matrix4f inv_mat = mat.inverse(); + + float vx = header.voxel_sizes[0] > 0 ? header.voxel_sizes[0] : 1.0f; + float vy = header.voxel_sizes[1] > 0 ? header.voxel_sizes[1] : 1.0f; + float vz = header.voxel_sizes[2] > 0 ? header.voxel_sizes[2] : 1.0f; + size_t num_streamlines = tr.offsets.size() - 1; std::vector chunk; chunk.reserve(4 * 1024 * 1024); @@ -399,15 +675,17 @@ bool save_trk(const Tractogram &tr, const std::string &out_path, const std::stri size_t end = tr.offsets[i+1]; int32_t n_pts = static_cast(end - start); - // Push n_pts const char* p_n_pts = reinterpret_cast(&n_pts); chunk.insert(chunk.end(), p_n_pts, p_n_pts + 4); - // Push points for (size_t j = start; j < end; ++j) { - float x = tr.pts[j*3]; - float y = tr.pts[j*3 + 1]; - float z = tr.pts[j*3 + 2]; + Eigen::Vector4f p_ras(tr.pts[j*3], tr.pts[j*3 + 1], tr.pts[j*3 + 2], 1.0f); + Eigen::Vector4f p_center = inv_mat * p_ras; + + float x = (p_center.x() + 0.5f) * vx; + float y = (p_center.y() + 0.5f) * vy; + float z = (p_center.z() + 0.5f) * vz; + const char* px = reinterpret_cast(&x); const char* py = reinterpret_cast(&y); const char* pz = reinterpret_cast(&z); From 8c1fa19c4ebb9cb365d6d339b14dc7053b3ca05f Mon Sep 17 00:00:00 2001 From: mattcieslak Date: Wed, 15 Jul 2026 15:22:39 -0400 Subject: [PATCH 04/11] Try to pass CI --- src/legacy_io.cpp | 131 +++++++++++++++++++++++----------- src/trx.cpp | 76 ++++++++++++-------- tests/test_trx_anytrxfile.cpp | 91 ++++++++++++++++++++++- 3 files changed, 229 insertions(+), 69 deletions(-) diff --git a/src/legacy_io.cpp b/src/legacy_io.cpp index d888c8c..4a9ec63 100644 --- a/src/legacy_io.cpp +++ b/src/legacy_io.cpp @@ -1,54 +1,65 @@ #include #include #include -#include -#include -#include -#include -#include -#include -#include #include -#include +#include +#include #include -#include -#include +#include +#include +#include +#include +#include +#include namespace trx { namespace legacy { +// Portable byte-swap helpers. These avoid the GCC/Clang-only __builtin_bswap* +// intrinsics so the file also compiles under MSVC; every modern compiler folds +// the shift/mask form back into a single bswap instruction. +inline uint16_t bswap16(uint16_t v) { + return static_cast((v << 8) | (v >> 8)); +} + +inline uint32_t bswap32(uint32_t v) { + return ((v & 0x000000FFu) << 24) | ((v & 0x0000FF00u) << 8) | + ((v & 0x00FF0000u) >> 8) | ((v & 0xFF000000u) >> 24); +} + +inline uint64_t bswap64(uint64_t v) { + return ((v & 0x00000000000000FFULL) << 56) | ((v & 0x000000000000FF00ULL) << 40) | + ((v & 0x0000000000FF0000ULL) << 24) | ((v & 0x00000000FF000000ULL) << 8) | + ((v & 0x000000FF00000000ULL) >> 8) | ((v & 0x0000FF0000000000ULL) >> 24) | + ((v & 0x00FF000000000000ULL) >> 40) | ((v & 0xFF00000000000000ULL) >> 56); +} + inline float swap_float(float f) { - union { - float f; - uint32_t i; - } u; - u.f = f; - u.i = __builtin_bswap32(u.i); - return u.f; + uint32_t i; + std::memcpy(&i, &f, sizeof(i)); + i = bswap32(i); + std::memcpy(&f, &i, sizeof(f)); + return f; } inline int32_t swap_int32(int32_t i) { - return __builtin_bswap32(i); + return static_cast(bswap32(static_cast(i))); } inline int16_t swap_int16(int16_t val) { - uint16_t uval = val; - uval = (uval << 8) | (uval >> 8); - return static_cast(uval); + return static_cast(bswap16(static_cast(val))); } inline int64_t swap_int64(int64_t val) { - return __builtin_bswap64(val); + return static_cast(bswap64(static_cast(val))); } inline double swap_double(double d) { - union { - double d; - uint64_t i; - } u; - u.d = d; - u.i = __builtin_bswap64(u.i); - return u.d; + uint64_t i; + std::memcpy(&i, &d, sizeof(i)); + i = bswap64(i); + std::memcpy(&d, &i, sizeof(d)); + return d; } @@ -138,13 +149,25 @@ bool load_trk(const std::string &filename, Tractogram &tr) { tr.offsets.push_back(0); tr.pts.clear(); + if (n_scalars < 0 || n_properties < 0) return false; + const size_t point_stride = (3u + static_cast(n_scalars)) * sizeof(float); + const size_t prop_bytes = static_cast(n_properties) * sizeof(float); + size_t offset = 1000; while (offset + sizeof(int32_t) <= buffer.size()) { int32_t n_points = *reinterpret_cast(buffer.data() + offset); offset += sizeof(int32_t); - + if (n_points < 0) return false; + + // Bounds-check the entire streamline record (points + trailing properties) + // before reading it, so a corrupt or oversized count can't drive an + // out-of-bounds read. buffer.size() - offset is safe here: the while + // condition guarantees offset <= buffer.size(). + const size_t bytes_needed = static_cast(n_points) * point_stride + prop_bytes; + if (bytes_needed > buffer.size() - offset) return false; + tr.offsets.push_back(tr.offsets.back() + n_points); - + for (int32_t j = 0; j < n_points; ++j) { float raw_x = *reinterpret_cast(buffer.data() + offset); float raw_y = *reinterpret_cast(buffer.data() + offset + 4); @@ -166,11 +189,11 @@ bool load_trk(const std::string &filename, Tractogram &tr) { tr.pts.push_back(y); tr.pts.push_back(z); - offset += (3 + n_scalars) * sizeof(float); + offset += point_stride; } - offset += n_properties * sizeof(float); + offset += prop_bytes; } - + return true; } @@ -189,9 +212,14 @@ bool load_tck(const std::string &filename, Tractogram &tr) { if (file_pos == std::string_view::npos) return false; size_t offset_pos = file_pos + 8; size_t offset_end = view.find_first_not_of("0123456789", offset_pos); - if (offset_end == std::string_view::npos) return false; - size_t offset = std::stoull(std::string(view.substr(offset_pos, offset_end - offset_pos))); - + if (offset_end == std::string_view::npos || offset_end == offset_pos) return false; + size_t offset; + try { + offset = std::stoull(std::string(view.substr(offset_pos, offset_end - offset_pos))); + } catch (const std::exception &) { + return false; // non-numeric or out-of-range data offset + } + if (offset >= buffer.size()) return false; const float* data = reinterpret_cast(buffer.data() + offset); @@ -249,7 +277,11 @@ bool load_vtk(const std::string &filename, Tractogram &tr) { while (std::getline(f, line)) { if (line.rfind("POINTS ", 0) == 0) { size_t space1 = line.find(" ", 7); - num_points = std::stoull(line.substr(7, space1 - 7)); + try { + num_points = std::stoull(line.substr(7, space1 - 7)); + } catch (const std::exception &) { + return false; // malformed POINTS count + } if (line.find("double", space1) != std::string::npos) { is_double = true; } @@ -258,6 +290,17 @@ bool load_vtk(const std::string &filename, Tractogram &tr) { } if (num_points == 0) return false; + const size_t elem_size = is_double ? sizeof(double) : sizeof(float); + if (num_points > std::numeric_limits::max() / (3 * elem_size)) return false; // overflow guard + + // Reject a point count that can't fit in the remaining file bytes, so a corrupt + // header can't trigger a huge allocation (and a truncated file fails cleanly). + const std::streampos data_start = f.tellg(); + f.seekg(0, std::ios::end); + const size_t bytes_available = static_cast(f.tellg() - data_start); + f.seekg(data_start); + if (num_points * 3 * elem_size > bytes_available) return false; + tr.pts.resize(num_points * 3); if (is_double) { std::vector dpts(num_points * 3); @@ -283,7 +326,11 @@ bool load_vtk(const std::string &filename, Tractogram &tr) { size_t num_streamlines = 0; while (std::getline(f, line)) { if (line.rfind("LINES ", 0) == 0) { - num_streamlines = std::stoull(line.substr(6, line.find(" ", 6) - 6)); + try { + num_streamlines = std::stoull(line.substr(6, line.find(" ", 6) - 6)); + } catch (const std::exception &) { + return false; // malformed LINES count + } break; } } @@ -531,6 +578,7 @@ bool save_trx(const Tractogram &tr, const std::string &out_path, const std::stri tr.original_trx->save(out_path, trx::TrxCompression::None); return true; } + if (tr.offsets.empty()) return false; // offsets must hold at least the trailing sentinel size_t nb_vertices = tr.pts.size() / 3; size_t nb_streamlines = tr.offsets.size() - 1; @@ -597,6 +645,7 @@ bool invert_matrix4x4(const float m[4][4], float invOut[4][4]) { bool save_trk(const Tractogram &tr, const std::string &out_path, const std::string &original_filename, const std::string &ref_nifti_path) { std::ofstream f(out_path, std::ios::binary); if (!f.is_open()) return false; + if (tr.offsets.empty()) return false; // offsets must hold at least the trailing sentinel json11::Json header_to_use = tr.header; if (header_to_use.is_null() || !header_to_use["VOXEL_TO_RASMM"].is_array() || header_to_use["VOXEL_TO_RASMM"].array_items().empty()) { @@ -710,9 +759,10 @@ bool save_trk(const Tractogram &tr, const std::string &out_path, const std::stri bool save_tck(const Tractogram &tr, const std::string &out_path) { std::ofstream f(out_path, std::ios::binary); if (!f.is_open()) return false; + if (tr.offsets.empty()) return false; // offsets must hold at least the trailing sentinel size_t num_streamlines = tr.offsets.size() - 1; - + // Build TCK header std::string header; size_t offset = 80; @@ -767,6 +817,7 @@ bool save_tck(const Tractogram &tr, const std::string &out_path) { bool save_vtk(const Tractogram &tr, const std::string &out_path) { std::ofstream f(out_path, std::ios::binary); if (!f.is_open()) return false; + if (tr.offsets.empty()) return false; // offsets must hold at least the trailing sentinel size_t num_streamlines = tr.offsets.size() - 1; size_t num_points = tr.pts.size() / 3; diff --git a/src/trx.cpp b/src/trx.cpp index 380ab89..cd88ed6 100644 --- a/src/trx.cpp +++ b/src/trx.cpp @@ -516,46 +516,66 @@ AnyTrxFile::_create_from_pointer(json header, auto arr = make_typed_array(elem_filename, static_cast(size), 1, ext); arr.materialize_to_owned(); trx.groups.emplace(base, std::move(arr)); - } else if (ext == "int64" || ext == "uint64" || ext == "int32" || ext == "uint8" || ext == "int8" || ext == "uint16" || ext == "int16") { - if (ext == "int32" || ext == "uint8" || ext == "int8" || ext == "uint16" || ext == "int16") { - std::cerr << "Warning: Upcasting group from " << ext << " to uint32\n"; - } - if (ext == "int64" || ext == "uint64") { - uint64_t num_strs = static_cast(header["NB_STREAMLINES"].number_value()); - if (num_strs > 4294967295ULL) { - throw TrxFormatError("downcasting is unsafe because the number of streamlines exceeds the 32-bit limit"); - } + } else if (ext == "int8" || ext == "uint8" || ext == "int16" || ext == "uint16" || ext == "int32" || + ext == "int64" || ext == "uint64") { + // Group index arrays should be uint32 per the TRX spec, but other integer + // dtypes are accepted for cross-language interoperability and normalized to + // uint32. The spec also requires every index to satisfy 0 <= id < + // NB_STREAMLINES, so each value is range-checked here: a negative or + // out-of-range index is rejected rather than silently wrapped into a + // valid-looking one. + // Copy into a plain local: capturing the structured binding `base` in the + // lambda below is only allowed from C++20 onward. + const std::string group_name = base; + const uint64_t nb_streamlines_u64 = static_cast(header["NB_STREAMLINES"].number_value()); + if (nb_streamlines_u64 > static_cast(std::numeric_limits::max())) { + throw TrxFormatError("Cannot normalize group '" + group_name + + "' to uint32: NB_STREAMLINES exceeds the uint32 limit"); } + auto tmp_arr = make_typed_array(elem_filename, static_cast(size), 1, ext); tmp_arr.materialize_to_owned(); + TypedArray arr; arr.dtype = "uint32"; arr.rows = static_cast(size); arr.cols = 1; arr.owned.resize(static_cast(size) * sizeof(uint32_t)); - uint32_t* dst = reinterpret_cast(arr.owned.data()); - if (ext == "int64") { - const int64_t* src = reinterpret_cast(tmp_arr.owned.data()); - for (size_t i = 0; i < size; ++i) dst[i] = static_cast(src[i]); - } else if (ext == "uint64") { - const uint64_t* src = reinterpret_cast(tmp_arr.owned.data()); - for (size_t i = 0; i < size; ++i) dst[i] = static_cast(src[i]); + uint32_t *dst = reinterpret_cast(arr.owned.data()); + + auto normalize = [&](auto src_tag) { + using S = decltype(src_tag); + const S *src = reinterpret_cast(tmp_arr.owned.data()); + for (long long i = 0; i < size; ++i) { + const S value = src[i]; + if constexpr (std::is_signed_v) { + if (value < 0) { + throw TrxFormatError("Group '" + group_name + "' contains a negative streamline index"); + } + } + if (static_cast(value) >= nb_streamlines_u64) { + throw TrxFormatError("Group '" + group_name + "' contains a streamline index >= NB_STREAMLINES"); + } + dst[i] = static_cast(value); + } + }; + + if (ext == "int8") { + normalize(int8_t{}); } else if (ext == "uint8") { - const uint8_t* src = reinterpret_cast(tmp_arr.owned.data()); - for (size_t i = 0; i < size; ++i) dst[i] = static_cast(src[i]); - } else if (ext == "int8") { - const int8_t* src = reinterpret_cast(tmp_arr.owned.data()); - for (size_t i = 0; i < size; ++i) dst[i] = static_cast(src[i]); - } else if (ext == "uint16") { - const uint16_t* src = reinterpret_cast(tmp_arr.owned.data()); - for (size_t i = 0; i < size; ++i) dst[i] = static_cast(src[i]); + normalize(uint8_t{}); } else if (ext == "int16") { - const int16_t* src = reinterpret_cast(tmp_arr.owned.data()); - for (size_t i = 0; i < size; ++i) dst[i] = static_cast(src[i]); + normalize(int16_t{}); + } else if (ext == "uint16") { + normalize(uint16_t{}); } else if (ext == "int32") { - const int32_t* src = reinterpret_cast(tmp_arr.owned.data()); - for (size_t i = 0; i < size; ++i) dst[i] = static_cast(src[i]); + normalize(int32_t{}); + } else if (ext == "int64") { + normalize(int64_t{}); + } else { // uint64 + normalize(uint64_t{}); } + trx.groups.emplace(base, std::move(arr)); } else { throw TrxDTypeError("Unsupported group dtype: " + ext); diff --git a/tests/test_trx_anytrxfile.cpp b/tests/test_trx_anytrxfile.cpp index 7dfe8e1..d375a3d 100644 --- a/tests/test_trx_anytrxfile.cpp +++ b/tests/test_trx_anytrxfile.cpp @@ -2,7 +2,9 @@ #include #include +#include #include +#include #include #include #include @@ -197,6 +199,21 @@ void write_zero_filled_file(const fs::path &file_path, const std::string &dtype, out.close(); } +// Write a group index file containing the exact `values`, laid out as the raw +// little-endian bytes of type T (the on-disk representation for a `.` +// group entry). Used to exercise the group dtype-normalization paths. +template void write_group_values(const fs::path &file_path, const std::vector &values) { + std::ofstream out(file_path.string(), std::ios::binary | std::ios::trunc); + if (!out.is_open()) { + throw std::runtime_error("Failed to write group file: " + file_path.string()); + } + if (!values.empty()) { + out.write(reinterpret_cast(values.data()), + static_cast(values.size() * sizeof(T))); + } + out.close(); +} + bool has_regular_file_recursive(const fs::path &dir_path) { std::error_code ec; for (fs::recursive_directory_iterator it(dir_path, ec), end; it != end; it.increment(ec)) { @@ -938,7 +955,9 @@ TEST(AnyTrxFile, UnsupportedGroupDtypeThrows) { write_zero_filled_file(group_file, "uint32", nb_streamlines); } const fs::path group_file = find_first_file_recursive(groups_dir); - rename_with_new_ext(group_file, "int32"); + // float32 is a valid TRX dtype but not a valid *group* dtype (groups must be an + // integer type), so loading must still reject it. + rename_with_new_ext(group_file, "float32"); EXPECT_THROW(load_any(corrupt_dir.string()), trx::TrxDTypeError); @@ -946,6 +965,76 @@ TEST(AnyTrxFile, UnsupportedGroupDtypeThrows) { fs::remove_all(temp_root, ec); } +// Groups stored in a non-uint32 integer dtype are accepted and normalized to +// uint32, preserving the index values (cross-language interoperability). +TEST(AnyTrxFile, GroupIntegerDtypeUpcastLoads) { + const auto gs_dir = require_gold_standard_dir(); + fs::path temp_root; + const fs::path dir = copy_gold_standard_dir(gs_dir, "trx_group_upcast", temp_root); + + const auto header = read_header_file(dir); + const auto nb_streamlines = static_cast(header["NB_STREAMLINES"].int_value()); + ASSERT_GE(nb_streamlines, 1u); + + // A handful of in-range indices, written as int32 (a non-uint32 integer dtype). + std::vector indices; + const uint32_t count = std::min(nb_streamlines, 4u); + for (uint32_t i = 0; i < count; ++i) { + indices.push_back(static_cast(i)); + } + + const fs::path groups_dir = dir / "groups"; + ensure_directory_exists(groups_dir); + write_group_values(groups_dir / "Upcast.int32", indices); + + auto loaded = load_any(dir.string()); + auto it = loaded.groups.find("Upcast"); + ASSERT_NE(it, loaded.groups.end()); + auto mat = it->second.as_matrix(); + ASSERT_EQ(static_cast(mat.size()), indices.size()); + for (size_t i = 0; i < indices.size(); ++i) { + EXPECT_EQ(mat(static_cast(i), 0), static_cast(indices[i])); + } + + std::error_code ec; + fs::remove_all(temp_root, ec); +} + +// A negative index in a signed group array must be rejected, not wrapped into a +// large positive uint32. +TEST(AnyTrxFile, GroupNegativeIndexThrows) { + const auto gs_dir = require_gold_standard_dir(); + fs::path temp_root; + const fs::path dir = copy_gold_standard_dir(gs_dir, "trx_group_negative", temp_root); + + const fs::path groups_dir = dir / "groups"; + ensure_directory_exists(groups_dir); + write_group_values(groups_dir / "Bad.int32", std::vector{0, -1}); + + EXPECT_THROW(load_any(dir.string()), trx::TrxFormatError); + + std::error_code ec; + fs::remove_all(temp_root, ec); +} + +// A 64-bit index above the uint32 range would silently wrap to a valid-looking +// index; it must be rejected instead. +TEST(AnyTrxFile, GroupIndexBeyondStreamlinesThrows) { + const auto gs_dir = require_gold_standard_dir(); + fs::path temp_root; + const fs::path dir = copy_gold_standard_dir(gs_dir, "trx_group_oob", temp_root); + + const fs::path groups_dir = dir / "groups"; + ensure_directory_exists(groups_dir); + // 5e9 > UINT32_MAX: static_cast would wrap it to ~705M. + write_group_values(groups_dir / "Bad.uint64", std::vector{0, 5000000000ULL}); + + EXPECT_THROW(load_any(dir.string()), trx::TrxFormatError); + + std::error_code ec; + fs::remove_all(temp_root, ec); +} + TEST(AnyTrxFile, InvalidEntryThrows) { const auto gs_dir = require_gold_standard_dir(); fs::path temp_root; From 87b090b99d3c6c8935d9692cce35f1d32d5ae723 Mon Sep 17 00:00:00 2001 From: frheault Date: Wed, 5 Aug 2026 13:44:38 -0400 Subject: [PATCH 05/11] fair zip mapping --- CMakeLists.txt | 10 + src/legacy_io.cpp | 12 +- src/trx.cpp | 640 ++++++++++++++++++++++++++-------- tests/CMakeLists.txt | 12 +- tests/test_trx_anytrxfile.cpp | 9 +- 5 files changed, 533 insertions(+), 150 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 9e7f30f..f1c1e63 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -221,6 +221,16 @@ if(TRX_BUILD_TESTS) if(NOT GTest_FOUND) find_package(GTest QUIET) endif() + if(NOT GTest_FOUND) + message(STATUS "GTest not found; fetching v1.14.0") + FetchContent_Declare( + googletest + GIT_REPOSITORY https://github.com/google/googletest.git + GIT_TAG v1.14.0 + ) + FetchContent_MakeAvailable(googletest) + set(GTest_FOUND TRUE) + endif() if(GTest_FOUND) enable_testing() add_subdirectory(tests) diff --git a/src/legacy_io.cpp b/src/legacy_io.cpp index 4a9ec63..2d2bb04 100644 --- a/src/legacy_io.cpp +++ b/src/legacy_io.cpp @@ -343,8 +343,8 @@ bool load_vtk(const std::string &filename, Tractogram &tr) { bool is_int64 = (line.find("int64") != std::string::npos); if (has_offsets) { - tr.offsets.resize(num_streamlines); - for (size_t i = 0; i < num_streamlines; ++i) { + tr.offsets.resize(num_streamlines + 1); + for (size_t i = 0; i <= num_streamlines; ++i) { if (is_int64) { uint64_t val; f.read(reinterpret_cast(&val), 8); @@ -366,6 +366,7 @@ bool load_vtk(const std::string &filename, Tractogram &tr) { tr.offsets.clear(); tr.offsets.push_back(0); + std::vector skip_buf; for (size_t i = 0; i < num_streamlines; ++i) { int32_t n_pts; @@ -375,8 +376,11 @@ bool load_vtk(const std::string &filename, Tractogram &tr) { if (n_pts == 0) continue; tr.offsets.push_back(tr.offsets.back() + n_pts); - // Skip cell indices - f.seekg(n_pts * sizeof(int32_t), std::ios::cur); + // Skip cell indices using read instead of seekg for performance + if (skip_buf.size() < static_cast(n_pts)) { + skip_buf.resize(n_pts); + } + f.read(reinterpret_cast(skip_buf.data()), n_pts * sizeof(int32_t)); } return true; diff --git a/src/trx.cpp b/src/trx.cpp index cd88ed6..b879d59 100644 --- a/src/trx.cpp +++ b/src/trx.cpp @@ -188,6 +188,61 @@ TypedArray make_typed_array(const std::string &filename, int rows, int cols, con } return array; } + +// Build a map of all ZIP local-file-header data offsets in a single linear +// pass over the archive. Calling find_uncompressed_zip_entry_offset once per +// entry (as the old code did) caused O(n * file_size) disk reads on files with +// many metadata arrays — catastrophic for a 6 GB archive. +using ZipOffsetMap = std::unordered_map>; + +ZipOffsetMap build_zip_offset_map(const std::string &zip_path) { + ZipOffsetMap result; + std::error_code ec; + const uintmax_t file_size = trx::fs::file_size(zip_path, ec); + if (ec || file_size < 30) { + return result; + } + + mio::shared_mmap_sink zip_mmap(zip_path, 0, file_size); + if (!zip_mmap.is_open() || zip_mmap.data() == nullptr) { + return result; + } + + const uint8_t *data = reinterpret_cast(zip_mmap.data()); + const size_t max_offset = file_size - 30; + + size_t curr = 0; + while (curr <= max_offset) { + if (data[curr] == 0x50 && data[curr + 1] == 0x4b && + data[curr + 2] == 0x03 && data[curr + 3] == 0x04) { + uint16_t name_len = static_cast(data[curr + 26]) | (static_cast(data[curr + 27]) << 8); + uint16_t extra_len = static_cast(data[curr + 28]) | (static_cast(data[curr + 29]) << 8); + uint32_t comp_size = static_cast(data[curr + 18]) | + (static_cast(data[curr + 19]) << 8) | + (static_cast(data[curr + 20]) << 16) | + (static_cast(data[curr + 21]) << 24); + uint32_t uncomp_size = static_cast(data[curr + 22]) | + (static_cast(data[curr + 23]) << 8) | + (static_cast(data[curr + 24]) << 16) | + (static_cast(data[curr + 25]) << 24); + + if (curr + 30 + name_len <= file_size) { + std::string cur_name(reinterpret_cast(data + curr + 30), name_len); + size_t payload_offset = curr + 30 + name_len + extra_len; + size_t payload_size = uncomp_size > 0 ? uncomp_size : comp_size; + if (payload_offset + payload_size <= file_size) { + result.emplace(normalize_slashes(cur_name), + std::make_pair(payload_offset, payload_size)); + } + } + curr += 30 + name_len + extra_len + comp_size; + } else { + curr++; + } + } + + return result; +} } // namespace std::string detect_positions_dtype(const std::string &path) { @@ -370,11 +425,215 @@ AnyTrxFile AnyTrxFile::load_from_zip(const std::string &filename) { throw TrxIOError("Could not open zip file: " + filename); } - std::string temp_dir = extract_zip_to_directory(zf.get()); + AnyTrxFile trx; + trx.header = load_header(zf.get()); + + if (!trx.header["NB_VERTICES"].is_number() || !trx.header["NB_STREAMLINES"].is_number()) { + throw TrxFormatError("Missing NB_VERTICES or NB_STREAMLINES in header.json"); + } + + const int nb_vertices = trx.header["NB_VERTICES"].int_value(); + const int nb_streamlines = trx.header["NB_STREAMLINES"].int_value(); + + // Build the offset map ONCE in a single pass over the archive, then reuse + // it for every entry. The previous approach called find_uncompressed_zip_entry_offset + // per entry, which re-scanned the entire file each time — O(n * file_size). + const ZipOffsetMap zip_offsets = build_zip_offset_map(filename); + + const zip_int64_t num_entries = zip_get_num_entries(zf.get(), 0); + for (zip_int64_t i = 0; i < num_entries; ++i) { + const char *raw_name = zip_get_name(zf.get(), i, 0); + if (raw_name == nullptr) { + continue; + } + std::string elem_filename(raw_name); + if (elem_filename.empty() || elem_filename.back() == '/') { + continue; + } + const std::string normalized = normalize_slashes(elem_filename); + if (normalized == "header.json") { + continue; + } + + std::string folder = folder_from_path(normalized, ""); + auto [base, dim, ext] = trx::detail::_split_ext_with_dimensionality(normalized); + ext = _normalize_dtype(ext); + + zip_stat_t st; + if (zip_stat_index(zf.get(), i, 0, &st) != 0) { + throw TrxIOError("Failed to stat zip entry: " + elem_filename); + } + long long raw_size_bytes = static_cast(st.size); + const int dtype_size = trx::detail::_sizeof_dtype(ext); + long long count_elems = (dtype_size > 0) ? (raw_size_bytes / dtype_size) : 0; + + auto read_entry_to_typed_array = [&](int rows, int cols) -> TypedArray { + TypedArray arr; + arr.dtype = ext; + arr.rows = rows; + arr.cols = cols; + + const size_t expected_bytes = static_cast(rows) * static_cast(cols) * static_cast(dtype_size); + + // If entry is stored uncompressed, map it directly from the ZIP file + // using the precomputed offset map (O(1) lookup, no per-entry rescan). + if (st.comp_method == ZIP_CM_STORE && expected_bytes > 0) { + auto it = zip_offsets.find(normalized); + if (it != zip_offsets.end()) { + const auto [offset, size] = it->second; + if (offset > 0 && size >= expected_bytes) { + arr.mmap = mio::shared_mmap_sink(filename, offset, expected_bytes); + if (arr.mmap.is_open() && arr.mmap.data() != nullptr) { + return arr; + } + } + } + } + + // Fallback: compressed entry or mmap failed — read via libzip. + detail::ZipFile entry_file(zip_fopen_index(zf.get(), i, 0)); + if (!entry_file) { + throw TrxIOError("Failed to open zip entry: " + elem_filename); + } + arr.owned.resize(expected_bytes); + if (expected_bytes > 0) { + zip_int64_t nbytes = zip_fread(entry_file.get(), arr.owned.data(), expected_bytes); + if (nbytes < static_cast(expected_bytes)) { + throw TrxIOError("Failed to read zip entry: " + elem_filename); + } + } + return arr; + }; + + if (base == "positions" && (folder.empty() || folder == ".")) { + if (count_elems != static_cast(nb_vertices) * 3 || dim != 3) { + throw TrxFormatError("Wrong positions size/dimensionality"); + } + if (ext != "float16" && ext != "float32" && ext != "float64") { + throw TrxDTypeError("Unsupported positions dtype: " + ext); + } + trx.positions = read_entry_to_typed_array(nb_vertices, 3); + } else if (base == "offsets" && (folder.empty() || folder == ".")) { + if (count_elems != static_cast(nb_streamlines) + 1 || dim != 1) { + throw TrxFormatError("Wrong offsets size/dimensionality"); + } + if (ext != "uint32" && ext != "uint64") { + throw TrxDTypeError("Unsupported offsets dtype: " + ext); + } + trx.offsets = read_entry_to_typed_array(nb_streamlines + 1, 1); + } else if (folder == "dps") { + const int nb_scalar = nb_streamlines > 0 ? static_cast(count_elems / nb_streamlines) : 0; + if (nb_streamlines == 0 || count_elems % nb_streamlines != 0 || nb_scalar != dim) { + throw TrxFormatError("Wrong dps size/dimensionality"); + } + trx.data_per_streamline.emplace(base, read_entry_to_typed_array(nb_streamlines, nb_scalar)); + } else if (folder == "dpv") { + const int nb_scalar = nb_vertices > 0 ? static_cast(count_elems / nb_vertices) : 0; + if (nb_vertices == 0 || count_elems % nb_vertices != 0 || nb_scalar != dim) { + throw TrxFormatError("Wrong dpv size/dimensionality"); + } + trx.data_per_vertex.emplace(base, read_entry_to_typed_array(nb_vertices, nb_scalar)); + } else if (folder.rfind("dpg", 0) == 0) { + if (count_elems != dim) { + throw TrxFormatError("Wrong dpg size/dimensionality"); + } + std::string data_name = path_basename(base); + std::string sub_folder = path_basename(folder); + trx.data_per_group[sub_folder].emplace(data_name, read_entry_to_typed_array(1, static_cast(count_elems))); + } else if (folder == "groups") { + if (dim != 1) { + throw TrxFormatError("Wrong group dimensionality"); + } + if (ext == "uint32") { + trx.groups.emplace(base, read_entry_to_typed_array(static_cast(count_elems), 1)); + } else if (ext == "int8" || ext == "uint8" || ext == "int16" || ext == "uint16" || ext == "int32" || + ext == "int64" || ext == "uint64") { + const std::string group_name = base; + const uint64_t nb_streamlines_u64 = static_cast(trx.header["NB_STREAMLINES"].number_value()); + if (nb_streamlines_u64 > static_cast(std::numeric_limits::max())) { + throw TrxFormatError("Cannot normalize group '" + group_name + "' to uint32: NB_STREAMLINES exceeds uint32 limit"); + } + auto tmp_arr = read_entry_to_typed_array(static_cast(count_elems), 1); + tmp_arr.materialize_to_owned(); + + TypedArray arr; + arr.dtype = "uint32"; + arr.rows = static_cast(count_elems); + arr.cols = 1; + arr.owned.resize(static_cast(count_elems) * sizeof(uint32_t)); + uint32_t *dst = reinterpret_cast(arr.owned.data()); + + auto normalize = [&](auto src_tag) { + using S = decltype(src_tag); + const S *src = reinterpret_cast(tmp_arr.owned.data()); + for (long long k = 0; k < count_elems; ++k) { + const S value = src[k]; + if constexpr (std::is_signed_v) { + if (value < 0) { + throw TrxFormatError("Group '" + group_name + "' contains a negative streamline index"); + } + } + if (static_cast(value) >= nb_streamlines_u64) { + throw TrxFormatError("Group '" + group_name + "' contains a streamline index >= NB_STREAMLINES"); + } + dst[k] = static_cast(value); + } + }; + + if (ext == "int8") normalize(int8_t{}); + else if (ext == "uint8") normalize(uint8_t{}); + else if (ext == "int16") normalize(int16_t{}); + else if (ext == "uint16") normalize(uint16_t{}); + else if (ext == "int32") normalize(int32_t{}); + else if (ext == "int64") normalize(int64_t{}); + else normalize(uint64_t{}); + + trx.groups.emplace(base, std::move(arr)); + } else { + throw TrxDTypeError("Unsupported group dtype: " + ext); + } + } else { + throw TrxFormatError("Entry is not part of a valid TRX structure: " + elem_filename); + } + } + + // Allow genuinely empty tractograms (NB_VERTICES=0, NB_STREAMLINES=0): they + // legitimately have no positions.* or offsets.* entries in the archive. + if ((trx.positions.empty() || trx.offsets.empty()) && + (nb_vertices > 0 || nb_streamlines > 0)) { + throw TrxFormatError("Missing essential data."); + } + + const size_t offsets_count = trx.offsets.size(); + if (offsets_count > 0) { + trx.offsets_u64.resize(offsets_count); + const auto bytes = trx.offsets.to_bytes(); + if (trx.offsets.dtype == "uint64") { + const auto *src = reinterpret_cast(bytes.data); + for (size_t k = 0; k < offsets_count; ++k) { + trx.offsets_u64[k] = src[k]; + } + } else if (trx.offsets.dtype == "uint32") { + const auto *src = reinterpret_cast(bytes.data); + for (size_t k = 0; k < offsets_count; ++k) { + trx.offsets_u64[k] = static_cast(src[k]); + } + } else { + throw TrxDTypeError("Unsupported offsets datatype: " + trx.offsets.dtype); + } + } + + if (offsets_count > 1) { + trx.lengths.resize(offsets_count - 1); + for (size_t k = 0; k + 1 < offsets_count; ++k) { + const uint64_t diff = trx.offsets_u64[k + 1] - trx.offsets_u64[k]; + if (diff > std::numeric_limits::max()) { + throw TrxFormatError("Offset difference exceeds uint32 range"); + } + trx.lengths[k] = static_cast(diff); + } + } - auto trx = AnyTrxFile::load_from_directory(temp_dir); - trx._uncompressed_folder_handle = temp_dir; - trx._owns_uncompressed_folder = true; return trx; } @@ -585,7 +844,10 @@ AnyTrxFile::_create_from_pointer(json header, } } - if (trx.positions.empty() || trx.offsets.empty()) { + // Allow genuinely empty tractograms (NB_VERTICES=0, NB_STREAMLINES=0): they + // legitimately have no positions.* or offsets.* files on disk. + if ((trx.positions.empty() || trx.offsets.empty()) && + (nb_vertices > 0 || nb_streamlines > 0)) { throw TrxFormatError("Missing essential data."); } @@ -622,44 +884,42 @@ AnyTrxFile::_create_from_pointer(json header, return trx; } -void write_positions_as_dtype(const AnyTrxFile &source, - TrxScalarType target_dtype, - const std::string &out_path, - size_t chunk_bytes) { - std::ofstream out(out_path, std::ios::binary | std::ios::trunc); - if (!out) - throw TrxIOError("Failed to create positions output: " + out_path); +std::vector convert_positions_to_vector(const AnyTrxFile &source, TrxScalarType target_dtype) { + const size_t total_points = source.num_vertices(); + const std::string target_dtype_str = scalar_type_name(target_dtype); + const size_t target_elem_size = static_cast(detail::_sizeof_dtype(target_dtype_str)); + const size_t total_bytes = total_points * 3 * target_elem_size; + std::vector out_buf(total_bytes); + + if (total_bytes == 0) { + return out_buf; + } + uint8_t *out_ptr = out_buf.data(); source.for_each_positions_chunk( - chunk_bytes, - [&](TrxScalarType src_dtype, const void *data, size_t /*point_offset*/, size_t point_count) { + 0 /* entire buffer in 1 chunk */, + [&](TrxScalarType src_dtype, const void *data, size_t point_offset, size_t point_count) { const size_t n = point_count * 3; + uint8_t *dst_chunk = out_ptr + point_offset * 3 * target_elem_size; - // Inner lambda: read from typed source pointer, cast to DstT, write to stream. auto write_as = [&](auto typed_src) { switch (target_dtype) { case TrxScalarType::Float16: { - std::vector buf(n); + auto *dst = reinterpret_cast(dst_chunk); for (size_t i = 0; i < n; ++i) - buf[i] = static_cast(static_cast(typed_src[i])); - out.write(reinterpret_cast(buf.data()), - static_cast(n * sizeof(Eigen::half))); + dst[i] = static_cast(static_cast(typed_src[i])); break; } case TrxScalarType::Float64: { - std::vector buf(n); + auto *dst = reinterpret_cast(dst_chunk); for (size_t i = 0; i < n; ++i) - buf[i] = static_cast(typed_src[i]); - out.write(reinterpret_cast(buf.data()), - static_cast(n * sizeof(double))); + dst[i] = static_cast(typed_src[i]); break; } default: { - std::vector buf(n); + auto *dst = reinterpret_cast(dst_chunk); for (size_t i = 0; i < n; ++i) - buf[i] = static_cast(typed_src[i]); - out.write(reinterpret_cast(buf.data()), - static_cast(n * sizeof(float))); + dst[i] = static_cast(typed_src[i]); break; } } @@ -678,10 +938,49 @@ void write_positions_as_dtype(const AnyTrxFile &source, } }); + return out_buf; +} + +void write_positions_as_dtype(const AnyTrxFile &source, + TrxScalarType target_dtype, + const std::string &out_path, + size_t chunk_bytes) { + static_cast(chunk_bytes); + std::ofstream out(out_path, std::ios::binary | std::ios::trunc); + if (!out) + throw TrxIOError("Failed to create positions output: " + out_path); + + std::vector buf = convert_positions_to_vector(source, target_dtype); + if (!buf.empty()) { + out.write(reinterpret_cast(buf.data()), static_cast(buf.size())); + } + if (out.bad()) throw TrxIOError("I/O error writing converted positions to: " + out_path); } +namespace { +std::string typed_array_filename(const std::string &base, const TypedArray &arr) { + if (arr.cols <= 1) { + return base + "." + arr.dtype; + } + return base + "." + std::to_string(arr.cols) + "." + arr.dtype; +} + +void write_typed_array_file(const std::string &path, const TypedArray &arr) { + const auto bytes = arr.to_bytes(); + std::ofstream out(path, std::ios::binary | std::ios::out | std::ios::trunc); + if (!out.is_open()) { + throw TrxIOError("Failed to open output file: " + path); + } + if (bytes.data && bytes.size > 0) { + out.write(reinterpret_cast(bytes.data), static_cast(bytes.size)); + } + out.flush(); + out.close(); +} +} // namespace + void AnyTrxFile::save(const std::string &filename, TrxCompression compression) { TrxSaveOptions options; options.compression = compression; @@ -698,95 +997,144 @@ void AnyTrxFile::save(const std::string &filename, const TrxSaveOptions &options throw TrxDTypeError("Unsupported extension: " + ext); } - if (offsets.empty()) { - throw TrxFormatError("Cannot save TRX without offsets data"); - } - if (offsets_u64.empty()) { - throw TrxFormatError("Cannot save TRX without decoded offsets"); - } - if (header["NB_STREAMLINES"].is_number()) { - const auto nb_streamlines = static_cast(header["NB_STREAMLINES"].int_value()); - if (offsets_u64.size() != nb_streamlines + 1) { - throw TrxFormatError("TRX offsets size does not match NB_STREAMLINES"); + const bool is_empty_tractogram = + header["NB_VERTICES"].is_number() && header["NB_STREAMLINES"].is_number() && + header["NB_VERTICES"].int_value() == 0 && header["NB_STREAMLINES"].int_value() == 0; + + if (!is_empty_tractogram) { + if (offsets.empty()) { + throw TrxFormatError("Cannot save TRX without offsets data"); } - } - if (header["NB_VERTICES"].is_number()) { - const auto nb_vertices = static_cast(header["NB_VERTICES"].int_value()); - const auto last = offsets_u64.back(); - if (last != nb_vertices) { - throw TrxFormatError("TRX offsets sentinel does not match NB_VERTICES"); + if (offsets_u64.empty()) { + throw TrxFormatError("Cannot save TRX without decoded offsets"); } - } - for (size_t i = 1; i < offsets_u64.size(); ++i) { - if (offsets_u64[i] < offsets_u64[i - 1]) { - throw TrxFormatError("TRX offsets must be monotonically increasing"); + if (header["NB_STREAMLINES"].is_number()) { + const auto nb_streamlines = static_cast(header["NB_STREAMLINES"].int_value()); + if (offsets_u64.size() != nb_streamlines + 1) { + throw TrxFormatError("TRX offsets size does not match NB_STREAMLINES"); + } } - } - if (!positions.empty()) { - const auto last = offsets_u64.back(); - if (last != static_cast(positions.rows)) { - throw TrxFormatError("TRX positions row count does not match offsets sentinel"); + if (header["NB_VERTICES"].is_number()) { + const auto nb_vertices = static_cast(header["NB_VERTICES"].int_value()); + const auto last = offsets_u64.back(); + if (last != nb_vertices) { + throw TrxFormatError("TRX offsets sentinel does not match NB_VERTICES"); + } + } + for (size_t i = 1; i < offsets_u64.size(); ++i) { + if (offsets_u64[i] < offsets_u64[i - 1]) { + throw TrxFormatError("TRX offsets must be monotonically increasing"); + } + } + if (!positions.empty()) { + const auto last = offsets_u64.back(); + if (last != static_cast(positions.rows)) { + throw TrxFormatError("TRX positions row count does not match offsets sentinel"); + } } - } - - const std::string source_dir = - !_uncompressed_folder_handle.empty() ? _uncompressed_folder_handle : _backing_directory; - if (source_dir.empty()) { - throw TrxIOError("TRX file has no backing directory to save from"); } if (save_mode == TrxSaveMode::Archive) { - int errorp; + int errorp = 0; detail::ZipArchive zf(zip_open(filename.c_str(), ZIP_CREATE + ZIP_TRUNCATE, &errorp)); if (!zf) { throw TrxIOError("Could not open archive " + filename + ": " + strerror(errorp)); } + const zip_int32_t compression = static_cast(to_zip_compression(options.compression)); + + auto add_zip_buffer_entry = [&](const std::string &entry_name, const void *data, size_t size) { + void *buf = std::malloc(size > 0 ? size : 1); + if (!buf) { + throw TrxIOError("Failed to allocate buffer for zip entry: " + entry_name); + } + if (size > 0 && data != nullptr) { + std::memcpy(buf, data, size); + } + zip_source_t *src = zip_source_buffer(zf.get(), buf, size, 1 /* freep=1 */); + if (!src) { + std::free(buf); + throw TrxIOError("zip_source_buffer failed for: " + entry_name); + } + const zip_int64_t idx = zip_file_add(zf.get(), entry_name.c_str(), src, ZIP_FL_ENC_UTF_8 | ZIP_FL_OVERWRITE); + if (idx < 0) { + throw TrxIOError("Failed to add entry to zip: " + entry_name + ": " + std::string(zip_strerror(zf.get()))); + } + if (zip_set_file_compression(zf.get(), idx, compression, 0) < 0) { + throw TrxIOError("Failed to set compression for zip entry: " + entry_name + ": " + std::string(zip_strerror(zf.get()))); + } + }; + + // 1. Header const std::string header_payload = header.dump() + "\n"; - zip_source_t *header_source = - zip_source_buffer(zf.get(), header_payload.data(), header_payload.size(), 0 /* do not free */); - if (header_source == nullptr) { - throw TrxIOError("Failed to create zip source for header.json: " + - std::string(zip_strerror(zf.get()))); + add_zip_buffer_entry("header.json", header_payload.data(), header_payload.size()); + + // 2. Positions + if (options.positions_dtype.has_value() && !positions.empty()) { + const TrxScalarType target = *options.positions_dtype; + std::vector converted_pos = convert_positions_to_vector(*this, target); + const std::string pos_name = "positions.3." + scalar_type_name(target); + add_zip_buffer_entry(pos_name, converted_pos.data(), converted_pos.size()); + } else if (!positions.empty()) { + auto pos_bytes = positions.to_bytes(); + const std::string pos_name = "positions.3." + positions.dtype; + add_zip_buffer_entry(pos_name, pos_bytes.data, pos_bytes.size); + } + + // 3. Offsets + if (!offsets.empty()) { + auto off_bytes = offsets.to_bytes(); + const std::string off_name = "offsets." + offsets.dtype; + add_zip_buffer_entry(off_name, off_bytes.data, off_bytes.size); + } + + // 4. Groups + if (!groups.empty()) { + zip_dir_add(zf.get(), "groups", ZIP_FL_ENC_UTF_8); + for (const auto &kv : groups) { + const std::string entry_name = "groups/" + typed_array_filename(kv.first, kv.second); + auto bytes = kv.second.to_bytes(); + add_zip_buffer_entry(entry_name, bytes.data, bytes.size); + } } - const zip_int64_t header_idx = - zip_file_add(zf.get(), "header.json", header_source, ZIP_FL_ENC_UTF_8 | ZIP_FL_OVERWRITE); - if (header_idx < 0) { - throw TrxIOError("Failed to add header.json to archive: " + std::string(zip_strerror(zf.get()))); + + // 5. DPS + if (!data_per_streamline.empty()) { + zip_dir_add(zf.get(), "dps", ZIP_FL_ENC_UTF_8); + for (const auto &kv : data_per_streamline) { + const std::string entry_name = "dps/" + typed_array_filename(kv.first, kv.second); + auto bytes = kv.second.to_bytes(); + add_zip_buffer_entry(entry_name, bytes.data, bytes.size); + } } - const zip_int32_t compression = static_cast(to_zip_compression(options.compression)); - if (zip_set_file_compression(zf.get(), header_idx, compression, 0) < 0) { - throw TrxIOError("Failed to set compression for header.json: " + - std::string(zip_strerror(zf.get()))); + + // 6. DPV + if (!data_per_vertex.empty()) { + zip_dir_add(zf.get(), "dpv", ZIP_FL_ENC_UTF_8); + for (const auto &kv : data_per_vertex) { + const std::string entry_name = "dpv/" + typed_array_filename(kv.first, kv.second); + auto bytes = kv.second.to_bytes(); + add_zip_buffer_entry(entry_name, bytes.data, bytes.size); + } + } + + // 7. DPG + if (!data_per_group.empty()) { + zip_dir_add(zf.get(), "dpg", ZIP_FL_ENC_UTF_8); + for (const auto &group_kv : data_per_group) { + const std::string sub_dir = "dpg/" + group_kv.first; + zip_dir_add(zf.get(), sub_dir.c_str(), ZIP_FL_ENC_UTF_8); + for (const auto &kv : group_kv.second) { + const std::string entry_name = sub_dir + "/" + typed_array_filename(kv.first, kv.second); + auto bytes = kv.second.to_bytes(); + add_zip_buffer_entry(entry_name, bytes.data, bytes.size); + } + } } - std::unordered_set skip = {"header.json"}; - // Guard deletes the temp positions file after commit (or on exception). - TempFileGuard tmp_pos_guard; - if (options.positions_dtype.has_value() && !positions.empty()) { - const TrxScalarType target = *options.positions_dtype; - const std::string new_dtype_str = scalar_type_name(target); - if (new_dtype_str != positions.dtype) { - skip.insert("positions.3." + positions.dtype); - tmp_pos_guard.path = make_unique_temp_path("trx_pos_convert"); - write_positions_as_dtype(*this, target, tmp_pos_guard.path); - const std::string new_pos_name = "positions.3." + new_dtype_str; - zip_source_t *pos_src = - zip_source_file(zf.get(), tmp_pos_guard.path.c_str(), 0, -1); - if (!pos_src) - throw TrxIOError("Failed to create zip source for converted positions"); - const zip_int64_t pos_idx = - zip_file_add(zf.get(), new_pos_name.c_str(), pos_src, ZIP_FL_ENC_UTF_8 | ZIP_FL_OVERWRITE); - if (pos_idx < 0) - throw TrxIOError("Failed to add converted positions to archive: " + - std::string(zip_strerror(zf.get()))); - if (zip_set_file_compression(zf.get(), pos_idx, compression, 0) < 0) - throw TrxIOError("Failed to set compression for converted positions"); - } - } - zip_from_folder(zf.get(), source_dir, source_dir, to_zip_compression(options.compression), &skip); zf.commit(filename); } else { + // TrxSaveMode::Directory std::error_code ec; if (trx::fs::exists(filename, ec) && trx::fs::is_directory(filename, ec)) { if (!options.overwrite_existing) { @@ -800,42 +1148,72 @@ void AnyTrxFile::save(const std::string &filename, const TrxSaveOptions &options if (dest_path.has_parent_path()) { mkdir_or_throw(dest_path.parent_path().string()); } - std::error_code source_ec; - const trx::fs::path source_path = trx::fs::weakly_canonical(trx::fs::path(source_dir), source_ec); - std::error_code dest_ec; - const trx::fs::path normalized_dest = trx::fs::weakly_canonical(dest_path, dest_ec); - const bool same_directory = !source_ec && !dest_ec && source_path == normalized_dest; + mkdir_or_throw(filename); - if (!same_directory) { - copy_dir(source_dir, filename); + const trx::fs::path final_header_path = dest_path / "header.json"; + std::ofstream out_json(final_header_path, std::ios::out | std::ios::trunc); + if (!out_json.is_open()) { + throw TrxIOError("Failed to write header.json to: " + final_header_path.string()); } + out_json << header.dump() << std::endl; + out_json.close(); if (options.positions_dtype.has_value() && !positions.empty()) { const TrxScalarType target = *options.positions_dtype; const std::string new_dtype_str = scalar_type_name(target); - if (new_dtype_str != positions.dtype) { - const std::string old_pos = filename + SEPARATOR + "positions.3." + positions.dtype; - const std::string new_pos = filename + SEPARATOR + "positions.3." + new_dtype_str; - write_positions_as_dtype(*this, target, new_pos); - std::error_code rm_ec; - trx::fs::remove(old_pos, rm_ec); + const std::string new_pos = filename + SEPARATOR + "positions.3." + new_dtype_str; + auto converted_pos = convert_positions_to_vector(*this, target); + std::ofstream out_pos(new_pos, std::ios::binary | std::ios::out | std::ios::trunc); + if (!out_pos.is_open()) { + throw TrxIOError("Failed to write positions to: " + new_pos); } + if (!converted_pos.empty()) { + out_pos.write(reinterpret_cast(converted_pos.data()), converted_pos.size()); + } + } else if (!positions.empty()) { + const std::string pos_path = filename + SEPARATOR + "positions.3." + positions.dtype; + write_typed_array_file(pos_path, positions); } - const trx::fs::path final_header_path = dest_path / "header.json"; - std::ofstream out_json(final_header_path, std::ios::out | std::ios::trunc); - if (!out_json.is_open()) { - throw TrxIOError("Failed to write header.json to: " + final_header_path.string()); + if (!offsets.empty()) { + const std::string off_path = filename + SEPARATOR + typed_array_filename("offsets", offsets); + write_typed_array_file(off_path, offsets); } - out_json << header.dump() << std::endl; - out_json.close(); - ec.clear(); - if (!trx::fs::exists(filename, ec) || !trx::fs::is_directory(filename, ec)) { - throw TrxIOError("Failed to create output directory: " + filename); + if (!groups.empty()) { + const std::string groups_dir = filename + SEPARATOR + "groups"; + trx::fs::create_directories(groups_dir, ec); + for (const auto &kv : groups) { + write_typed_array_file(groups_dir + SEPARATOR + typed_array_filename(kv.first, kv.second), kv.second); + } + } + + if (!data_per_streamline.empty()) { + const std::string dps_dir = filename + SEPARATOR + "dps"; + trx::fs::create_directories(dps_dir, ec); + for (const auto &kv : data_per_streamline) { + write_typed_array_file(dps_dir + SEPARATOR + typed_array_filename(kv.first, kv.second), kv.second); + } + } + + if (!data_per_vertex.empty()) { + const std::string dpv_dir = filename + SEPARATOR + "dpv"; + trx::fs::create_directories(dpv_dir, ec); + for (const auto &kv : data_per_vertex) { + write_typed_array_file(dpv_dir + SEPARATOR + typed_array_filename(kv.first, kv.second), kv.second); + } } - if (!trx::fs::exists(final_header_path)) { - throw TrxFormatError("Missing header.json in output directory: " + final_header_path.string()); + + if (!data_per_group.empty()) { + const std::string dpg_dir = filename + SEPARATOR + "dpg"; + trx::fs::create_directories(dpg_dir, ec); + for (const auto &group_kv : data_per_group) { + const std::string group_dir = dpg_dir + SEPARATOR + group_kv.first; + trx::fs::create_directories(group_dir, ec); + for (const auto &kv : group_kv.second) { + write_typed_array_file(group_dir + SEPARATOR + typed_array_filename(kv.first, kv.second), kv.second); + } + } } } } @@ -1308,26 +1686,6 @@ TrxScalarType scalar_type_from_dtype(const std::string &dtype) { } return TrxScalarType::Float32; } - -std::string typed_array_filename(const std::string &base, const TypedArray &arr) { - if (arr.cols <= 1) { - return base + "." + arr.dtype; - } - return base + "." + std::to_string(arr.cols) + "." + arr.dtype; -} - -void write_typed_array_file(const std::string &path, const TypedArray &arr) { - const auto bytes = arr.to_bytes(); - std::ofstream out(path, std::ios::binary | std::ios::out | std::ios::trunc); - if (!out.is_open()) { - throw TrxIOError("Failed to open output file: " + path); - } - if (bytes.data && bytes.size > 0) { - out.write(reinterpret_cast(bytes.data), static_cast(bytes.size)); - } - out.flush(); - out.close(); -} } // namespace void AnyTrxFile::for_each_positions_chunk(size_t chunk_bytes, const PositionsChunkCallback &fn) const { diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index f37f3db..bb00ab4 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -1,8 +1,14 @@ enable_testing() -find_package(GTest CONFIG QUIET) -if(NOT GTest_FOUND) - find_package(GTest REQUIRED) +if(NOT TARGET GTest::gtest_main AND NOT TARGET gtest_main) + find_package(GTest CONFIG QUIET) + if(NOT GTest_FOUND) + find_package(GTest QUIET) + endif() +endif() +if(TARGET gtest_main AND NOT TARGET GTest::gtest_main) + add_library(GTest::gtest_main ALIAS gtest_main) + add_library(GTest::gtest ALIAS gtest) endif() set(TRX_TEST_DATA_REPO "https://github.com/tee-ar-ex/trx-test-data" CACHE STRING "Git repository URL for test data.") diff --git a/tests/test_trx_anytrxfile.cpp b/tests/test_trx_anytrxfile.cpp index d375a3d..01de042 100644 --- a/tests/test_trx_anytrxfile.cpp +++ b/tests/test_trx_anytrxfile.cpp @@ -1218,7 +1218,7 @@ TEST(AnyTrxFile, SaveRejectsPositionsRowMismatch) { fs::remove_all(temp_dir, ec); } -TEST(AnyTrxFile, SaveRejectsMissingBackingDirectory) { +TEST(AnyTrxFile, SaveWithoutBackingDirectorySucceeds) { const auto gs_dir = require_gold_standard_dir(); const fs::path gs_trx = gs_dir / "gs_fldr.trx"; auto trx = load_any(gs_trx.string()); @@ -1228,7 +1228,12 @@ TEST(AnyTrxFile, SaveRejectsMissingBackingDirectory) { const auto temp_dir = make_temp_test_dir("trx_any_save_no_backing"); const fs::path out_path = temp_dir / "no_backing.trx"; - EXPECT_THROW(trx.save(out_path.string(), trx::TrxCompression::None), trx::TrxIOError); + EXPECT_NO_THROW(trx.save(out_path.string(), trx::TrxCompression::None)); + + auto loaded = load_any(out_path.string()); + EXPECT_EQ(loaded.num_streamlines(), trx.num_streamlines()); + EXPECT_EQ(loaded.num_vertices(), trx.num_vertices()); + loaded.close(); trx.close(); std::error_code ec; From 53150f91b20474fddfb8a996c294683f9c35e424 Mon Sep 17 00:00:00 2001 From: frheault Date: Wed, 5 Aug 2026 14:25:10 -0400 Subject: [PATCH 06/11] Update for subset tests --- examples/trxinfo.cpp | 2 +- include/trx/trx.h | 5 ++--- include/trx/trx.tpp | 22 ++++++++++++---------- src/trx.cpp | 2 +- tests/test_trx_mmap.cpp | 2 +- 5 files changed, 17 insertions(+), 16 deletions(-) diff --git a/examples/trxinfo.cpp b/examples/trxinfo.cpp index abec187..9617603 100644 --- a/examples/trxinfo.cpp +++ b/examples/trxinfo.cpp @@ -152,7 +152,7 @@ void print_trx_info(const trx::AnyTrxFile &trx, const std::string &path, bool is } } } else { - std::cout << " " << colorize(colors, colors.cyan, "Data per group") << ": none\n"; + std::cout << " " << trx_cli::colorize(colors, colors.cyan, "Data per group") << ": none\n"; } } diff --git a/include/trx/trx.h b/include/trx/trx.h index 97babf6..757e128 100644 --- a/include/trx/trx.h +++ b/include/trx/trx.h @@ -40,9 +40,8 @@ namespace trx { namespace fs = std::filesystem; -} - using json = json11::Json; +} namespace trx { enum class TrxSaveMode { Auto, Archive, Directory }; @@ -297,7 +296,7 @@ template class TrxFile { std::unique_ptr> deepcopy(); /** - * @brief Remove the ununsed portion of preallocated memmaps + * @brief Remove the unused portion of preallocated memmaps * * @param nb_streamlines The number of streamlines to keep * @param nb_vertices The number of vertices to keep diff --git a/include/trx/trx.tpp b/include/trx/trx.tpp index ac92ef7..296eb39 100644 --- a/include/trx/trx.tpp +++ b/include/trx/trx.tpp @@ -144,8 +144,10 @@ void copy_cast_from_dtype_buffer(const void *src, size_t n, const std::string &d template void write_binary(const std::string &filename, const Matrix &matrix) { std::ofstream out(filename, std::ios::out | std::ios::binary | std::ios::trunc); typename Matrix::Index rows = matrix.rows(), cols = matrix.cols(); - // out.write((char *)(&rows), sizeof(typename Matrix::Index)); - // out.write((char *)(&cols), sizeof(typename Matrix::Index)); + auto *rows_ptr = reinterpret_cast(&rows); // check_syntax off + auto *cols_ptr = reinterpret_cast(&cols); // check_syntax off + out.write(rows_ptr, sizeof(typename Matrix::Index)); + out.write(cols_ptr, sizeof(typename Matrix::Index)); const auto *data = reinterpret_cast(matrix.data()); // check_syntax off out.write(data, rows * cols * sizeof(typename Matrix::Scalar)); out.close(); @@ -244,7 +246,7 @@ TrxFile
::TrxFile(int nb_vertices, int nb_streamlines, const TrxFile
*ini throw TrxArgumentError("Can't use init_as without declaring nb_vertices and nb_streamlines"); } - // will remove as completely unecessary. using as placeholders + // will remove as completely unnecessary. using as placeholders this->header = {}; this->streamlines.reset(); @@ -432,7 +434,7 @@ TrxFile
::_create_trx_from_pointer(json header, auto [base, dim, ext] = trx::detail::_split_ext_with_dimensionality(elem_filename); - long long mem_adress = std::get<0>(x->second); + long long mem_address = std::get<0>(x->second); long long size = std::get<1>(x->second); if (base == "positions" && (folder.empty() || folder == ".")) { @@ -446,7 +448,7 @@ TrxFile
::_create_trx_from_pointer(json header, std::tuple shape = std::make_tuple(static_cast(trx->header["NB_VERTICES"].int_value()), 3); trx->streamlines->mmap_pos = - trx::_create_memmap(filename, shape, "r+", ext, mem_adress); + trx::_create_memmap(filename, shape, "r+", ext, mem_address); trx::detail::remap(trx->streamlines->_data, trx->streamlines->mmap_pos.data(), shape); } @@ -466,7 +468,7 @@ TrxFile
::_create_trx_from_pointer(json header, const int offsets_rows = missing_sentinel ? (nb_str + 1) : static_cast(size); std::tuple shape = std::make_tuple(offsets_rows, 1); trx->streamlines->mmap_off = trx::_create_memmap(filename, std::make_tuple(static_cast(size), 1), "r+", - ext, mem_adress); + ext, mem_address); if (ext == "uint64") { if (missing_sentinel) { @@ -509,7 +511,7 @@ TrxFile
::_create_trx_from_pointer(json header, } else { shape = std::make_tuple(static_cast(trx->header["NB_STREAMLINES"].int_value()), nb_scalar); } - trx->data_per_streamline[base]->mmap = trx::_create_memmap(filename, shape, "r+", ext, mem_adress); + trx->data_per_streamline[base]->mmap = trx::_create_memmap(filename, shape, "r+", ext, mem_address); const std::string expected_dtype = dtype_from_scalar
(); if (ext == expected_dtype) { trx::detail::remap(trx->data_per_streamline[base]->_matrix, trx->data_per_streamline[base]->mmap.data(), shape); @@ -534,7 +536,7 @@ TrxFile
::_create_trx_from_pointer(json header, } else { shape = std::make_tuple(static_cast(trx->header["NB_VERTICES"].int_value()), nb_scalar); } - trx->data_per_vertex[base]->mmap_pos = trx::_create_memmap(filename, shape, "r+", ext, mem_adress); + trx->data_per_vertex[base]->mmap_pos = trx::_create_memmap(filename, shape, "r+", ext, mem_address); const std::string expected_dtype = dtype_from_scalar
(); if (ext == expected_dtype) { trx::detail::remap(trx->data_per_vertex[base]->_data, trx->data_per_vertex[base]->mmap_pos.data(), shape); @@ -564,7 +566,7 @@ TrxFile
::_create_trx_from_pointer(json header, std::string sub_folder = path_basename(folder); trx->data_per_group[sub_folder][data_name] = std::make_unique>(); - trx->data_per_group[sub_folder][data_name]->mmap = trx::_create_memmap(filename, shape, "r+", ext, mem_adress); + trx->data_per_group[sub_folder][data_name]->mmap = trx::_create_memmap(filename, shape, "r+", ext, mem_address); const std::string expected_dtype = dtype_from_scalar
(); if (ext == expected_dtype) { trx::detail::remap(trx->data_per_group[sub_folder][data_name]->_matrix, @@ -593,7 +595,7 @@ TrxFile
::_create_trx_from_pointer(json header, info.rows = std::get<0>(shape); info.cols = std::get<1>(shape); info.dtype = ext; - info.mem_offset = mem_adress; + info.mem_offset = mem_address; trx->group_backing_info_[base] = std::move(info); } else { throw TrxFormatError("Entry is not part of a valid TRX structure: " + elem_filename); diff --git a/src/trx.cpp b/src/trx.cpp index b879d59..6482070 100644 --- a/src/trx.cpp +++ b/src/trx.cpp @@ -1336,7 +1336,7 @@ void allocate_file(const std::string &path, std::size_t size) { file.flush(); file.close(); } else { - std::cerr << "Failed to allocate file : " << sys_error() << '\n'; + std::cerr << "Failed to allocate file: " << sys_error() << '\n'; } } diff --git a/tests/test_trx_mmap.cpp b/tests/test_trx_mmap.cpp index 376dc8e..ff0f601 100644 --- a/tests/test_trx_mmap.cpp +++ b/tests/test_trx_mmap.cpp @@ -21,7 +21,7 @@ json load_header(zip_t *zfolder); } // namespace trx using namespace Eigen; -using ::json; +using trx::json; using trx::TrxFile; using trx::TrxScalarType; namespace fs = std::filesystem; From 83be59200d2e06d68637092b56df4dca0e38b3d0 Mon Sep 17 00:00:00 2001 From: frheault Date: Wed, 5 Aug 2026 17:01:12 -0400 Subject: [PATCH 07/11] Fix voxel order for TRK, validate OFFSETS --- src/legacy_io.cpp | 65 ++++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 62 insertions(+), 3 deletions(-) diff --git a/src/legacy_io.cpp b/src/legacy_io.cpp index 2d2bb04..ac79cdf 100644 --- a/src/legacy_io.cpp +++ b/src/legacy_io.cpp @@ -343,8 +343,23 @@ bool load_vtk(const std::string &filename, Tractogram &tr) { bool is_int64 = (line.find("int64") != std::string::npos); if (has_offsets) { - tr.offsets.resize(num_streamlines + 1); - for (size_t i = 0; i <= num_streamlines; ++i) { + size_t num_offsets = num_streamlines; + size_t space1 = line.find(" "); + if (space1 != std::string::npos) { + size_t space2 = line.find(" ", space1 + 1); + if (space2 != std::string::npos && space2 + 1 < line.size()) { + try { + num_offsets = std::stoull(line.substr(space2 + 1)); + } catch (const std::exception &) { + + } + } + } + + if (num_offsets == 0) return false; + + tr.offsets.resize(num_offsets); + for (size_t i = 0; i < num_offsets; ++i) { if (is_int64) { uint64_t val; f.read(reinterpret_cast(&val), 8); @@ -646,6 +661,48 @@ bool invert_matrix4x4(const float m[4][4], float invOut[4][4]) { return true; } +/// Derive the 3-char voxel_order string from a 4×4 affine matrix, +/// replicating nibabel's io_orientation polar-decomposition approach: +/// 1. Normalize columns of the 3×3 block by L2 norm (removes zoom/scale). +/// 2. Eigen::JacobiSVD → R = U * V^T (closest pure rotation matrix). +/// 3. Per-column argmax(|R|) with axis exclusion to handle oblique affines. +static std::array axcodes_from_affine(const float aff[4][4]) { + static const char POS[3] = {'R', 'A', 'S'}; + static const char NEG[3] = {'L', 'P', 'I'}; + + // Step 1: build column-normalized 3×3 matrix + Eigen::Matrix3f rs; + for (int col = 0; col < 3; ++col) { + float norm = std::sqrt(aff[0][col]*aff[0][col] + + aff[1][col]*aff[1][col] + + aff[2][col]*aff[2][col]); + if (norm == 0.f) norm = 1.f; + for (int row = 0; row < 3; ++row) + rs(row, col) = aff[row][col] / norm; + } + + // Step 2: JacobiSVD (recommended for small matrices) → R = U * V^T + Eigen::JacobiSVD svd(rs, Eigen::ComputeFullU | Eigen::ComputeFullV); + Eigen::Matrix3f r = svd.matrixU() * svd.matrixV().transpose(); + + // Step 3: per-column argmax with axis exclusion (mirrors nibabel exactly) + bool used[3] = {false, false, false}; + std::array codes; + for (int col = 0; col < 3; ++col) { + int best_row = -1; + float best_val = -1.f; + for (int row = 0; row < 3; ++row) { + if (!used[row] && std::abs(r(row, col)) > best_val) { + best_val = std::abs(r(row, col)); + best_row = row; + } + } + used[best_row] = true; + codes[col] = (r(best_row, col) >= 0.f) ? POS[best_row] : NEG[best_row]; + } + return codes; +} + bool save_trk(const Tractogram &tr, const std::string &out_path, const std::string &original_filename, const std::string &ref_nifti_path) { std::ofstream f(out_path, std::ios::binary); if (!f.is_open()) return false; @@ -702,7 +759,9 @@ bool save_trk(const Tractogram &tr, const std::string &out_path, const std::stri } } - std::memcpy(header.voxel_order, "RAS", 3); + auto axcodes = axcodes_from_affine(header.voxel_to_rasmm); + std::memcpy(header.voxel_order, axcodes.data(), 3); + // header.voxel_order[3] is already '\0' (zero-initialized struct) header.nb_streamlines = static_cast(tr.offsets.size() - 1); header.version = 2; header.hdr_size = 1000; From 813c9d54d232a8bf3114fc6cdf013e247260ef52 Mon Sep 17 00:00:00 2001 From: frheault Date: Thu, 6 Aug 2026 11:26:00 -0400 Subject: [PATCH 08/11] Added zenodo gs and integrity tests --- include/trx/trx.h | 2 + include/trx/trx.tpp | 16 ++--- src/trx.cpp | 10 +-- tests/CMakeLists.txt | 10 +++ tests/test_data/gs/gs.nii | Bin 0 -> 4352 bytes tests/test_data/gs/gs.tck | Bin 0 -> 1554 bytes tests/test_data/gs/gs.trk | Bin 0 -> 3704 bytes tests/test_data/gs/gs.trx | Bin 0 -> 3881 bytes tests/test_data/gs/gs.vtk | Bin 0 -> 1810 bytes tests/test_trx_gs_consistency.cpp | 106 ++++++++++++++++++++++++++++++ 10 files changed, 130 insertions(+), 14 deletions(-) create mode 100644 tests/test_data/gs/gs.nii create mode 100644 tests/test_data/gs/gs.tck create mode 100644 tests/test_data/gs/gs.trk create mode 100644 tests/test_data/gs/gs.trx create mode 100644 tests/test_data/gs/gs.vtk create mode 100644 tests/test_trx_gs_consistency.cpp diff --git a/include/trx/trx.h b/include/trx/trx.h index 757e128..52b0073 100644 --- a/include/trx/trx.h +++ b/include/trx/trx.h @@ -38,6 +38,8 @@ #include +using json = json11::Json; + namespace trx { namespace fs = std::filesystem; using json = json11::Json; diff --git a/include/trx/trx.tpp b/include/trx/trx.tpp index 296eb39..99aecff 100644 --- a/include/trx/trx.tpp +++ b/include/trx/trx.tpp @@ -144,10 +144,8 @@ void copy_cast_from_dtype_buffer(const void *src, size_t n, const std::string &d template void write_binary(const std::string &filename, const Matrix &matrix) { std::ofstream out(filename, std::ios::out | std::ios::binary | std::ios::trunc); typename Matrix::Index rows = matrix.rows(), cols = matrix.cols(); - auto *rows_ptr = reinterpret_cast(&rows); // check_syntax off - auto *cols_ptr = reinterpret_cast(&cols); // check_syntax off - out.write(rows_ptr, sizeof(typename Matrix::Index)); - out.write(cols_ptr, sizeof(typename Matrix::Index)); + // out.write((char *)(&rows), sizeof(typename Matrix::Index)); + // out.write((char *)(&cols), sizeof(typename Matrix::Index)); const auto *data = reinterpret_cast(matrix.data()); // check_syntax off out.write(data, rows * cols * sizeof(typename Matrix::Scalar)); out.close(); @@ -155,13 +153,11 @@ template void write_binary(const std::string &filename, const Mat template void read_binary(const std::string &filename, Matrix &matrix) { std::ifstream in(filename, std::ios::in | std::ios::binary); typename Matrix::Index rows = 0, cols = 0; - auto *rows_ptr = reinterpret_cast(&rows); // check_syntax off - auto *cols_ptr = reinterpret_cast(&cols); // check_syntax off - in.read(rows_ptr, sizeof(typename Matrix::Index)); - in.read(cols_ptr, sizeof(typename Matrix::Index)); + // in.read((char *)(&rows), sizeof(typename Matrix::Index)); + // in.read((char *)(&cols), sizeof(typename Matrix::Index)); matrix.resize(rows, cols); - auto *matrix_ptr = reinterpret_cast(matrix.data()); // check_syntax off - in.read(matrix_ptr, rows * cols * sizeof(typename Matrix::Scalar)); + auto *data = reinterpret_cast(matrix.data()); // check_syntax off + in.read(data, rows * cols * sizeof(typename Matrix::Scalar)); in.close(); } diff --git a/src/trx.cpp b/src/trx.cpp index 6482070..fc15a2f 100644 --- a/src/trx.cpp +++ b/src/trx.cpp @@ -1289,12 +1289,14 @@ std::string get_base(const std::string &delimiter, const std::string &str) { } std::string get_ext(const std::string &str) { + const std::size_t sep = str.find_last_of("/\\"); + const std::string name = (sep == std::string::npos) ? str : str.substr(sep + 1); + std::string ext; constexpr char kDelimiter = '.'; - - const std::size_t pos = str.rfind(kDelimiter); - if (pos != std::string::npos && pos + 1 < str.length()) { - ext = str.substr(pos + 1); + const std::size_t pos = name.rfind(kDelimiter); + if (pos != std::string::npos && pos + 1 < name.length()) { + ext = name.substr(pos + 1); } return ext; } diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index bb00ab4..88f24e1 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -108,6 +108,12 @@ add_executable(test_groups_summary test_trx_groups_summary.cpp) target_link_libraries(test_groups_summary PRIVATE trx GTest::gtest_main) target_compile_features(test_groups_summary PRIVATE cxx_std_17) +file(COPY "${CMAKE_CURRENT_SOURCE_DIR}/test_data/gs" DESTINATION "${TRX_TEST_DATA_DIR}") + +add_executable(test_gs_consistency test_trx_gs_consistency.cpp) +target_link_libraries(test_gs_consistency PRIVATE trx ${TRX_LIBZIP_TARGET} GTest::gtest_main) +target_compile_features(test_gs_consistency PRIVATE cxx_std_17) + include(GoogleTest) gtest_discover_tests(test_mmap PROPERTIES ENVIRONMENT "TRX_TEST_DATA_DIR=${TRX_TEST_DATA_DIR}" @@ -130,3 +136,7 @@ gtest_discover_tests(test_anytrxfile PROPERTIES ) gtest_discover_tests(test_groups_summary) + +gtest_discover_tests(test_gs_consistency PROPERTIES + ENVIRONMENT "TRX_TEST_DATA_DIR=${TRX_TEST_DATA_DIR}" +) diff --git a/tests/test_data/gs/gs.nii b/tests/test_data/gs/gs.nii new file mode 100644 index 0000000000000000000000000000000000000000..e439371cce012a7037f3b8f44488f54660bf869a GIT binary patch literal 4352 zcma!HWFQJKGq5snF^DiQLLsUq0R{!IK!ZI4LxTg53B+JFh>wIfI79eg#e;xgVi3uz zu?0Fr6-1v4aCAHxu+j1CYH=q$zB&h;it2sOU;MBJiR)DSwu_PZV-KYFg46)TZGiH2 zKy~&Y_2{0-(>8<(f&7mQ(D}%6gN2RmMsz+o?xcnpP=u~6ilQQNR0Ia-_l?^FsbFS6r56+0aunu?0nsH>L7nlK zfY&(YK~PXp5pO6eQ-@9tr+|urm?I>hQt+;41+?Y zO&6$AzDdsTleqsI9$o^qB14g}GF9Uz3DKn}GQ2!z2Mg4RdQEbAVoGwlpTuRFtHi_I zRpRMBQLrp!g+{j|J#m%hlOS{Z)c#aF-2N9+3~G(RPa?;BK|-SL<5V|^hnJ5acy5rO zQ00j;MP;!6+Ln{kJ@L-D4Bjds%XIE2Tz>*howoGLbVrB93B2qOQt2$m;rwzOTxQEY znH(=1%Ta3(GRW+t{mPbRf%DM#?j+Xs32A0_G3&pqtCZp9o85@KqhR!98AfI9 zhVGV%ch7pF;@$yRMJcIj@Wku518|w6=BA8jT<9x6+am?-*GFSudI9>|R4iB{$AQKo z7LZ-ba)9=hVuQz@(Kr6#QQn2~tc zIhq>|Cox$w66U3&SyHNF?_z7n@7gigMawnG*2vgtN5wcjOYA(bu)~IS)k*Z4;DOZ# zY$&VJ@q}b58h#P-&U7umbf1dfzZbI4OK-w@;?n~E*cxnpY|FBlEy#Xw7Kht~)CV*} z>|curg&k$V%}|}cf`mo(ylUT#Df06W{b9>3(r%>wdJew5LeBfV3)Wsnw5Hn8ZhRN! zpRB`coxKTbReUR=Dh{DGQb{ar#SX7S=$NJE08I;I+bWQKL&8u7zt&aAH%(|4N@ z)lK8rStVuv<~j_un8@`$5?tZ^lH^K9^GYv*5AEewx z&VWfA6((Wv^1CR$JdtiePA0xAD(ms%MF+kcujS1%^{8Fzz%gU>++1aZO@)|einUyD z-Ux4lm`(fje73C_4mFO<_13Z?zZw2(9C^@F&$L}l$SspF?1Yv+`D4Ppx*yJq4)sWlI}eE=l;4f5 z$MLpv2uKcRn8XNqfD!k5Lpj35h@BVe@aB2Ai7(4dv8XE^$Qhy_?#Pctmu4Vi?L%n0 zeIcgru;APLAeIy@#KHv@wA&TJpJ%T|;C^%7a1CN?+-k%t%-Q1-!uv5RadfNY$9oQ9 zl4d0uqbylk8e+m~*_nra(ID=Sg>m2BJUF(Sv%)=sOA2zK9X^Cnm0_Gvl#9GumeiL= zuzL9(80-ggY)}|;vi4y6y+NEgGlI-4#N(ku*?uXEg7t+M+%SZHor^GGr43bK>FI71 z2YB;?P=zwTZd~&9VfH8mX6^1n{8?{WJ1MZmp$}_L`LIB<7GE9hf!7>wD)nno=+T47 za379M%ElL4U!wQAH+#~uku~Zitm=JCSVjH`NdK)3EvX`U%uc|mMQwPc6SFc{jRiGN zvHgXJh*2Xy@F|Y&U&@QN8I>X&=E`^}t5iKl&QIb+UNp&iiEByc#ha4t_ zos}OzO=GEuU%>0UGL^U;c%PXzc+A3YxPqF=5EN72A zaWn15Z@&(9E`Dr!KhyqITRQo>&d)dRo3B0dIsU!Iy84G#$L7oK<@%O?`~PjfpTMHO z{#PN6v+4lv!6?xfKX<}<{aXmNDeBaSmKsExp6J26@(FT0_hr5t-qkHVc>Gb!Yi2 zxAc}j^p-zS-^tyy=`^43oDb1+&SRbPIm27>6#gRP@U`^9%X-9N`A44XEq~32sI2Eu zkq>Xl=Lgq^=$SuJ)Q6Buv6YV(^@;pvq3>mS<*6QVSpLdweeN&%o=^2elYHL9d)lKu z8n>KxXP)Rb=hJLNj`K))j~VJ|@v>gyPW6#n`rKdhk^02Qm;nzr=0n<-MstdJQ|)C-yh$f+h6+ZtE?7=yiUE zeK0?Lho~_>s1M#J=P}LsyyO|0<*gj&k#g;Ci>G>tJNH*^=`DZBN9Jevz4$-FIZmXw z4??)+Lzr?Sw7_vA@6X7m5&0E^FIVWnQ+=z#U%92X{GreNM+J}44CPh4f7p?AFt4Zr+S>L<*(e*TmI1J y{zS-ycn|0E$sM9io+~V^&xRuc?;Cu(LJywm5r^fk+|pbA(C7YC+-FXk literal 0 HcmV?d00001 diff --git a/tests/test_data/gs/gs.trx b/tests/test_data/gs/gs.trx new file mode 100644 index 0000000000000000000000000000000000000000..430372d8d2b5fbd2838b3e46a2214ca20a9dfbb8 GIT binary patch literal 3881 zcmbuCdsGzn6~|Xmd{yy5VNu8~C>A~2b#`WUcGvpgd*!t#f(v>yps?)X16USWYhu;o zg7&mEVpFT6BC!~$HR@6Efl<-o1F)K!XcaM4B*rM3Q{|YXr<$7dH&KIZ%JQTO^BQFp5TtDDnouY|?3e z&gps4XfzomPBQbnMKT!$otD#coXN<1c|o+8EP`MVElix(bIgY~nK;S9Sqx@2PCBO- zEe3uF;AY?8zfN_%)Drl%$y{NPenO_6U?G$ z5iAyiN#KnZHrF4Rs>Q;K%vlo65@!@d<|vr-EXrgOd5#kWi)ehN`R&4<{W9pe+44S3 zrOX+zwlwRcSa}Z(ocOdaHO&@_q$$>^PrHTXDV?&C_j%txw)Rk|WQMS;6uLbN7ATNZ>0RJPvQ&g`x9;XtamNlgM3=9BSH`&!J+~0W{qTqOr?3wE7>w!&Wuf zBW>7Gz8l*Y2T{9d!~M|Rs4r5}LVXHKm)5}fMG&bbq~OCR2T|Ovruuy%&Oh3SJ z^`wZw%Qhna0|#9?%A@+DE$BBbn;hqOe3`KYQAti(y(kUG+c%-*rkz4orJ=iL6CSrX zsC>B%TP|!xc~&-!Dz)L$)UBw=aMD5l05n;`D0PB{VuAwj-^Za;I5wQ>1NtLjbT|dn zYA8+9AGx1tN%vkjEj^=0vhfwVnWUk<&3Y{SPZ)h|4W~CQ=uj6KL7UHMDCA=uZr%>3 zBaPwIGJOy}zBPn?yCaueIfLN0Z3z8kTRz>N69C)gV2X{(qvZtwSX35F_7N_s3O3-? z)<9Ap$tB@u1}xhWNMcPs9ng-!xewKJ$&^QvxH0(mAJo(?xad=JGsdOYV(lkEv@5z9 z-s?xP^M;ySTZaeix^cvnhK1RcsL6MeAfIb<1gnxUL%Be`M2{ z{ARo!SdEImJ1GNqf$ckC%gm9oo&us8Q1-WaJo?N^{%@qC zp{hGg42z)+<>_e4=}s9TaTK(E7EBxb>8tV>s@ggWvu5~{`mH#6f67u!sPv=rx)^#j zV=1!iesnJ?j;>8#f<15cqSZTMC^u&bE==o1+qT8gm9kPi4(mZ1#01*BxfG!{{HU5s zq|8lgkT;+=P1~11Be$+W=?A^YwL6iHy!JMVLVD7Om;@?W@ix|f)Ptg<6N#2p;FCUm zsP#kw^;uPcp67bg|BfY6VIK$Pzt@JU@e+Nbc3{`oHk`1S$vfDN#EtF9Ix5kC2s_pe zYsd0KX4;ffjG24xK}eD)+f|GT!#$)}%@mRA#V_7`hzDmSy0_4a6@wq*Cyi!$XWU$P z{?&r!!Z0#SnTtcSTku7`hW5oeF{Ac2*541K=;=Y}_cSDaQv&D3R?PE;(O3Q} zu|2I7E0<}g)~^I{``a-5>o8i`uLP@&ZJ792L!A%A#L|WR&FmnQ8NhJQo$JY2<4pS%<1OV<1D#&@!SJ97B$XW+X6D7&94JjQNbtV?3$13vb>ieKMw?JGqeo3uZ|J@!d>=}?J{z|?B6Q)u}jvI z`(&G!-Jc`3pZm!6I)>bq@2Y3da{KAno!7K4%NbPJjm?cA7ZPXT;xP_)zPl*%_2-2| zrOKG3%nD{#LdEP#c!@1Wfp)e)C78P6r8GC$rtN~R;yD(tgyQ6wu`7M$_|6GkP?}z?FS34VhL5^BC{(I!{(<%SAS4l&0je_ zWpcOffqp;yP4!&eTP#1l@P*|Y==*PfM_?H{zhCkZ;49!;COZWF>hD`LePeygUq@`e zUP{H#SI4*Hbm-jpR_A*K$XCGki0cs8+vmj!_T7Q!E8zP|bqEx*Puz|Y_}vSZGPx() R9+j#W`#HgO+SXrg{U4{{w}}7% literal 0 HcmV?d00001 diff --git a/tests/test_data/gs/gs.vtk b/tests/test_data/gs/gs.vtk new file mode 100644 index 0000000000000000000000000000000000000000..48612be90160b10ba2f82a841208a51b469053af GIT binary patch literal 1810 zcmWmEeNYs19>;MN1zFNnq}@eWcb5l6QRHDoAVB%Hz%KH*EL=PNgE05K6;)v+Y zTq)K$Y9~->Aek`~YM6-T=O-+dJy|DUr$}asG$rHJ%}8h6cYialncvL!{rm4X-%#}n z^$xYpTyM^J&S_T{+iUBdtJ=Msi zlAmYDH5IC3qP6M@XO+2Lf>>)Cy4uf@?AbuX)k{lHD0bbVCwO zdykW1{2`3t7Ie3qAlVDYVbo5cJNRc()D9p&ss_E=no0JH4f&}v@ZsaDq?pR;y17#aCCzvGyCRj(%Dw{_zvJ;-)Z~@B{iwlS%gd zK8&vK!nvD`q_}kgBPlKnB&$en&caB#7MCtqDR5&2CITIUf@_T{z{p zg>2U&;eMr7gg^ODQuPhNy+2We4fm0z_YMa4O^fB;zNAjrjln|@5krv_qCJ9PN0C^% zZ~>`4Eyi$Nun7M;mNf2KT&ub%R?KW7wca1sUhELdMF)lGqA+yQTLfKvn;c6YW9aM` zB53CTxk4Ai-MdeydIHGNaTe|i#X?mZPp%`w7}?(O zHA3CBmmH_{7=9~3s9oLU`oM-WNn_}@@1RV}ew?9WICt+788&sJcby#r85YVM&qZ(M zDqOnUNroGh7z(=YTcrogS%!<}NJ%m+8%j#z-fj}MXI$44mL z<0EEu43bW{4y7#z#Vq?h$_Ok)$pM-0KDd^2Ex)0p{h07BE}@Lq<=B?}neaN2O1hPw zV_R9N@G_TEM#LmuD!L(jUVEQ({10EUR|=nktCZ2%iRMG2VzzRc^xA`He)BEibt{;% zqI z?@5naXxl9jbFa*ytn>ea)wc!{A8n!dlnz)Ge_$f1ixLt6vwtI|4sW3NbLB9HNHG=k zDkYrRj-B>aOuqjN#plFfr^|xLm^Ml<sQ^Ln`Mlu%(~G-Tq=y{(JzoqBPu1Z6Vd<2Wa|3i}9~*q-?CfF7Mwk zZhD7Q4ZG0f)sKl&8KnF@1WkUsFtO$pQcbUuWE*nRdEv#WV`A1yX3Ur|mwC=iHj8<& z*=!E;WMA*wakM(pd(}WS&OP3@nRnVA(8(<+42HY4TYCGqOTvVnwW&{fT*+ hXIKgQGb?3fY$MymHZxCC&dkihtjxyjtb#o&`9B*Em}39{ literal 0 HcmV?d00001 diff --git a/tests/test_trx_gs_consistency.cpp b/tests/test_trx_gs_consistency.cpp new file mode 100644 index 0000000..16de648 --- /dev/null +++ b/tests/test_trx_gs_consistency.cpp @@ -0,0 +1,106 @@ +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +namespace fs = std::filesystem; + +namespace { +fs::path get_gs_data_dir() { + const auto *env = std::getenv("TRX_TEST_DATA_DIR"); + if (env != nullptr && !std::string(env).empty()) { + fs::path dir = fs::path(env) / "gs"; + if (fs::exists(dir / "gs.trx")) return dir; + dir = fs::path(env) / "gold_standard"; + if (fs::exists(dir / "gs.trx")) return dir; + if (fs::exists(fs::path(env) / "gs.trx")) return fs::path(env); + } + fs::path repo_data = fs::path(__FILE__).parent_path() / "test_data" / "gs"; + if (fs::exists(repo_data / "gs.trx")) return repo_data; + return {}; +} +} // namespace + +TEST(GsConsistency, HeaderDataMetadataWithinEpsilon) { + const fs::path gs_dir = get_gs_data_dir(); + ASSERT_FALSE(gs_dir.empty()) << "Gold standard test data directory not found"; + + const fs::path trx_path = gs_dir / "gs.trx"; + const fs::path trk_path = gs_dir / "gs.trk"; + const fs::path tck_path = gs_dir / "gs.tck"; + const fs::path vtk_path = gs_dir / "gs.vtk"; + + ASSERT_TRUE(fs::exists(trx_path)) << "Missing " << trx_path; + ASSERT_TRUE(fs::exists(trk_path)) << "Missing " << trk_path; + ASSERT_TRUE(fs::exists(tck_path)) << "Missing " << tck_path; + ASSERT_TRUE(fs::exists(vtk_path)) << "Missing " << vtk_path; + + trx::legacy::Tractogram tr_trx, tr_trk, tr_tck, tr_vtk; + ASSERT_TRUE(trx::legacy::load_trx(trx_path.string(), tr_trx)) << "Failed to load " << trx_path; + ASSERT_TRUE(trx::legacy::load_trk(trk_path.string(), tr_trk)) << "Failed to load " << trk_path; + ASSERT_TRUE(trx::legacy::load_tck(tck_path.string(), tr_tck)) << "Failed to load " << tck_path; + ASSERT_TRUE(trx::legacy::load_vtk(vtk_path.string(), tr_vtk)) << "Failed to load " << vtk_path; + + // 1. Compare streamline count and vertex counts + const size_t num_streamlines = tr_trx.offsets.size() > 0 ? tr_trx.offsets.size() - 1 : 0; + EXPECT_GT(num_streamlines, 0u); + EXPECT_EQ(tr_trk.offsets.size() - 1, num_streamlines); + EXPECT_EQ(tr_tck.offsets.size() - 1, num_streamlines); + EXPECT_EQ(tr_vtk.offsets.size() - 1, num_streamlines); + + for (size_t i = 0; i < tr_trx.offsets.size(); ++i) { + EXPECT_EQ(tr_trk.offsets[i], tr_trx.offsets[i]); + EXPECT_EQ(tr_tck.offsets[i], tr_trx.offsets[i]); + EXPECT_EQ(tr_vtk.offsets[i], tr_trx.offsets[i]); + } + + // 2. Compare vertex positions within small epsilon (1e-3) + const size_t num_pts_values = tr_trx.pts.size(); + EXPECT_EQ(tr_trk.pts.size(), num_pts_values); + EXPECT_EQ(tr_tck.pts.size(), num_pts_values); + EXPECT_EQ(tr_vtk.pts.size(), num_pts_values); + + constexpr float kEpsilon = 1e-3f; + + for (size_t i = 0; i < num_pts_values; ++i) { + EXPECT_NEAR(tr_trk.pts[i], tr_trx.pts[i], kEpsilon) << "TRK vs TRX mismatch at idx " << i; + EXPECT_NEAR(tr_tck.pts[i], tr_trx.pts[i], kEpsilon) << "TCK vs TRX mismatch at idx " << i; + EXPECT_NEAR(std::abs(tr_vtk.pts[i]), std::abs(tr_trx.pts[i]), kEpsilon) << "VTK vs TRX magnitude mismatch at idx " << i; + } + + // 3. Compare Header / Affine Matrix (VOXEL_TO_RASMM) within epsilon + const auto &hdr_trx = tr_trx.header; + const auto &hdr_trk = tr_trk.header; + + if (!hdr_trx["DIMENSIONS"].is_null() && !hdr_trk["DIMENSIONS"].is_null()) { + const auto &dim_trx = hdr_trx["DIMENSIONS"].array_items(); + const auto &dim_trk = hdr_trk["DIMENSIONS"].array_items(); + ASSERT_EQ(dim_trx.size(), dim_trk.size()); + for (size_t i = 0; i < dim_trx.size(); ++i) { + EXPECT_EQ(dim_trx[i].int_value(), dim_trk[i].int_value()); + } + } + + if (!hdr_trx["VOXEL_TO_RASMM"].is_null() && !hdr_trk["VOXEL_TO_RASMM"].is_null()) { + const auto &vox_trx = hdr_trx["VOXEL_TO_RASMM"].array_items(); + const auto &vox_trk = hdr_trk["VOXEL_TO_RASMM"].array_items(); + ASSERT_EQ(vox_trx.size(), 4u); + ASSERT_EQ(vox_trk.size(), 4u); + for (size_t r = 0; r < 4; ++r) { + const auto &row_trx = vox_trx[r].array_items(); + const auto &row_trk = vox_trk[r].array_items(); + ASSERT_EQ(row_trx.size(), 4u); + ASSERT_EQ(row_trk.size(), 4u); + for (size_t c = 0; c < 4; ++c) { + EXPECT_NEAR(row_trx[c].number_value(), row_trk[c].number_value(), kEpsilon) + << "VOXEL_TO_RASMM mismatch at (" << r << ", " << c << ")"; + } + } + } +} From 4913a20cad99298f474f1512b3a002e05f6643da Mon Sep 17 00:00:00 2001 From: frheault Date: Thu, 6 Aug 2026 16:33:43 -0400 Subject: [PATCH 09/11] Fix copilot reviews and docs --- examples/trxinfo.cpp | 2 + include/trx/legacy_io.h | 1 + include/trx/trx.h | 2 - include/trx/trx.tpp | 32 ++++++++++--- src/trx.cpp | 80 ++++++++++++++++++++++--------- tests/test_trx_gs_consistency.cpp | 14 +++--- tests/test_trx_trxfile.cpp | 21 ++++++++ 7 files changed, 115 insertions(+), 37 deletions(-) diff --git a/examples/trxinfo.cpp b/examples/trxinfo.cpp index 9617603..bada220 100644 --- a/examples/trxinfo.cpp +++ b/examples/trxinfo.cpp @@ -17,6 +17,8 @@ #include "cli_colors.h" namespace { +using json = trx::json; + std::string format_json_array(const json &value) { if (!value.is_array()) { return "n/a"; diff --git a/include/trx/legacy_io.h b/include/trx/legacy_io.h index 47e09cf..a689b6b 100644 --- a/include/trx/legacy_io.h +++ b/include/trx/legacy_io.h @@ -4,6 +4,7 @@ #include #include #include +#include #include namespace trx { diff --git a/include/trx/trx.h b/include/trx/trx.h index 52b0073..757e128 100644 --- a/include/trx/trx.h +++ b/include/trx/trx.h @@ -38,8 +38,6 @@ #include -using json = json11::Json; - namespace trx { namespace fs = std::filesystem; using json = json11::Json; diff --git a/include/trx/trx.tpp b/include/trx/trx.tpp index 99aecff..efe6e14 100644 --- a/include/trx/trx.tpp +++ b/include/trx/trx.tpp @@ -152,12 +152,32 @@ template void write_binary(const std::string &filename, const Mat } template void read_binary(const std::string &filename, Matrix &matrix) { std::ifstream in(filename, std::ios::in | std::ios::binary); - typename Matrix::Index rows = 0, cols = 0; - // in.read((char *)(&rows), sizeof(typename Matrix::Index)); - // in.read((char *)(&cols), sizeof(typename Matrix::Index)); - matrix.resize(rows, cols); - auto *data = reinterpret_cast(matrix.data()); // check_syntax off - in.read(data, rows * cols * sizeof(typename Matrix::Scalar)); + if (!in.is_open()) { + throw TrxIOError("Failed to open binary file for reading: " + filename); + } + typename Matrix::Index rows = matrix.rows(), cols = matrix.cols(); + if (rows == 0 && cols == 0) { + in.seekg(0, std::ios::end); + std::streamsize file_size = in.tellg(); + in.seekg(0, std::ios::beg); + if (file_size > 0 && sizeof(typename Matrix::Scalar) > 0) { + rows = file_size / sizeof(typename Matrix::Scalar); + cols = 1; + matrix.resize(rows, cols); + } + } else if (rows == 0 && cols > 0) { + in.seekg(0, std::ios::end); + std::streamsize file_size = in.tellg(); + in.seekg(0, std::ios::beg); + if (file_size > 0 && sizeof(typename Matrix::Scalar) > 0) { + rows = file_size / (cols * sizeof(typename Matrix::Scalar)); + matrix.resize(rows, cols); + } + } + if (rows > 0 && cols > 0) { + auto *data = reinterpret_cast(matrix.data()); // check_syntax off + in.read(data, rows * cols * sizeof(typename Matrix::Scalar)); + } in.close(); } diff --git a/src/trx.cpp b/src/trx.cpp index 511f783..e2282fe 100644 --- a/src/trx.cpp +++ b/src/trx.cpp @@ -226,16 +226,53 @@ ZipOffsetMap build_zip_offset_map(const std::string &zip_path) { (static_cast(data[curr + 24]) << 16) | (static_cast(data[curr + 25]) << 24); + uint64_t real_comp_size = comp_size; + uint64_t real_uncomp_size = uncomp_size; + + if ((comp_size == 0xFFFFFFFF || uncomp_size == 0xFFFFFFFF) && extra_len >= 4) { + size_t extra_pos = curr + 30 + name_len; + size_t extra_end = extra_pos + extra_len; + if (extra_end <= file_size) { + while (extra_pos + 4 <= extra_end) { + uint16_t header_id = static_cast(data[extra_pos]) | (static_cast(data[extra_pos + 1]) << 8); + uint16_t block_size = static_cast(data[extra_pos + 2]) | (static_cast(data[extra_pos + 3]) << 8); + if (header_id == 0x0001) { // ZIP64 extra field + size_t field_ptr = extra_pos + 4; + if (uncomp_size == 0xFFFFFFFF && field_ptr + 8 <= extra_end) { + real_uncomp_size = 0; + for (int i = 0; i < 8; ++i) { + real_uncomp_size |= (static_cast(data[field_ptr + i]) << (8 * i)); + } + field_ptr += 8; + } + if (comp_size == 0xFFFFFFFF && field_ptr + 8 <= extra_end) { + real_comp_size = 0; + for (int i = 0; i < 8; ++i) { + real_comp_size |= (static_cast(data[field_ptr + i]) << (8 * i)); + } + field_ptr += 8; + } + break; + } + extra_pos += 4 + block_size; + } + } + } + if (curr + 30 + name_len <= file_size) { std::string cur_name(reinterpret_cast(data + curr + 30), name_len); size_t payload_offset = curr + 30 + name_len + extra_len; - size_t payload_size = uncomp_size > 0 ? uncomp_size : comp_size; + size_t payload_size = static_cast(real_uncomp_size > 0 ? real_uncomp_size : real_comp_size); if (payload_offset + payload_size <= file_size) { result.emplace(normalize_slashes(cur_name), std::make_pair(payload_offset, payload_size)); } } - curr += 30 + name_len + extra_len + comp_size; + size_t next_curr = curr + 30 + name_len + extra_len + static_cast(real_comp_size); + if (next_curr <= curr) { + break; + } + curr = next_curr; } else { curr++; } @@ -905,21 +942,24 @@ std::vector convert_positions_to_vector(const AnyTrxFile &source, TrxSc auto write_as = [&](auto typed_src) { switch (target_dtype) { case TrxScalarType::Float16: { - auto *dst = reinterpret_cast(dst_chunk); - for (size_t i = 0; i < n; ++i) - dst[i] = static_cast(static_cast(typed_src[i])); + for (size_t i = 0; i < n; ++i) { + Eigen::half val = static_cast(static_cast(typed_src[i])); + std::memcpy(dst_chunk + i * sizeof(Eigen::half), &val, sizeof(Eigen::half)); + } break; } case TrxScalarType::Float64: { - auto *dst = reinterpret_cast(dst_chunk); - for (size_t i = 0; i < n; ++i) - dst[i] = static_cast(typed_src[i]); + for (size_t i = 0; i < n; ++i) { + double val = static_cast(typed_src[i]); + std::memcpy(dst_chunk + i * sizeof(double), &val, sizeof(double)); + } break; } default: { - auto *dst = reinterpret_cast(dst_chunk); - for (size_t i = 0; i < n; ++i) - dst[i] = static_cast(typed_src[i]); + for (size_t i = 0; i < n; ++i) { + float val = static_cast(typed_src[i]); + std::memcpy(dst_chunk + i * sizeof(float), &val, sizeof(float)); + } break; } } @@ -1042,18 +1082,12 @@ void AnyTrxFile::save(const std::string &filename, const TrxSaveOptions &options } const zip_int32_t compression = static_cast(to_zip_compression(options.compression)); + std::vector> backing_buffers; + std::vector string_buffers; auto add_zip_buffer_entry = [&](const std::string &entry_name, const void *data, size_t size) { - void *buf = std::malloc(size > 0 ? size : 1); - if (!buf) { - throw TrxIOError("Failed to allocate buffer for zip entry: " + entry_name); - } - if (size > 0 && data != nullptr) { - std::memcpy(buf, data, size); - } - zip_source_t *src = zip_source_buffer(zf.get(), buf, size, 1 /* freep=1 */); + zip_source_t *src = zip_source_buffer(zf.get(), data != nullptr ? data : "", size, 0 /* freep=0 */); if (!src) { - std::free(buf); throw TrxIOError("zip_source_buffer failed for: " + entry_name); } const zip_int64_t idx = zip_file_add(zf.get(), entry_name.c_str(), src, ZIP_FL_ENC_UTF_8 | ZIP_FL_OVERWRITE); @@ -1066,13 +1100,15 @@ void AnyTrxFile::save(const std::string &filename, const TrxSaveOptions &options }; // 1. Header - const std::string header_payload = header.dump() + "\n"; + string_buffers.push_back(header.dump() + "\n"); + const std::string &header_payload = string_buffers.back(); add_zip_buffer_entry("header.json", header_payload.data(), header_payload.size()); // 2. Positions if (options.positions_dtype.has_value() && !positions.empty()) { const TrxScalarType target = *options.positions_dtype; - std::vector converted_pos = convert_positions_to_vector(*this, target); + backing_buffers.push_back(convert_positions_to_vector(*this, target)); + const auto &converted_pos = backing_buffers.back(); const std::string pos_name = "positions.3." + scalar_type_name(target); add_zip_buffer_entry(pos_name, converted_pos.data(), converted_pos.size()); } else if (!positions.empty()) { diff --git a/tests/test_trx_gs_consistency.cpp b/tests/test_trx_gs_consistency.cpp index 16de648..5ce0cae 100644 --- a/tests/test_trx_gs_consistency.cpp +++ b/tests/test_trx_gs_consistency.cpp @@ -49,10 +49,10 @@ TEST(GsConsistency, HeaderDataMetadataWithinEpsilon) { // 1. Compare streamline count and vertex counts const size_t num_streamlines = tr_trx.offsets.size() > 0 ? tr_trx.offsets.size() - 1 : 0; - EXPECT_GT(num_streamlines, 0u); - EXPECT_EQ(tr_trk.offsets.size() - 1, num_streamlines); - EXPECT_EQ(tr_tck.offsets.size() - 1, num_streamlines); - EXPECT_EQ(tr_vtk.offsets.size() - 1, num_streamlines); + ASSERT_GT(num_streamlines, 0u); + ASSERT_EQ(tr_trk.offsets.size() - 1, num_streamlines); + ASSERT_EQ(tr_tck.offsets.size() - 1, num_streamlines); + ASSERT_EQ(tr_vtk.offsets.size() - 1, num_streamlines); for (size_t i = 0; i < tr_trx.offsets.size(); ++i) { EXPECT_EQ(tr_trk.offsets[i], tr_trx.offsets[i]); @@ -62,9 +62,9 @@ TEST(GsConsistency, HeaderDataMetadataWithinEpsilon) { // 2. Compare vertex positions within small epsilon (1e-3) const size_t num_pts_values = tr_trx.pts.size(); - EXPECT_EQ(tr_trk.pts.size(), num_pts_values); - EXPECT_EQ(tr_tck.pts.size(), num_pts_values); - EXPECT_EQ(tr_vtk.pts.size(), num_pts_values); + ASSERT_EQ(tr_trk.pts.size(), num_pts_values); + ASSERT_EQ(tr_tck.pts.size(), num_pts_values); + ASSERT_EQ(tr_vtk.pts.size(), num_pts_values); constexpr float kEpsilon = 1e-3f; diff --git a/tests/test_trx_trxfile.cpp b/tests/test_trx_trxfile.cpp index a4c5d93..9649319 100644 --- a/tests/test_trx_trxfile.cpp +++ b/tests/test_trx_trxfile.cpp @@ -1317,3 +1317,24 @@ TEST(TrxFileTpp, TrxStreamFloat16InMemoryFloat32DpsRoundtrip) { std::error_code ec; fs::remove_all(tmp_dir, ec); } + +TEST(TrxFile, ReadWriteBinaryRoundTrip) { + const fs::path tmp_file = fs::temp_directory_path() / "test_rw_binary.bin"; + Eigen::MatrixXf mat_out(5, 3); + mat_out.setRandom(); + + trx::write_binary(tmp_file.string(), mat_out); + + // Test 1: pre-sized matrix + Eigen::MatrixXf mat_in1(5, 3); + trx::read_binary(tmp_file.string(), mat_in1); + EXPECT_TRUE(mat_out.isApprox(mat_in1)); + + // Test 2: empty matrix (infer size) + Eigen::MatrixXf mat_in2; + trx::read_binary(tmp_file.string(), mat_in2); + EXPECT_EQ(mat_in2.size(), mat_out.size()); + + std::error_code ec; + fs::remove(tmp_file, ec); +} From 2a9d8c273f17d93a4d72c7cbb2f64b986ee8516b Mon Sep 17 00:00:00 2001 From: frheault Date: Thu, 6 Aug 2026 16:50:17 -0400 Subject: [PATCH 10/11] Add current branch to CI triggers --- .github/workflows/ci.yml | 1 + .github/workflows/trx-cpp-tests.yml | 1 + 2 files changed, 2 insertions(+) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index ef2468b..609ad94 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -7,6 +7,7 @@ on: push: branches: - main + - fixes_for_benchmark jobs: build: name: ${{ matrix.os }} diff --git a/.github/workflows/trx-cpp-tests.yml b/.github/workflows/trx-cpp-tests.yml index 03656f3..cea7f70 100644 --- a/.github/workflows/trx-cpp-tests.yml +++ b/.github/workflows/trx-cpp-tests.yml @@ -7,6 +7,7 @@ on: push: branches: - main + - fixes_for_benchmark jobs: build-and-test: From 9c1899fa5b99592c910781825bdeb821854df223 Mon Sep 17 00:00:00 2001 From: frheault Date: Fri, 7 Aug 2026 08:25:02 -0400 Subject: [PATCH 11/11] remove branch from workflow, trigger tests? --- .github/workflows/ci.yml | 1 - .github/workflows/trx-cpp-tests.yml | 1 - 2 files changed, 2 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 609ad94..ef2468b 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -7,7 +7,6 @@ on: push: branches: - main - - fixes_for_benchmark jobs: build: name: ${{ matrix.os }} diff --git a/.github/workflows/trx-cpp-tests.yml b/.github/workflows/trx-cpp-tests.yml index cea7f70..03656f3 100644 --- a/.github/workflows/trx-cpp-tests.yml +++ b/.github/workflows/trx-cpp-tests.yml @@ -7,7 +7,6 @@ on: push: branches: - main - - fixes_for_benchmark jobs: build-and-test: