diff --git a/CMakeLists.txt b/CMakeLists.txt index 123ee78..0c12e65 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 @@ -220,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/examples/trxinfo.cpp b/examples/trxinfo.cpp index abec187..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"; @@ -152,7 +154,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/legacy_io.h b/include/trx/legacy_io.h new file mode 100644 index 0000000..a689b6b --- /dev/null +++ b/include/trx/legacy_io.h @@ -0,0 +1,63 @@ +#ifndef TRX_LEGACY_IO_H +#define TRX_LEGACY_IO_H + +#include +#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 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); + +} // namespace legacy +} // namespace trx + +#endif // TRX_LEGACY_IO_H 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..efe6e14 100644 --- a/include/trx/trx.tpp +++ b/include/trx/trx.tpp @@ -152,14 +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; - 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)); - matrix.resize(rows, cols); - auto *matrix_ptr = reinterpret_cast(matrix.data()); // check_syntax off - in.read(matrix_ptr, 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(); } @@ -244,7 +262,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 +450,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 +464,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 +484,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 +527,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 +552,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 +582,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 +611,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/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 new file mode 100644 index 0000000..ac79cdf --- /dev/null +++ b/src/legacy_io.cpp @@ -0,0 +1,942 @@ +#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) { + 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 static_cast(bswap32(static_cast(i))); +} + +inline int16_t swap_int16(int16_t val) { + return static_cast(bswap16(static_cast(val))); +} + +inline int64_t swap_int64(int64_t val) { + return static_cast(bswap64(static_cast(val))); +} + +inline double swap_double(double d) { + uint64_t i; + std::memcpy(&i, &d, sizeof(i)); + i = bswap64(i); + std::memcpy(&d, &i, sizeof(d)); + return d; +} + + + +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(); + + 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); + 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 += point_stride; + } + offset += prop_bytes; + } + + 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 || 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); + 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); + 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; + } + break; + } + } + 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); + 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) { + try { + num_streamlines = std::stoull(line.substr(6, line.find(" ", 6) - 6)); + } catch (const std::exception &) { + return false; // malformed LINES count + } + 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) { + 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); + 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); + std::vector skip_buf; + + 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 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; +} + +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; + } + 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; + + 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 = header_to_use; + + 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; + } +} + +// 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; +} + +/// 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; + 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()) { + 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); + + // 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 (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 (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) { + 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]); + } + } + + 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; + + 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); + + 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); + + const char* p_n_pts = reinterpret_cast(&n_pts); + chunk.insert(chunk.end(), p_n_pts, p_n_pts + 4); + + for (size_t j = start; j < end; ++j) { + 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); + 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; + 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; + 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; + 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; + + // 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 diff --git a/src/trx.cpp b/src/trx.cpp index d53d207..e2282fe 100644 --- a/src/trx.cpp +++ b/src/trx.cpp @@ -188,6 +188,98 @@ 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); + + 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 = 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)); + } + } + 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++; + } + } + + return result; +} } // namespace std::string detect_positions_dtype(const std::string &path) { @@ -370,11 +462,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; } @@ -512,18 +808,83 @@ 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 == "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()); + + 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") { + 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 { // uint64 + normalize(uint64_t{}); + } + + 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); } } - 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."); } @@ -560,44 +921,45 @@ 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); - 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))); + 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: { - std::vector buf(n); - 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))); + 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: { - std::vector buf(n); - 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))); + 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; } } @@ -616,10 +978,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; @@ -636,95 +1037,140 @@ 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 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()))); - } - 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()))); - } 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()))); - } + std::vector> backing_buffers; + std::vector string_buffers; - std::unordered_set skip = {"header.json"}; - // Guard deletes the temp positions file after commit (or on exception). - TempFileGuard tmp_pos_guard; + auto add_zip_buffer_entry = [&](const std::string &entry_name, const void *data, size_t size) { + zip_source_t *src = zip_source_buffer(zf.get(), data != nullptr ? data : "", size, 0 /* freep=0 */); + if (!src) { + 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 + 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; - 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); + 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()) { + 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); + } + } + + // 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); + } + } + + // 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); + } + } + } + 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) { @@ -738,42 +1184,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 (!trx::fs::exists(final_header_path)) { - throw TrxFormatError("Missing header.json in output directory: " + final_header_path.string()); + + 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 (!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); + } + } } } } @@ -901,7 +1377,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'; } } @@ -1251,26 +1727,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..88f24e1 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.") @@ -102,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}" @@ -124,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 0000000..e439371 Binary files /dev/null and b/tests/test_data/gs/gs.nii differ diff --git a/tests/test_data/gs/gs.tck b/tests/test_data/gs/gs.tck new file mode 100644 index 0000000..1210080 Binary files /dev/null and b/tests/test_data/gs/gs.tck differ diff --git a/tests/test_data/gs/gs.trk b/tests/test_data/gs/gs.trk new file mode 100644 index 0000000..481b25d Binary files /dev/null and b/tests/test_data/gs/gs.trk differ diff --git a/tests/test_data/gs/gs.trx b/tests/test_data/gs/gs.trx new file mode 100644 index 0000000..430372d Binary files /dev/null and b/tests/test_data/gs/gs.trx differ diff --git a/tests/test_data/gs/gs.vtk b/tests/test_data/gs/gs.vtk new file mode 100644 index 0000000..48612be Binary files /dev/null and b/tests/test_data/gs/gs.vtk differ diff --git a/tests/test_trx_anytrxfile.cpp b/tests/test_trx_anytrxfile.cpp index 7dfe8e1..01de042 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; @@ -1129,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()); @@ -1139,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; diff --git a/tests/test_trx_gs_consistency.cpp b/tests/test_trx_gs_consistency.cpp new file mode 100644 index 0000000..5ce0cae --- /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; + 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]); + 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(); + 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; + + 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 << ")"; + } + } + } +} 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; 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); +}