diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index da058ae..3adc137 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -273,21 +273,37 @@ jobs: if: runner.os == 'Windows' uses: ilammy/msvc-dev-cmd@v1 - - name: Configure trx-cpp + - name: Configure trx-cpp (macOS) + if: runner.os == 'macOS' run: | + BREW_PREFIX="$(brew --prefix)" cmake -S . -B build \ -G Ninja \ -DTRX_USE_CONAN=OFF \ - -DTRX_BUILD_TESTS=ON \ - -DTRX_BUILD_EXAMPLES=ON \ + -DTRX_BUILD_TESTS=OFF \ + -DTRX_BUILD_EXAMPLES=OFF \ -DTRX_ENABLE_NIFTI=ON \ -DCMAKE_BUILD_TYPE=Release \ - ${{ runner.os == 'Windows' && format('-DCMAKE_TOOLCHAIN_FILE={0}/vcpkg/scripts/buildsystems/vcpkg.cmake', github.workspace) || '' }} + -DCMAKE_PREFIX_PATH="${BREW_PREFIX}" \ + -DCMAKE_INSTALL_PREFIX="${GITHUB_WORKSPACE}/install" + + - name: Configure trx-cpp (Windows) + if: runner.os == 'Windows' + run: > + cmake -S . -B build + -G Ninja + -DTRX_USE_CONAN=OFF + -DTRX_BUILD_TESTS=OFF + -DTRX_BUILD_EXAMPLES=OFF + -DTRX_ENABLE_NIFTI=ON + -DCMAKE_BUILD_TYPE=Release + -DCMAKE_TOOLCHAIN_FILE=${{ github.workspace }}/vcpkg/scripts/buildsystems/vcpkg.cmake + -DCMAKE_INSTALL_PREFIX=${{ github.workspace }}/install - name: Build and install trx-cpp run: cmake --build build --config Release --target install - - name: Configure downstream consumer + - name: Generate downstream consumer project shell: bash run: | consumer="${GITHUB_WORKSPACE}/ci-consumer" @@ -311,10 +327,23 @@ jobs: } EOF - cmake -S "$consumer" -B "$consumer/build" \ + - name: Configure downstream consumer (macOS) + if: runner.os == 'macOS' + run: | + BREW_PREFIX="$(brew --prefix)" + cmake -S ci-consumer -B ci-consumer/build \ -G Ninja \ -DCMAKE_BUILD_TYPE=Release \ - -DCMAKE_PREFIX_PATH="${GITHUB_WORKSPACE}/build" + -DCMAKE_PREFIX_PATH="${GITHUB_WORKSPACE}/install;${BREW_PREFIX}" + + - name: Configure downstream consumer (Windows) + if: runner.os == 'Windows' + run: > + cmake -S ci-consumer -B ci-consumer/build + -G Ninja + -DCMAKE_BUILD_TYPE=Release + -DCMAKE_PREFIX_PATH=${{ github.workspace }}/install + -DCMAKE_TOOLCHAIN_FILE=${{ github.workspace }}/vcpkg/scripts/buildsystems/vcpkg.cmake - name: Build downstream consumer - run: cmake --build ci-consumer/build --config Release \ No newline at end of file + run: cmake --build ci-consumer/build --config Release diff --git a/CMakeLists.txt b/CMakeLists.txt index 9b3a96f..1543485 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -177,6 +177,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 @@ -202,6 +203,12 @@ target_include_directories(trx PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/src ) +if(WIN32) + target_compile_definitions(trx PUBLIC NOMINMAX) + if(MSVC) + target_compile_definitions(trx PUBLIC _CRT_SECURE_NO_WARNINGS) + endif() +endif() # Fallback for system libzip packages that don't expose include dirs via # their CMake targets. Only needed for distro packages with broken configs. get_target_property(_trx_libzip_includes ${TRX_LIBZIP_TARGET} INTERFACE_INCLUDE_DIRECTORIES) @@ -223,6 +230,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) @@ -272,6 +289,12 @@ if(TRX_ENABLE_NIFTI) ${TRX_EIGEN3_TARGET} ${TRX_ZLIB_TARGET} ) + if(WIN32) + target_compile_definitions(trx-nifti PUBLIC NOMINMAX) + if(MSVC) + target_compile_definitions(trx-nifti PUBLIC _CRT_SECURE_NO_WARNINGS) + endif() + endif() endif() # ── Examples ──────────────────────────────────────────────────────────────── diff --git a/bench/CMakeLists.txt b/bench/CMakeLists.txt index 304ba4a..d5d8bc1 100644 --- a/bench/CMakeLists.txt +++ b/bench/CMakeLists.txt @@ -51,5 +51,5 @@ else() endif() add_executable(bench_trx_realdata bench_trx_realdata.cpp) -target_link_libraries(bench_trx_realdata PRIVATE trx benchmark::benchmark) +target_link_libraries(bench_trx_realdata PRIVATE trx ${TRX_LIBZIP_TARGET} benchmark::benchmark) target_compile_features(bench_trx_realdata PRIVATE cxx_std_17) diff --git a/bench/bench_trx_realdata.cpp b/bench/bench_trx_realdata.cpp index bee558d..65edb30 100644 --- a/bench/bench_trx_realdata.cpp +++ b/bench/bench_trx_realdata.cpp @@ -1,29 +1,30 @@ // Benchmark TRX streaming workloads for realistic datasets. #include #include +#include #include #include #include +#include #include #include -#include -#include -#include #include +#include #include +#include #include #include +#include #include #include -#include #include #include +#include #include #include #include #include -#include #if defined(__unix__) || defined(__APPLE__) #include @@ -64,7 +65,7 @@ enum class GroupScenario : int { None = 0, Bundles = 1, Connectome = 2 }; constexpr size_t kBundleCount = 80; constexpr std::array kConnectomeAtlasSizes = {80, 400, 1000}; -constexpr size_t kConnectomeTotalGroups = 1480; // sum of atlas sizes +constexpr size_t kConnectomeTotalGroups = 1480; // sum of atlas sizes std::string make_temp_path(const std::string &prefix) { static std::atomic counter{0}; @@ -133,8 +134,7 @@ double get_current_rss_kb() { #if defined(__APPLE__) struct mach_task_basic_info info; mach_msg_type_number_t count = MACH_TASK_BASIC_INFO_COUNT; - if (task_info(mach_task_self(), MACH_TASK_BASIC_INFO, - reinterpret_cast(&info), &count) != KERN_SUCCESS) { + if (task_info(mach_task_self(), MACH_TASK_BASIC_INFO, reinterpret_cast(&info), &count) != KERN_SUCCESS) { return 0.0; } return static_cast(info.resident_size) / 1024.0; @@ -194,17 +194,11 @@ bool is_core_profile() { return raw && std::string(raw) == "core"; } -bool include_bundles_in_core_profile() { - return parse_env_bool("TRX_BENCH_CORE_INCLUDE_BUNDLES", false); -} +bool include_bundles_in_core_profile() { return parse_env_bool("TRX_BENCH_CORE_INCLUDE_BUNDLES", false); } -size_t core_dpv_max_streamlines() { - return parse_env_size("TRX_BENCH_CORE_DPV_MAX_STREAMLINES", 1000000); -} +size_t core_dpv_max_streamlines() { return parse_env_size("TRX_BENCH_CORE_DPV_MAX_STREAMLINES", 1000000); } -size_t core_zip_max_streamlines() { - return parse_env_size("TRX_BENCH_CORE_ZIP_MAX_STREAMLINES", 1000000); -} +size_t core_zip_max_streamlines() { return parse_env_size("TRX_BENCH_CORE_ZIP_MAX_STREAMLINES", 1000000); } std::vector group_cases_for_benchmarks() { std::vector groups = {static_cast(GroupScenario::None)}; @@ -236,13 +230,13 @@ size_t group_count_for(GroupScenario scenario) { std::size_t buffer_bytes_for_streamlines(std::size_t streamlines) { std::size_t base_bytes; if (streamlines >= 5000000) { - base_bytes = 2ULL * 1024ULL * 1024ULL * 1024ULL; // 2 GB + base_bytes = 2ULL * 1024ULL * 1024ULL * 1024ULL; // 2 GB } else if (streamlines >= 1000000) { - base_bytes = 256ULL * 1024ULL * 1024ULL; // 256 MB + base_bytes = 256ULL * 1024ULL * 1024ULL; // 256 MB } else { - base_bytes = 16ULL * 1024ULL * 1024ULL; // 16 MB + base_bytes = 16ULL * 1024ULL * 1024ULL; // 16 MB } - + // Allow scaling buffer sizes for slower storage (HDD, NFS) to amortize I/O latency const size_t multiplier = std::max(1, parse_env_size("TRX_BENCH_BUFFER_MULTIPLIER", 1)); return base_bytes * multiplier; @@ -362,18 +356,16 @@ void assign_groups_to_trx(trx::TrxFile &trx, GroupScenario scenario, size_ } } -std::unique_ptr> build_prefix_subset_trx(size_t streamlines, - GroupScenario scenario, - bool add_dps, - bool add_dpv) { +std::unique_ptr> +build_prefix_subset_trx(size_t streamlines, GroupScenario scenario, bool add_dps, bool add_dpv) { if (g_reference_trx_path.empty()) { throw std::runtime_error("Reference TRX path not set."); } auto ref_trx = trx::load(g_reference_trx_path); const size_t ref_count = ref_trx->num_streamlines(); if (streamlines > ref_count) { - throw std::runtime_error("Requested " + std::to_string(streamlines) + - " streamlines but reference only has " + std::to_string(ref_count)); + throw std::runtime_error("Requested " + std::to_string(streamlines) + " streamlines but reference only has " + + std::to_string(ref_count)); } const auto ids = build_prefix_ids(streamlines); @@ -402,7 +394,8 @@ std::unique_ptr> build_prefix_subset_trx(size_t streamlines, std::vector dpv(n_verts); std::mt19937 rng(12345); std::uniform_real_distribution dist(-1.0f, 1.0f); - for (auto &v : dpv) v = dist(rng); + for (auto &v : dpv) + v = dist(rng); out->add_dpv_from_vector("dpv_random", "float32", dpv); } else { out->data_per_vertex.clear(); @@ -504,10 +497,7 @@ void register_cleanup(const std::string &path) { } } -TrxWriteStats run_trx_file_size(size_t streamlines, - bool add_dps, - bool add_dpv, - zip_uint32_t compression) { +TrxWriteStats run_trx_file_size(size_t streamlines, bool add_dps, bool add_dpv, zip_uint32_t compression) { const size_t progress_every = parse_env_size("TRX_BENCH_LOG_PROGRESS_EVERY", 0); const bool collect_rss = std::getenv("TRX_RSS_SAMPLES_PATH") != nullptr; @@ -550,7 +540,7 @@ TrxWriteStats run_trx_file_size(size_t streamlines, } trx::TrxSaveOptions save_opts; - save_opts.compression_standard = compression; + save_opts.compression = (compression == ZIP_CM_DEFLATE) ? trx::TrxCompression::Deflate : trx::TrxCompression::None; const auto start = std::chrono::steady_clock::now(); trx_subset->save(out_path, save_opts); const auto end = std::chrono::steady_clock::now(); @@ -599,9 +589,11 @@ static std::pair parse_trx_array_dims(const std::string &filenam std::istringstream ss(filename); std::string tok; while (std::getline(ss, tok, '.')) { - if (!tok.empty()) parts.push_back(tok); + if (!tok.empty()) + parts.push_back(tok); } - if (parts.size() < 2) return {1, 4}; + if (parts.size() < 2) + return {1, 4}; const std::string dtype_str = parts.back(); const size_t elem_size = static_cast(trx::detail::_sizeof_dtype(dtype_str)); if (parts.size() >= 3) { @@ -625,9 +617,11 @@ static void truncate_file_to(const std::string &path, off_t byte_size) { // Truncate every regular file in dir to row_count rows based on the per-file dtype/ncols. static void truncate_array_dir(const std::string &dir_path, size_t row_count) { std::error_code ec; - if (!trx::fs::exists(dir_path, ec)) return; + if (!trx::fs::exists(dir_path, ec)) + return; for (const auto &entry : trx::fs::directory_iterator(dir_path, ec)) { - if (ec || !entry.is_regular_file()) continue; + if (ec || !entry.is_regular_file()) + continue; const auto [ncols, elem_size] = parse_trx_array_dims(entry.path().filename().string()); truncate_file_to(entry.path().string(), static_cast(row_count * ncols * elem_size)); } @@ -644,18 +638,17 @@ static void write_synthetic_dpv_to_dir(const std::string &temp_dir, size_t n_ver if (!f.is_open()) { throw std::runtime_error("Cannot open DPV output file: " + dpv_path); } - constexpr size_t kChunkSize = 1024ULL * 1024ULL; // 1M floats = 4 MB per chunk + constexpr size_t kChunkSize = 1024ULL * 1024ULL; // 1M floats = 4 MB per chunk std::vector chunk(kChunkSize); std::mt19937 rng(12345); std::uniform_real_distribution dist(-1.0f, 1.0f); size_t remaining = n_vertices; while (remaining > 0) { - const size_t to_write = std::min(kChunkSize, remaining); + const size_t to_write = (std::min)(kChunkSize, remaining); for (size_t i = 0; i < to_write; ++i) { chunk[i] = dist(rng); } - f.write(reinterpret_cast(chunk.data()), - static_cast(to_write * sizeof(float))); + f.write(reinterpret_cast(chunk.data()), static_cast(to_write * sizeof(float))); remaining -= to_write; } } @@ -687,10 +680,13 @@ TrxOnDisk build_trx_file_on_disk_single(size_t streamlines, { std::error_code ec; for (const auto &entry : trx::fs::directory_iterator(temp_dir, ec)) { - if (ec) break; + if (ec) + break; const std::string fn = entry.path().filename().string(); - if (fn.rfind("positions", 0) == 0) pos_path = entry.path().string(); - else if (fn.rfind("offsets", 0) == 0) off_path = entry.path().string(); + if (fn.rfind("positions", 0) == 0) + pos_path = entry.path().string(); + else if (fn.rfind("offsets", 0) == 0) + off_path = entry.path().string(); } } if (pos_path.empty() || off_path.empty()) { @@ -699,8 +695,7 @@ TrxOnDisk build_trx_file_on_disk_single(size_t streamlines, // Read vertex_cutoff directly from the offsets file so the dtype (uint32 vs uint64) // is always respected, regardless of how the Eigen map interprets the mmap width. - const auto [off_ncols, off_elem] = - parse_trx_array_dims(trx::fs::path(off_path).filename().string()); + const auto [off_ncols, off_elem] = parse_trx_array_dims(trx::fs::path(off_path).filename().string()); size_t vertex_cutoff = 0; { std::ifstream ofs(off_path, std::ios::binary); @@ -722,8 +717,7 @@ TrxOnDisk build_trx_file_on_disk_single(size_t streamlines, if (!is_full_reference) { // Truncate positions and offsets to the prefix boundary. - const auto [pos_ncols, pos_elem] = - parse_trx_array_dims(trx::fs::path(pos_path).filename().string()); + const auto [pos_ncols, pos_elem] = parse_trx_array_dims(trx::fs::path(pos_path).filename().string()); truncate_file_to(pos_path, static_cast(vertex_cutoff * pos_ncols * pos_elem)); truncate_file_to(off_path, static_cast((streamlines + 1) * off_ncols * off_elem)); @@ -746,8 +740,9 @@ TrxOnDisk build_trx_file_on_disk_single(size_t streamlines, std::ifstream in(header_path); std::string raw((std::istreambuf_iterator(in)), {}); std::string parse_err; - json hdr = json::parse(raw, parse_err); - if (!parse_err.empty()) throw std::runtime_error("header.json parse error: " + parse_err); + trx::json hdr = trx::json::parse(raw, parse_err); + if (!parse_err.empty()) + throw std::runtime_error("header.json parse error: " + parse_err); hdr = trx::_json_set(hdr, "NB_STREAMLINES", static_cast(streamlines)); hdr = trx::_json_set(hdr, "NB_VERTICES", static_cast(vertex_cutoff)); std::ofstream out(header_path, std::ios::trunc); @@ -792,16 +787,16 @@ TrxOnDisk build_trx_file_on_disk_single(size_t streamlines, std::vector dpv_data(vertex_cutoff); std::mt19937 rng(12345); std::uniform_real_distribution dist(-1.0f, 1.0f); - for (auto &v : dpv_data) v = dist(rng); + for (auto &v : dpv_data) + v = dist(rng); trx->add_dpv_from_vector("dpv_random", "float32", dpv_data); } assign_groups_to_trx(*trx, scenario, streamlines); - const std::string out_path = - out_path_override.empty() ? make_temp_path("trx_input") : out_path_override; + const std::string out_path = out_path_override.empty() ? make_temp_path("trx_input") : out_path_override; trx::TrxSaveOptions save_opts; - save_opts.compression_standard = compression; + save_opts.compression = (compression == ZIP_CM_DEFLATE) ? trx::TrxCompression::Deflate : trx::TrxCompression::None; trx->save(out_path, save_opts); const size_t total_vertices = trx->num_vertices(); trx->close(); @@ -816,11 +811,8 @@ TrxOnDisk build_trx_file_on_disk_single(size_t streamlines, } } -TrxOnDisk build_trx_file_on_disk(size_t streamlines, - GroupScenario scenario, - bool add_dps, - bool add_dpv, - zip_uint32_t compression) { +TrxOnDisk build_trx_file_on_disk( + size_t streamlines, GroupScenario scenario, bool add_dps, bool add_dpv, zip_uint32_t compression) { return build_trx_file_on_disk_single(streamlines, scenario, add_dps, add_dpv, compression); } @@ -839,8 +831,8 @@ void build_slabs(std::vector> &mins, std::vector(i) / static_cast(kSlabCount - 1); const float center_z = kFov.min_z + t * z_range; - const float min_z = std::max(kFov.min_z, center_z - kSlabThicknessMm * 0.5f); - const float max_z = std::min(kFov.max_z, center_z + kSlabThicknessMm * 0.5f); + const float min_z = (std::max)(kFov.min_z, center_z - kSlabThicknessMm * 0.5f); + const float max_z = (std::min)(kFov.max_z, center_z + kSlabThicknessMm * 0.5f); mins.push_back({kFov.min_x, kFov.min_y, min_z}); maxs.push_back({kFov.max_x, kFov.max_y, max_z}); } @@ -857,9 +849,7 @@ struct KeyHash { using Key = std::tuple; size_t operator()(const Key &key) const { size_t h = 0; - auto hash_combine = [&](size_t v) { - h ^= v + 0x9e3779b97f4a7c15ULL + (h << 6) + (h >> 2); - }; + auto hash_combine = [&](size_t v) { h ^= v + 0x9e3779b97f4a7c15ULL + (h << 6) + (h >> 2); }; hash_combine(std::hash{}(std::get<0>(key))); hash_combine(std::hash{}(std::get<1>(key))); hash_combine(std::hash{}(std::get<2>(key))); @@ -871,10 +861,8 @@ struct KeyHash { void maybe_write_query_timings(const ScenarioParams &scenario, const std::vector &timings_ms) { static std::mutex mutex; static std::unordered_set seen; - const KeyHash::Key key{scenario.streamlines, - static_cast(scenario.scenario), - scenario.add_dps ? 1 : 0, - scenario.add_dpv ? 1 : 0}; + const KeyHash::Key key{ + scenario.streamlines, static_cast(scenario.scenario), scenario.add_dps ? 1 : 0, scenario.add_dpv ? 1 : 0}; std::lock_guard lock(mutex); if (!seen.insert(key).second) { @@ -932,8 +920,7 @@ static void BM_TrxFileSize_Float16(benchmark::State &state) { return; } log_bench_start("BM_TrxFileSize_Float16", - "streamlines=" + std::to_string(streamlines) + - " group_case=" + std::to_string(state.range(1)) + + "streamlines=" + std::to_string(streamlines) + " group_case=" + std::to_string(state.range(1)) + " dps=" + std::to_string(static_cast(add_dps)) + " dpv=" + std::to_string(static_cast(add_dpv)) + " compression=" + std::to_string(static_cast(use_zip))); @@ -947,10 +934,9 @@ static void BM_TrxFileSize_Float16(benchmark::State &state) { for (auto _ : state) { const double rss_iter_start = get_current_rss_kb(); const auto start = std::chrono::steady_clock::now(); - const auto on_disk = - build_trx_file_on_disk(streamlines, scenario, add_dps, add_dpv, compression); + const auto on_disk = build_trx_file_on_disk(streamlines, scenario, add_dps, add_dpv, compression); const auto end = std::chrono::steady_clock::now(); - max_rss_delta_kb = std::max(max_rss_delta_kb, get_current_rss_kb() - rss_iter_start); + max_rss_delta_kb = (std::max)(max_rss_delta_kb, get_current_rss_kb() - rss_iter_start); const std::chrono::duration elapsed = end - start; total_build_ms += elapsed.count(); total_merge_ms += on_disk.shard_merge_ms; @@ -975,17 +961,14 @@ static void BM_TrxFileSize_Float16(benchmark::State &state) { } state.counters["file_bytes"] = total_file_bytes / static_cast(state.iterations()); - log_bench_end("BM_TrxFileSize_Float16", - "streamlines=" + std::to_string(streamlines)); + log_bench_end("BM_TrxFileSize_Float16", "streamlines=" + std::to_string(streamlines)); } // Pack a directory tree into a TRX zip archive using zip_source_file for every // file so libzip reads in chunks via the OS page cache rather than mapping the // data into the process address space. Avoids the RSS spike that occurs when // TrxFile::save() syncs and re-reads large DPS/DPV mmaps before archiving. -static void pack_dir_to_zip(const std::string &src_dir, - const std::string &out_path, - zip_uint32_t compression) { +static void pack_dir_to_zip(const std::string &src_dir, const std::string &out_path, zip_uint32_t compression) { int errorp = 0; zip_t *za = zip_open(out_path.c_str(), ZIP_CREATE | ZIP_TRUNCATE, &errorp); if (!za) { @@ -1041,8 +1024,7 @@ static void BM_TrxStream_TranslateWrite(benchmark::State &state) { const Key key{streamlines, static_cast(scenario), add_dps ? 1 : 0, add_dpv ? 1 : 0}; if (cache.find(key) == cache.end()) { state.PauseTiming(); - cache.emplace(key, - build_trx_file_on_disk(streamlines, scenario, add_dps, add_dpv, ZIP_CM_STORE)); + cache.emplace(key, build_trx_file_on_disk(streamlines, scenario, add_dps, add_dpv, ZIP_CM_STORE)); state.ResumeTiming(); } @@ -1066,8 +1048,8 @@ static void BM_TrxStream_TranslateWrite(benchmark::State &state) { const long s = static_cast(get_current_rss_kb()); long prev = peak_rss_kb.load(std::memory_order_relaxed); while (s > prev && - !peak_rss_kb.compare_exchange_weak(prev, s, - std::memory_order_relaxed, std::memory_order_relaxed)) {} + !peak_rss_kb.compare_exchange_weak(prev, s, std::memory_order_relaxed, std::memory_order_relaxed)) { + } std::this_thread::sleep_for(std::chrono::milliseconds(50)); } }); @@ -1095,16 +1077,18 @@ static void BM_TrxStream_TranslateWrite(benchmark::State &state) { continue; } if (entry.path().filename().string() == pos_fname) { - continue; // will be replaced with translated positions + continue; // will be replaced with translated positions } std::filesystem::copy_file(entry.path(), std::filesystem::path(out_dir) / entry.path().filename(), - std::filesystem::copy_options::overwrite_existing, ec); + std::filesystem::copy_options::overwrite_existing, + ec); } for (const char *sub : {"dps", "dpv", "groups", "dpg"}) { const auto src = std::filesystem::path(input_dir) / sub; if (std::filesystem::exists(src, ec)) { - std::filesystem::copy(src, std::filesystem::path(out_dir) / sub, + std::filesystem::copy(src, + std::filesystem::path(out_dir) / sub, std::filesystem::copy_options::recursive | std::filesystem::copy_options::overwrite_existing, ec); @@ -1122,40 +1106,40 @@ static void BM_TrxStream_TranslateWrite(benchmark::State &state) { if (!out_positions.is_open()) { throw std::runtime_error("Failed to open output positions file: " + positions_path); } - trx.for_each_positions_chunk(chunk_bytes, - [&](trx::TrxScalarType dtype, const void *data, size_t offset, size_t count) { - (void)offset; - if (progress_every > 0 && ((offset + count) % progress_every == 0)) { - std::cerr << "[trx-bench] progress translate points=" << (offset + count) - << " / " << total_points << std::endl; - } - const size_t total_vals = count * 3; - if (dtype == trx::TrxScalarType::Float16) { - const auto *src = reinterpret_cast(data); - std::vector tmp(total_vals); - for (size_t i = 0; i < total_vals; ++i) { - tmp[i] = static_cast(static_cast(src[i]) + 1.0f); - } - out_positions.write(reinterpret_cast(tmp.data()), - static_cast(tmp.size() * sizeof(Eigen::half))); - } else if (dtype == trx::TrxScalarType::Float32) { - const auto *src = reinterpret_cast(data); - std::vector tmp(total_vals); - for (size_t i = 0; i < total_vals; ++i) { - tmp[i] = src[i] + 1.0f; - } - out_positions.write(reinterpret_cast(tmp.data()), - static_cast(tmp.size() * sizeof(float))); - } else { - const auto *src = reinterpret_cast(data); - std::vector tmp(total_vals); - for (size_t i = 0; i < total_vals; ++i) { - tmp[i] = src[i] + 1.0; - } - out_positions.write(reinterpret_cast(tmp.data()), - static_cast(tmp.size() * sizeof(double))); - } - }); + trx.for_each_positions_chunk( + chunk_bytes, [&](trx::TrxScalarType dtype, const void *data, size_t offset, size_t count) { + (void)offset; + if (progress_every > 0 && ((offset + count) % progress_every == 0)) { + std::cerr << "[trx-bench] progress translate points=" << (offset + count) << " / " << total_points + << std::endl; + } + const size_t total_vals = count * 3; + if (dtype == trx::TrxScalarType::Float16) { + const auto *src = reinterpret_cast(data); + std::vector tmp(total_vals); + for (size_t i = 0; i < total_vals; ++i) { + tmp[i] = static_cast(static_cast(src[i]) + 1.0f); + } + out_positions.write(reinterpret_cast(tmp.data()), + static_cast(tmp.size() * sizeof(Eigen::half))); + } else if (dtype == trx::TrxScalarType::Float32) { + const auto *src = reinterpret_cast(data); + std::vector tmp(total_vals); + for (size_t i = 0; i < total_vals; ++i) { + tmp[i] = src[i] + 1.0f; + } + out_positions.write(reinterpret_cast(tmp.data()), + static_cast(tmp.size() * sizeof(float))); + } else { + const auto *src = reinterpret_cast(data); + std::vector tmp(total_vals); + for (size_t i = 0; i < total_vals; ++i) { + tmp[i] = src[i] + 1.0; + } + out_positions.write(reinterpret_cast(tmp.data()), + static_cast(tmp.size() * sizeof(double))); + } + }); out_positions.flush(); out_positions.close(); @@ -1181,7 +1165,7 @@ static void BM_TrxStream_TranslateWrite(benchmark::State &state) { rss_sampling.store(false, std::memory_order_relaxed); rss_sampler.join(); const double delta = static_cast(peak_rss_kb.load(std::memory_order_relaxed)) - rss_iter_start; - max_rss_delta_kb = std::max(max_rss_delta_kb, delta); + max_rss_delta_kb = (std::max)(max_rss_delta_kb, delta); } state.counters["streamlines"] = static_cast(streamlines); @@ -1241,7 +1225,8 @@ static void BM_TrxQueryAabb_Slabs(benchmark::State &state) { const auto &min_corner = dataset.slab_mins[i]; const auto &max_corner = dataset.slab_maxs[i]; const auto q_start = std::chrono::steady_clock::now(); - auto subset = dataset.trx->query_aabb(min_corner, max_corner, + auto subset = dataset.trx->query_aabb(min_corner, + max_corner, /*precomputed_aabbs=*/nullptr, /*build_cache_for_result=*/false, max_query_streamlines, @@ -1261,11 +1246,11 @@ static void BM_TrxQueryAabb_Slabs(benchmark::State &state) { std::sort(sorted.begin(), sorted.end()); const auto p50 = sorted[sorted.size() / 2]; const auto p95_idx = static_cast(std::ceil(0.95 * sorted.size())) - 1; - const auto p95 = sorted[std::min(p95_idx, sorted.size() - 1)]; + const auto p95 = sorted[(std::min)(p95_idx, sorted.size() - 1)]; state.counters["query_p50_ms"] = p50; state.counters["query_p95_ms"] = p95; - max_rss_delta_kb = std::max(max_rss_delta_kb, get_current_rss_kb() - rss_iter_start); + max_rss_delta_kb = (std::max)(max_rss_delta_kb, get_current_rss_kb() - rss_iter_start); ScenarioParams params; params.streamlines = streamlines; @@ -1298,12 +1283,10 @@ static void ApplySizeArgs(benchmark::internal::Benchmark *bench) { const auto counts_desc = streamlines_for_benchmarks(); const auto groups = group_cases_for_benchmarks(); for (const auto count : counts_desc) { - const std::vector dpv_flags = (!core_profile || count <= dpv_limit) - ? std::vector{0, 1} - : std::vector{0}; - const std::vector compression_flags = (!core_profile || count <= zip_limit) - ? std::vector{0, 1} - : std::vector{0}; + const std::vector dpv_flags = + (!core_profile || count <= dpv_limit) ? std::vector{0, 1} : std::vector{0}; + const std::vector compression_flags = + (!core_profile || count <= zip_limit) ? std::vector{0, 1} : std::vector{0}; for (const auto group_case : groups) { for (const auto dps : flags) { for (const auto dpv : dpv_flags) { @@ -1338,9 +1321,8 @@ static void ApplyQueryArgs(benchmark::internal::Benchmark *bench) { const auto groups = group_cases_for_benchmarks(); const auto counts_desc = streamlines_for_benchmarks(); for (const auto count : counts_desc) { - const std::vector dpv_flags = (!core_profile || count <= dpv_limit) - ? std::vector{0, 1} - : std::vector{0}; + const std::vector dpv_flags = + (!core_profile || count <= dpv_limit) ? std::vector{0, 1} : std::vector{0}; for (const auto group_case : groups) { for (const auto dps : flags) { for (const auto dpv : dpv_flags) { @@ -1352,26 +1334,18 @@ static void ApplyQueryArgs(benchmark::internal::Benchmark *bench) { bench->Iterations(1); } -BENCHMARK(BM_TrxFileSize_Float16) - ->Apply(ApplySizeArgs) - ->Unit(benchmark::kMillisecond); +BENCHMARK(BM_TrxFileSize_Float16)->Apply(ApplySizeArgs)->Unit(benchmark::kMillisecond); -BENCHMARK(BM_TrxStream_TranslateWrite) - ->Apply(ApplyStreamArgs) - ->UseManualTime() - ->Unit(benchmark::kMillisecond); +BENCHMARK(BM_TrxStream_TranslateWrite)->Apply(ApplyStreamArgs)->UseManualTime()->Unit(benchmark::kMillisecond); -BENCHMARK(BM_TrxQueryAabb_Slabs) - ->Apply(ApplyQueryArgs) - ->UseManualTime() - ->Unit(benchmark::kMillisecond); +BENCHMARK(BM_TrxQueryAabb_Slabs)->Apply(ApplyQueryArgs)->UseManualTime()->Unit(benchmark::kMillisecond); int main(int argc, char **argv) { // Parse custom flags before benchmark::Initialize bool verbose = false; bool show_help = false; std::string reference_trx; - + // First pass: detect custom flags for (int i = 1; i < argc; ++i) { const std::string arg = argv[i]; @@ -1381,10 +1355,10 @@ int main(int argc, char **argv) { show_help = true; } else if (arg == "--reference-trx" && i + 1 < argc) { reference_trx = argv[i + 1]; - ++i; // Skip next arg since it's the value + ++i; // Skip next arg since it's the value } } - + if (show_help) { std::cout << "\nCustom benchmark options:\n" << " --reference-trx PATH Path to reference TRX file for sampling (REQUIRED)\n" @@ -1396,22 +1370,23 @@ int main(int argc, char **argv) { << std::endl; return 0; } - + // Validate reference TRX path if (reference_trx.empty()) { std::cerr << "Error: --reference-trx flag is required\n" << "Usage: " << argv[0] << " --reference-trx [benchmark_options]\n" - << "Use --help-custom for more information\n" << std::endl; + << "Use --help-custom for more information\n" + << std::endl; return 1; } - + // Check if reference file exists std::error_code ec; if (!std::filesystem::exists(reference_trx, ec)) { std::cerr << "Error: Reference TRX file not found: " << reference_trx << std::endl; return 1; } - + // Set global reference path g_reference_trx_path = reference_trx; std::cerr << "[trx-bench] Using reference TRX: " << g_reference_trx_path << std::endl; @@ -1425,34 +1400,34 @@ int main(int argc, char **argv) { std::cerr << "[trx-bench] Reference: " << g_reference_streamline_count << " streamlines, dpv=" << (g_reference_has_dpv ? "yes" : "no") << std::endl; } - + // Enable verbose logging if requested if (verbose) { - setenv("TRX_BENCH_LOG", "1", 0); // Don't override if already set + setenv("TRX_BENCH_LOG", "1", 0); // Don't override if already set setenv("TRX_BENCH_CHILD_LOG", "1", 0); if (std::getenv("TRX_BENCH_LOG_PROGRESS_EVERY") == nullptr) { setenv("TRX_BENCH_LOG_PROGRESS_EVERY", "50000", 1); } std::cerr << "[trx-bench] Verbose mode enabled (progress every " - << parse_env_size("TRX_BENCH_LOG_PROGRESS_EVERY", 50000) - << " streamlines)\n" << std::endl; + << parse_env_size("TRX_BENCH_LOG_PROGRESS_EVERY", 50000) << " streamlines)\n" + << std::endl; } - + // Second pass: remove custom flags from argv before passing to benchmark::Initialize - std::vector filtered_argv; - filtered_argv.push_back(argv[0]); // Keep program name + std::vector filtered_argv; + filtered_argv.push_back(argv[0]); // Keep program name for (int i = 1; i < argc; ++i) { const std::string arg = argv[i]; if (arg == "--verbose" || arg == "-v" || arg == "--help-custom") { continue; } else if (arg == "--reference-trx") { - ++i; // Skip the next arg (the path value) + ++i; // Skip the next arg (the path value) continue; } filtered_argv.push_back(argv[i]); } int filtered_argc = static_cast(filtered_argv.size()); - + ::benchmark::Initialize(&filtered_argc, filtered_argv.data()); if (::benchmark::ReportUnrecognizedArguments(filtered_argc, filtered_argv.data())) { return 1; diff --git a/cmake/trx-cppConfig.cmake.in b/cmake/trx-cppConfig.cmake.in index 85e15bb..369a41e 100644 --- a/cmake/trx-cppConfig.cmake.in +++ b/cmake/trx-cppConfig.cmake.in @@ -2,8 +2,8 @@ include(CMakeFindDependencyMacro) find_dependency(Eigen3) +find_dependency(libzip) # mio and json11 are header-only and vendored; no package required -# libzip is a PRIVATE dependency — consumers do not need it include("${CMAKE_CURRENT_LIST_DIR}/trx-cppTargets.cmake") 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..dedf8e7 100644 --- a/include/trx/trx.h +++ b/include/trx/trx.h @@ -1,6 +1,12 @@ #ifndef TRX_H // include guard #define TRX_H +#ifdef _WIN32 +# ifndef NOMINMAX +# define NOMINMAX +# endif +#endif + #include #include #include @@ -40,9 +46,8 @@ namespace trx { namespace fs = std::filesystem; -} - using json = json11::Json; +} namespace trx { enum class TrxSaveMode { Auto, Archive, Directory }; @@ -297,7 +302,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 @@ -1183,7 +1188,7 @@ TrxFile
::compute_group_connectivity(ConnectivityMeasure measure, const std:: if (b == this->group_backing_info_.end()) { continue; } - const size_t expected_ids = static_cast(std::max(0, b->second.rows)) * static_cast(std::max(0, b->second.cols)); + const size_t expected_ids = static_cast((std::max)(0, b->second.rows)) * static_cast((std::max)(0, b->second.cols)); tmp_ids.resize(expected_ids); if (expected_ids > 0) { std::ifstream in(b->second.filename, std::ios::binary); diff --git a/include/trx/trx.tpp b/include/trx/trx.tpp index ac92ef7..e9b2244 100644 --- a/include/trx/trx.tpp +++ b/include/trx/trx.tpp @@ -54,8 +54,7 @@ inline std::string folder_from_path(const std::string &elem_filename, const std: return folder; } -template -void materialize_matrix_map_and_unmap(MMappedMatrix &mapped_matrix) { +template void materialize_matrix_map_and_unmap(MMappedMatrix &mapped_matrix) { const int rows = mapped_matrix._matrix.rows(); const int cols = mapped_matrix._matrix.cols(); const size_t n = static_cast(rows) * static_cast(cols); @@ -69,8 +68,7 @@ void materialize_matrix_map_and_unmap(MMappedMatrix &mapped_matrix) { mapped_matrix.mmap.unmap(); } -template -void materialize_sequence_data_and_unmap(ArraySequence &sequence) { +template void materialize_sequence_data_and_unmap(ArraySequence &sequence) { const int rows = sequence._data.rows(); const int cols = sequence._data.cols(); const size_t n = static_cast(rows) * static_cast(cols); @@ -84,8 +82,7 @@ void materialize_sequence_data_and_unmap(ArraySequence &sequence) { sequence.mmap_pos.unmap(); } -template -void copy_cast_buffer(const InT *src, size_t n, std::vector &dst) { +template void copy_cast_buffer(const InT *src, size_t n, std::vector &dst) { dst.resize(n); for (size_t i = 0; i < n; ++i) { dst[i] = static_cast(src[i]); @@ -152,14 +149,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(); } @@ -202,8 +217,7 @@ std::string _generate_filename_from_data(const Eigen::MatrixBase
&arr, std:: return new_filename; } -template -std::unique_ptr> TrxFile
::make_empty_like() const { +template std::unique_ptr> TrxFile
::make_empty_like() const { auto empty = std::make_unique>(); empty->header = _json_set(this->header, "NB_VERTICES", 0); empty->header = _json_set(empty->header, "NB_STREAMLINES", 0); @@ -244,7 +258,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(); @@ -363,11 +377,13 @@ std::unique_ptr> _initialize_empty_trx(int nb_streamlines, int nb_ve std::tuple dpv_shape = std::make_tuple(rows, cols); trx->data_per_vertex[x.first] = std::make_unique>(); trx->data_per_vertex[x.first]->mmap_pos = trx::_create_memmap(dpv_filename, dpv_shape, "w+", dpv_dtype); - trx::detail::remap(trx->data_per_vertex[x.first]->_data, trx->data_per_vertex[x.first]->mmap_pos.data(), rows, - cols); + trx::detail::remap( + trx->data_per_vertex[x.first]->_data, trx->data_per_vertex[x.first]->mmap_pos.data(), rows, cols); - trx::detail::remap(trx->data_per_vertex[x.first]->_offsets, trx->streamlines->_offsets.data(), - int(trx->streamlines->_offsets.rows()), int(trx->streamlines->_offsets.cols())); + trx::detail::remap(trx->data_per_vertex[x.first]->_offsets, + trx->streamlines->_offsets.data(), + int(trx->streamlines->_offsets.rows()), + int(trx->streamlines->_offsets.cols())); trx->data_per_vertex[x.first]->_lengths = trx->streamlines->_lengths; } @@ -393,8 +409,8 @@ std::unique_ptr> _initialize_empty_trx(int nb_streamlines, int nb_ve trx->data_per_streamline[x.first]->mmap = trx::_create_memmap(dps_filename, dps_shape, std::string("w+"), dps_dtype); - trx::detail::remap(trx->data_per_streamline[x.first]->_matrix, trx->data_per_streamline[x.first]->mmap.data(), - rows, cols); + trx::detail::remap( + trx->data_per_streamline[x.first]->_matrix, trx->data_per_streamline[x.first]->mmap.data(), rows, cols); } } @@ -432,21 +448,19 @@ 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 == ".")) { const auto nb_vertices = static_cast(trx->header["NB_VERTICES"].int_value()); const auto expected = nb_vertices * 3; if (size != expected || dim != 3) { - throw TrxFormatError("Wrong data size/dimensionality: size=" + std::to_string(size) + - " expected=" + std::to_string(expected) + " dim=" + std::to_string(dim) + - " filename=" + elem_filename); + throw TrxFormatError("Wrong data size/dimensionality: size=" + std::to_string(size) + " expected=" + + std::to_string(expected) + " dim=" + std::to_string(dim) + " filename=" + elem_filename); } 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->streamlines->mmap_pos = trx::_create_memmap(filename, shape, "r+", ext, mem_address); trx::detail::remap(trx->streamlines->_data, trx->streamlines->mmap_pos.data(), shape); } @@ -457,16 +471,15 @@ TrxFile
::_create_trx_from_pointer(json header, const auto expected = nb_streamlines + 1; const bool missing_sentinel = (size == nb_streamlines && dim == 1); if ((size != expected && !missing_sentinel) || dim != 1) { - throw TrxFormatError("Wrong offsets size/dimensionality: size=" + std::to_string(size) + - " expected=" + std::to_string(expected) + " dim=" + std::to_string(dim) + - " filename=" + elem_filename); + throw TrxFormatError("Wrong offsets size/dimensionality: size=" + std::to_string(size) + " expected=" + + std::to_string(expected) + " dim=" + std::to_string(dim) + " filename=" + elem_filename); } const int nb_str = static_cast(trx->header["NB_STREAMLINES"].int_value()); 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); + trx->streamlines->mmap_off = + trx::_create_memmap(filename, std::make_tuple(static_cast(size), 1), "r+", ext, mem_address); if (ext == "uint64") { if (missing_sentinel) { @@ -509,17 +522,17 @@ 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); materialize_matrix_map_and_unmap(*trx->data_per_streamline[base]); } else { const size_t n = static_cast(std::get<0>(shape)) * static_cast(std::get<1>(shape)); - copy_cast_from_dtype_buffer
(trx->data_per_streamline[base]->mmap.data(), n, ext, - trx->data_per_streamline[base]->_matrix_owned); - trx::detail::remap(trx->data_per_streamline[base]->_matrix, - trx->data_per_streamline[base]->_matrix_owned.data(), shape); + copy_cast_from_dtype_buffer
( + trx->data_per_streamline[base]->mmap.data(), n, ext, trx->data_per_streamline[base]->_matrix_owned); + trx::detail::remap( + trx->data_per_streamline[base]->_matrix, trx->data_per_streamline[base]->_matrix_owned.data(), shape); trx->data_per_streamline[base]->mmap.unmap(); } } @@ -534,19 +547,21 @@ 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); } else { const size_t n = static_cast(std::get<0>(shape)) * static_cast(std::get<1>(shape)); - copy_cast_from_dtype_buffer
(trx->data_per_vertex[base]->mmap_pos.data(), n, ext, - trx->data_per_vertex[base]->_data_owned); + copy_cast_from_dtype_buffer
( + trx->data_per_vertex[base]->mmap_pos.data(), n, ext, trx->data_per_vertex[base]->_data_owned); trx::detail::remap(trx->data_per_vertex[base]->_data, trx->data_per_vertex[base]->_data_owned.data(), shape); trx->data_per_vertex[base]->mmap_pos.unmap(); } - trx::detail::remap(trx->data_per_vertex[base]->_offsets, trx->streamlines->_offsets.data(), - int(trx->streamlines->_offsets.rows()), int(trx->streamlines->_offsets.cols())); + trx::detail::remap(trx->data_per_vertex[base]->_offsets, + trx->streamlines->_offsets.data(), + int(trx->streamlines->_offsets.rows()), + int(trx->streamlines->_offsets.cols())); trx->data_per_vertex[base]->_lengths = trx->streamlines->_lengths; } @@ -564,18 +579,22 @@ 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, - trx->data_per_group[sub_folder][data_name]->mmap.data(), shape); + trx->data_per_group[sub_folder][data_name]->mmap.data(), + shape); materialize_matrix_map_and_unmap(*trx->data_per_group[sub_folder][data_name]); } else { const size_t n = static_cast(std::get<0>(shape)) * static_cast(std::get<1>(shape)); - copy_cast_from_dtype_buffer
(trx->data_per_group[sub_folder][data_name]->mmap.data(), n, ext, + copy_cast_from_dtype_buffer
(trx->data_per_group[sub_folder][data_name]->mmap.data(), + n, + ext, trx->data_per_group[sub_folder][data_name]->_matrix_owned); trx::detail::remap(trx->data_per_group[sub_folder][data_name]->_matrix, - trx->data_per_group[sub_folder][data_name]->_matrix_owned.data(), shape); + trx->data_per_group[sub_folder][data_name]->_matrix_owned.data(), + shape); trx->data_per_group[sub_folder][data_name]->mmap.unmap(); } } @@ -593,7 +612,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); @@ -708,7 +727,9 @@ template std::unique_ptr> TrxFile
::deepcopy() { copy->data_per_group[group_kv.first][field.first]->mmap = _create_memmap(dpg_filename, dpg_shape, "w+", dpg_dtype); trx::detail::remap(copy->data_per_group[group_kv.first][field.first]->_matrix, - copy->data_per_group[group_kv.first][field.first]->mmap.data(), rows, cols); + copy->data_per_group[group_kv.first][field.first]->mmap.data(), + rows, + cols); copy->data_per_group[group_kv.first][field.first]->_matrix = field.second->_matrix; } } @@ -774,7 +795,8 @@ TrxFile
::_copy_fixed_arrays_from(TrxFile
*trx, int strs_start, int pts_s this->data_per_vertex[x.first]->_data.block( pts_start, 0, curr_pts_len, this->data_per_vertex[x.first]->_data.cols()) = trx->data_per_vertex[x.first]->_data.block(0, 0, curr_pts_len, trx->data_per_vertex[x.first]->_data.cols()); - trx::detail::remap(this->data_per_vertex[x.first]->_offsets, trx->data_per_vertex[x.first]->_offsets.data(), + trx::detail::remap(this->data_per_vertex[x.first]->_offsets, + trx->data_per_vertex[x.first]->_offsets.data(), static_cast(trx->data_per_vertex[x.first]->_offsets.rows()), static_cast(trx->data_per_vertex[x.first]->_offsets.cols())); this->data_per_vertex[x.first]->_lengths = trx->data_per_vertex[x.first]->_lengths; @@ -815,8 +837,7 @@ template void TrxFile
::close() { this->header = json(header_obj); } -template -TrxFile
::~TrxFile() { +template TrxFile
::~TrxFile() { // Release mmap-backed members before deleting temporary backing directory. this->streamlines.reset(); this->groups.clear(); @@ -953,7 +974,8 @@ void TrxFile
::resize(int nb_streamlines, int nb_vertices, bool delete_dpg) { trx->data_per_group[x.first][y.first]->mmap = _create_memmap(dpg_filename, dpg_shape, "w+", dpg_dtype); trx::detail::remap(trx->data_per_group[x.first][y.first]->_matrix, - trx->data_per_group[x.first][y.first]->mmap.data(), dpg_shape); + trx->data_per_group[x.first][y.first]->mmap.data(), + dpg_shape); // update values for (int i = 0; i < trx->data_per_group[x.first][y.first]->_matrix.rows(); ++i) { @@ -1055,12 +1077,10 @@ template std::unique_ptr> TrxFile
::load(const std: return TrxFile
::load_from_zip(path); } -template std::unique_ptr> load(const std::string &path) { - return TrxFile
::load(path); -} +template std::unique_ptr> load(const std::string &path) { return TrxFile
::load(path); } -inline std::unique_ptr> -load_float32_positions(const std::string &path, const LoadFloat32Options &options) { +inline std::unique_ptr> load_float32_positions(const std::string &path, + const LoadFloat32Options &options) { const TrxScalarType dtype = detect_positions_scalar_type(path, TrxScalarType::Float32); if (dtype == TrxScalarType::Float32) { return load(path); @@ -1293,31 +1313,41 @@ template void TrxFile
::save(const std::string &filename, const converted_pos_entry = "positions.3." + new_dtype_str; } } - write_trx_archive(filename, tmp_dir_name, options.compression, - converted_pos_path, converted_pos_entry, - skip.empty() ? nullptr : &skip); + write_trx_archive(filename, + tmp_dir_name, + options.compression, + converted_pos_path, + converted_pos_entry, + skip.empty() ? nullptr : &skip); // tmp_pos_guard destructor removes the temp file after archive is written. } else { std::error_code ec; if (!trx::fs::exists(tmp_dir_name, ec) || !trx::fs::is_directory(tmp_dir_name, ec)) { throw TrxIOError("Temporary TRX directory does not exist: " + tmp_dir_name); } - if (trx::fs::exists(filename, ec) && trx::fs::is_directory(filename, ec)) { - if (!options.overwrite_existing) { - throw TrxIOError("Output directory already exists: " + filename); + trx::fs::path dest_path(filename); + std::error_code source_ec, dest_ec; + const trx::fs::path source_path = trx::fs::weakly_canonical(trx::fs::path(tmp_dir_name), source_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; + + if (!same_directory) { + if (trx::fs::exists(filename, ec) && trx::fs::is_directory(filename, ec)) { + if (!options.overwrite_existing) { + throw TrxIOError("Output directory already exists: " + filename); + } + if (rm_dir(filename) != 0) { + throw TrxIOError("Could not remove existing directory " + filename); + } } - if (rm_dir(filename) != 0) { - throw TrxIOError("Could not remove existing directory " + filename); + if (dest_path.has_parent_path()) { + mkdir_or_throw(dest_path.parent_path().string()); + } + copy_dir(tmp_dir_name, filename); + ec.clear(); + if (!trx::fs::exists(filename, ec) || !trx::fs::is_directory(filename, ec)) { + throw TrxIOError("Failed to create output directory: " + filename); } - } - trx::fs::path dest_path(filename); - if (dest_path.has_parent_path()) { - mkdir_or_throw(dest_path.parent_path().string()); - } - copy_dir(tmp_dir_name, filename); - ec.clear(); - if (!trx::fs::exists(filename, ec) || !trx::fs::is_directory(filename, ec)) { - throw TrxIOError("Failed to create output directory: " + filename); } if (options.positions_dtype.has_value() && save_trx->streamlines) { @@ -1394,7 +1424,7 @@ void TrxFile
::add_dps_from_vector(const std::string &name, const std::string if (values.size() != nb_streamlines) { throw TrxFormatError("DPS values (" + std::to_string(values.size()) + ") do not match number of streamlines (" + - std::to_string(nb_streamlines) + ")"); + std::to_string(nb_streamlines) + ")"); } std::string dps_dirname = this->_uncompressed_folder_handle + SEPARATOR + "dps" + SEPARATOR; @@ -1434,13 +1464,16 @@ void TrxFile
::add_dps_from_vector(const std::string &name, const std::string // DT memory so the in-memory matrix uses the correct element size. if (dtype_norm == "float16") { auto *ptr = reinterpret_cast(matrix->mmap.data()); - for (size_t i = 0; i < n; ++i) ptr[i] = static_cast(values[i]); + for (size_t i = 0; i < n; ++i) + ptr[i] = static_cast(values[i]); } else if (dtype_norm == "float32") { auto *ptr = reinterpret_cast(matrix->mmap.data()); - for (size_t i = 0; i < n; ++i) ptr[i] = static_cast(values[i]); + for (size_t i = 0; i < n; ++i) + ptr[i] = static_cast(values[i]); } else { auto *ptr = reinterpret_cast(matrix->mmap.data()); - for (size_t i = 0; i < n; ++i) ptr[i] = static_cast(values[i]); + for (size_t i = 0; i < n; ++i) + ptr[i] = static_cast(values[i]); } copy_cast_from_dtype_buffer
(matrix->mmap.data(), n, dtype_norm, matrix->_matrix_owned); trx::detail::remap(matrix->_matrix, matrix->_matrix_owned.data(), shape); @@ -1482,7 +1515,7 @@ void TrxFile
::add_dpv_from_vector(const std::string &name, const std::string if (values.size() != nb_vertices) { throw TrxFormatError("DPV values (" + std::to_string(values.size()) + ") do not match number of vertices (" + - std::to_string(nb_vertices) + ")"); + std::to_string(nb_vertices) + ")"); } std::string dpv_dirname = this->_uncompressed_folder_handle + SEPARATOR + "dpv" + SEPARATOR; @@ -1522,13 +1555,16 @@ void TrxFile
::add_dpv_from_vector(const std::string &name, const std::string // DT memory so the in-memory matrix uses the correct element size. if (dtype_norm == "float16") { auto *ptr = reinterpret_cast(seq->mmap_pos.data()); - for (size_t i = 0; i < n; ++i) ptr[i] = static_cast(values[i]); + for (size_t i = 0; i < n; ++i) + ptr[i] = static_cast(values[i]); } else if (dtype_norm == "float32") { auto *ptr = reinterpret_cast(seq->mmap_pos.data()); - for (size_t i = 0; i < n; ++i) ptr[i] = static_cast(values[i]); + for (size_t i = 0; i < n; ++i) + ptr[i] = static_cast(values[i]); } else { auto *ptr = reinterpret_cast(seq->mmap_pos.data()); - for (size_t i = 0; i < n; ++i) ptr[i] = static_cast(values[i]); + for (size_t i = 0; i < n; ++i) + ptr[i] = static_cast(values[i]); } copy_cast_from_dtype_buffer
(seq->mmap_pos.data(), n, dtype_norm, seq->_data_owned); trx::detail::remap(seq->_data, seq->_data_owned.data(), shape); @@ -1536,7 +1572,9 @@ void TrxFile
::add_dpv_from_vector(const std::string &name, const std::string } if (this->streamlines && this->streamlines->_offsets.size() > 0) { - trx::detail::remap(seq->_offsets, this->streamlines->_offsets.data(), int(this->streamlines->_offsets.rows()), + trx::detail::remap(seq->_offsets, + this->streamlines->_offsets.data(), + int(this->streamlines->_offsets.rows()), int(this->streamlines->_offsets.cols())); seq->_lengths = this->streamlines->_lengths; } @@ -1599,8 +1637,7 @@ void TrxFile
::add_group_from_indices(const std::string &name, const std::vec this->groups[name] = std::move(group); } -template -void TrxFile
::set_voxel_to_rasmm(const Eigen::Matrix4f &affine) { +template void TrxFile
::set_voxel_to_rasmm(const Eigen::Matrix4f &affine) { std::vector> matrix(4, std::vector(4, 0.0f)); for (int i = 0; i < 4; ++i) { for (int j = 0; j < 4; ++j) { @@ -1813,7 +1850,7 @@ TrxStream::push_dps_from_vector(const std::string &name, const std::string &dtyp tmp.reserve(chunk_elems); size_t offset = 0; while (offset < values.size()) { - const size_t count = std::min(chunk_elems, values.size() - offset); + const size_t count = (std::min)(chunk_elems, values.size() - offset); tmp.clear(); for (size_t i = 0; i < count; ++i) { tmp.push_back(static_cast(values[offset + i])); @@ -1827,7 +1864,7 @@ TrxStream::push_dps_from_vector(const std::string &name, const std::string &dtyp tmp.reserve(chunk_elems); size_t offset = 0; while (offset < values.size()) { - const size_t count = std::min(chunk_elems, values.size() - offset); + const size_t count = (std::min)(chunk_elems, values.size() - offset); tmp.clear(); for (size_t i = 0; i < count; ++i) { tmp.push_back(static_cast(values[offset + i])); @@ -1841,7 +1878,7 @@ TrxStream::push_dps_from_vector(const std::string &name, const std::string &dtyp tmp.reserve(chunk_elems); size_t offset = 0; while (offset < values.size()) { - const size_t count = std::min(chunk_elems, values.size() - offset); + const size_t count = (std::min)(chunk_elems, values.size() - offset); tmp.clear(); for (size_t i = 0; i < count; ++i) { tmp.push_back(static_cast(values[offset + i])); @@ -1892,7 +1929,7 @@ TrxStream::push_dpv_from_vector(const std::string &name, const std::string &dtyp tmp.reserve(chunk_elems); size_t offset = 0; while (offset < values.size()) { - const size_t count = std::min(chunk_elems, values.size() - offset); + const size_t count = (std::min)(chunk_elems, values.size() - offset); tmp.clear(); for (size_t i = 0; i < count; ++i) { tmp.push_back(static_cast(values[offset + i])); @@ -1906,7 +1943,7 @@ TrxStream::push_dpv_from_vector(const std::string &name, const std::string &dtyp tmp.reserve(chunk_elems); size_t offset = 0; while (offset < values.size()) { - const size_t count = std::min(chunk_elems, values.size() - offset); + const size_t count = (std::min)(chunk_elems, values.size() - offset); tmp.clear(); for (size_t i = 0; i < count; ++i) { tmp.push_back(static_cast(values[offset + i])); @@ -1920,7 +1957,7 @@ TrxStream::push_dpv_from_vector(const std::string &name, const std::string &dtyp tmp.reserve(chunk_elems); size_t offset = 0; while (offset < values.size()) { - const size_t count = std::min(chunk_elems, values.size() - offset); + const size_t count = (std::min)(chunk_elems, values.size() - offset); tmp.clear(); for (size_t i = 0; i < count; ++i) { tmp.push_back(static_cast(values[offset + i])); @@ -1976,7 +2013,7 @@ inline void TrxStream::push_group_from_indices(const std::string &name, const st const size_t chunk_elems = std::max(1, metadata_buffer_max_bytes_ / sizeof(uint32_t)); size_t offset = 0; while (offset < indices.size()) { - const size_t count = std::min(chunk_elems, indices.size() - offset); + const size_t count = (std::min)(chunk_elems, indices.size() - offset); out.write(reinterpret_cast(indices.data() + offset), static_cast(count * sizeof(uint32_t))); offset += count; @@ -2078,9 +2115,7 @@ template void TrxStream::finalize(const std::string &filename, Trx cleanup_tmp(); } -inline void TrxStream::finalize(const std::string &filename, - TrxScalarType output_dtype, - TrxCompression compression) { +inline void TrxStream::finalize(const std::string &filename, TrxScalarType output_dtype, TrxCompression compression) { switch (output_dtype) { case TrxScalarType::Float16: finalize(filename, compression); @@ -2137,7 +2172,7 @@ inline void TrxStream::finalize_directory_impl(const std::string &directory, boo trx::fs::remove_all(directory, ec); ec.clear(); } - + // Create directory if it doesn't exist if (!trx::fs::exists(directory, ec)) { mkdir_or_throw(directory); @@ -2192,7 +2227,7 @@ inline void TrxStream::finalize_directory_impl(const std::string &directory, boo tmp.reserve(chunk); size_t idx = 0; while (idx < count) { - const size_t n = std::min(chunk, count - idx); + const size_t n = (std::min)(chunk, count - idx); tmp.clear(); for (size_t i = 0; i < n; ++i) { tmp.push_back(static_cast(values.values[idx + i])); @@ -2206,7 +2241,7 @@ inline void TrxStream::finalize_directory_impl(const std::string &directory, boo tmp.reserve(chunk); size_t idx = 0; while (idx < count) { - const size_t n = std::min(chunk, count - idx); + const size_t n = (std::min)(chunk, count - idx); tmp.clear(); for (size_t i = 0; i < n; ++i) { tmp.push_back(static_cast(values.values[idx + i])); @@ -2220,7 +2255,7 @@ inline void TrxStream::finalize_directory_impl(const std::string &directory, boo tmp.reserve(chunk); size_t idx = 0; while (idx < count) { - const size_t n = std::min(chunk, count - idx); + const size_t n = (std::min)(chunk, count - idx); tmp.clear(); for (size_t i = 0; i < n; ++i) { tmp.push_back(values.values[idx + i]); @@ -2283,9 +2318,7 @@ inline void TrxStream::finalize_directory_impl(const std::string &directory, boo cleanup_tmp(); } -inline void TrxStream::finalize_directory(const std::string &directory) { - finalize_directory_impl(directory, true); -} +inline void TrxStream::finalize_directory(const std::string &directory) { finalize_directory_impl(directory, true); } inline void TrxStream::finalize_directory_persistent(const std::string &directory) { finalize_directory_impl(directory, false); @@ -2504,7 +2537,7 @@ void TrxFile
::add_dpv_from_tsf(const std::string &name, const std::string &d } if (values.size() != nb_vertices) { throw TrxFormatError("TSF values (" + std::to_string(values.size()) + ") do not match number of vertices (" + - std::to_string(nb_vertices) + ")"); + std::to_string(nb_vertices) + ")"); } std::string dpv_dirname = this->_uncompressed_folder_handle + SEPARATOR + "dpv" + SEPARATOR; @@ -2535,7 +2568,8 @@ void TrxFile
::add_dpv_from_tsf(const std::string &name, const std::string &d seq->_data(i, 0) = static_cast
(values[static_cast(i)]); } - trx::detail::remap(seq->_offsets, this->streamlines->_offsets.data(), + trx::detail::remap(seq->_offsets, + this->streamlines->_offsets.data(), static_cast(this->streamlines->_offsets.rows()), static_cast(this->streamlines->_offsets.cols())); seq->_lengths = this->streamlines->_lengths; @@ -2689,8 +2723,7 @@ template std::ostream &operator<<(std::ostream &out, const TrxFile return out; } -template -std::vector> TrxFile
::build_streamline_aabbs() const { +template std::vector> TrxFile
::build_streamline_aabbs() const { std::vector> aabbs; if (!this->streamlines) { return aabbs; @@ -2727,8 +2760,7 @@ std::vector> TrxFile
::build_streamline_aabbs() co throw TrxFormatError("Offsets exceed positions row count in build_streamline_aabbs"); } if (end <= start) { - aabbs[i] = {Eigen::half(0), Eigen::half(0), Eigen::half(0), - Eigen::half(0), Eigen::half(0), Eigen::half(0)}; + aabbs[i] = {Eigen::half(0), Eigen::half(0), Eigen::half(0), Eigen::half(0), Eigen::half(0), Eigen::half(0)}; continue; } @@ -2743,16 +2775,20 @@ std::vector> TrxFile
::build_streamline_aabbs() co const float x = static_cast(this->streamlines->_data(static_cast(p), 0)); const float y = static_cast(this->streamlines->_data(static_cast(p), 1)); const float z = static_cast(this->streamlines->_data(static_cast(p), 2)); - min_x = std::min(min_x, x); - min_y = std::min(min_y, y); - min_z = std::min(min_z, z); - max_x = std::max(max_x, x); - max_y = std::max(max_y, y); - max_z = std::max(max_z, z); + min_x = (std::min)(min_x, x); + min_y = (std::min)(min_y, y); + min_z = (std::min)(min_z, z); + max_x = (std::max)(max_x, x); + max_y = (std::max)(max_y, y); + max_z = (std::max)(max_z, z); } - aabbs[i] = {static_cast(min_x), static_cast(min_y), static_cast(min_z), - static_cast(max_x), static_cast(max_y), static_cast(max_z)}; + aabbs[i] = {static_cast(min_x), + static_cast(min_y), + static_cast(min_z), + static_cast(max_x), + static_cast(max_y), + static_cast(max_z)}; } this->aabb_cache_ = aabbs; @@ -2768,13 +2804,12 @@ const std::vector> &TrxFile
::get_or_build_streaml } template -std::unique_ptr> TrxFile
::query_aabb( - const std::array &min_corner, - const std::array &max_corner, - const std::vector> *precomputed_aabbs, - bool build_cache_for_result, - size_t max_streamlines, - uint32_t rng_seed) const { +std::unique_ptr> TrxFile
::query_aabb(const std::array &min_corner, + const std::array &max_corner, + const std::vector> *precomputed_aabbs, + bool build_cache_for_result, + size_t max_streamlines, + uint32_t rng_seed) const { if (!this->streamlines) { return this->make_empty_like(); } @@ -2785,9 +2820,10 @@ std::unique_ptr> TrxFile
::query_aabb( } std::vector> aabbs_local; - const std::vector> &aabbs = precomputed_aabbs - ? *precomputed_aabbs - : (!this->aabb_cache_.empty() ? this->aabb_cache_ : (aabbs_local = this->build_streamline_aabbs())); + const std::vector> &aabbs = + precomputed_aabbs + ? *precomputed_aabbs + : (!this->aabb_cache_.empty() ? this->aabb_cache_ : (aabbs_local = this->build_streamline_aabbs())); if (aabbs.size() != nb_streamlines) { throw TrxArgumentError("AABB size does not match streamlines count"); } @@ -2811,9 +2847,8 @@ std::unique_ptr> TrxFile
::query_aabb( const float box_max_y = static_cast(box[4]); const float box_max_z = static_cast(box[5]); - if (box_min_x <= max_x && box_max_x >= min_x && - box_min_y <= max_y && box_max_y >= min_y && - box_min_z <= max_z && box_max_z >= min_z) { + if (box_min_x <= max_x && box_max_x >= min_x && box_min_y <= max_y && box_max_y >= min_y && box_min_z <= max_z && + box_max_z >= min_z) { selected.push_back(static_cast(i)); } } @@ -2829,13 +2864,9 @@ std::unique_ptr> TrxFile
::query_aabb( return this->subset_streamlines(selected, build_cache_for_result); } -template -void TrxFile
::invalidate_aabb_cache() const { - this->aabb_cache_.clear(); -} +template void TrxFile
::invalidate_aabb_cache() const { this->aabb_cache_.clear(); } -template -const MMappedMatrix
*TrxFile
::get_dps(const std::string &name) const { +template const MMappedMatrix
*TrxFile
::get_dps(const std::string &name) const { auto it = this->data_per_streamline.find(name); if (it == this->data_per_streamline.end()) { return nullptr; @@ -2843,8 +2874,7 @@ const MMappedMatrix
*TrxFile
::get_dps(const std::string &name) const { return it->second.get(); } -template -const ArraySequence
*TrxFile
::get_dpv(const std::string &name) const { +template const ArraySequence
*TrxFile
::get_dpv(const std::string &name) const { auto it = this->data_per_vertex.find(name); if (it == this->data_per_vertex.end()) { return nullptr; @@ -2852,8 +2882,7 @@ const ArraySequence
*TrxFile
::get_dpv(const std::string &name) const { return it->second.get(); } -template -const MMappedMatrix *TrxFile
::get_group_members(const std::string &name) const { +template const MMappedMatrix *TrxFile
::get_group_members(const std::string &name) const { auto it = this->groups.find(name); if (it == this->groups.end()) { return nullptr; @@ -2867,7 +2896,7 @@ const MMappedMatrix *TrxFile
::get_group_members(const std::string const int cols = b->second.cols; std::tuple shape = std::make_tuple(rows, cols); it->second = std::make_unique>(); - const size_t n = static_cast(std::max(0, rows)) * static_cast(std::max(0, cols)); + const size_t n = static_cast((std::max)(0, rows)) * static_cast((std::max)(0, cols)); it->second->_matrix_owned.resize(n); if (n > 0) { @@ -2885,7 +2914,8 @@ const MMappedMatrix *TrxFile
::get_group_members(const std::string return nullptr; } } - in.read(reinterpret_cast(it->second->_matrix_owned.data()), static_cast(n * sizeof(uint32_t))); + in.read(reinterpret_cast(it->second->_matrix_owned.data()), + static_cast(n * sizeof(uint32_t))); if (!in) { it->second.reset(); return nullptr; @@ -2898,8 +2928,7 @@ const MMappedMatrix *TrxFile
::get_group_members(const std::string return it->second.get(); } -template -void TrxFile
::ensure_all_groups_loaded() const { +template void TrxFile
::ensure_all_groups_loaded() const { std::vector names; names.reserve(this->groups.size()); for (const auto &kv : this->groups) { @@ -2910,8 +2939,7 @@ void TrxFile
::ensure_all_groups_loaded() const { } } -template -std::vector> TrxFile
::get_streamline(size_t streamline_index) const { +template std::vector> TrxFile
::get_streamline(size_t streamline_index) const { if (!this->streamlines || this->streamlines->_offsets.size() == 0) { throw TrxFormatError("TRX streamlines are not available"); } @@ -2920,7 +2948,8 @@ std::vector> TrxFile
::get_streamline(size_t streamline_ind throw std::out_of_range("Streamline index out of range"); } - const uint64_t start = static_cast(this->streamlines->_offsets(static_cast(streamline_index), 0)); + const uint64_t start = + static_cast(this->streamlines->_offsets(static_cast(streamline_index), 0)); const uint64_t end = static_cast(this->streamlines->_offsets(static_cast(streamline_index + 1), 0)); std::vector> points; @@ -2936,9 +2965,7 @@ std::vector> TrxFile
::get_streamline(size_t streamline_ind return points; } -template -template -void TrxFile
::for_each_streamline(Fn &&fn) const { +template template void TrxFile
::for_each_streamline(Fn &&fn) const { if (!this->streamlines || this->streamlines->_offsets.size() == 0) { return; } @@ -3035,8 +3062,7 @@ void TrxFile
::add_dpg_from_matrix(const std::string &group, values.push_back(matrix(i, j)); } } - add_dpg_from_vector(group, name, dtype, values, static_cast(matrix.rows()), - static_cast(matrix.cols())); + add_dpg_from_vector(group, name, dtype, values, static_cast(matrix.rows()), static_cast(matrix.cols())); } template @@ -3052,8 +3078,7 @@ const MMappedMatrix
*TrxFile
::get_dpg(const std::string &group, const st return field_it->second.get(); } -template -std::vector TrxFile
::list_dpg_groups() const { +template std::vector TrxFile
::list_dpg_groups() const { std::vector groups; groups.reserve(this->data_per_group.size()); for (const auto &kv : this->data_per_group) { @@ -3062,8 +3087,7 @@ std::vector TrxFile
::list_dpg_groups() const { return groups; } -template -std::vector TrxFile
::list_dpg_fields(const std::string &group) const { +template std::vector TrxFile
::list_dpg_fields(const std::string &group) const { std::vector fields; auto it = this->data_per_group.find(group); if (it == this->data_per_group.end()) { @@ -3076,8 +3100,7 @@ std::vector TrxFile
::list_dpg_fields(const std::string &group) return fields; } -template -void TrxFile
::remove_dpg(const std::string &group, const std::string &name) { +template void TrxFile
::remove_dpg(const std::string &group, const std::string &name) { auto group_it = this->data_per_group.find(group); if (group_it == this->data_per_group.end()) { return; @@ -3088,8 +3111,7 @@ void TrxFile
::remove_dpg(const std::string &group, const std::string &name) } } -template -void TrxFile
::remove_dpg_group(const std::string &group) { +template void TrxFile
::remove_dpg_group(const std::string &group) { this->data_per_group.erase(group); } @@ -3149,9 +3171,7 @@ std::unique_ptr> TrxFile
::subset_streamlines(const std::vector(end - start); } - auto out = std::make_unique>(static_cast(total_vertices), - static_cast(selected.size()), - this); + auto out = std::make_unique>(static_cast(total_vertices), static_cast(selected.size()), this); out->header = _json_set(this->header, "NB_VERTICES", static_cast(total_vertices)); out->header = _json_set(out->header, "NB_STREAMLINES", static_cast(selected.size())); @@ -3168,14 +3188,11 @@ std::unique_ptr> TrxFile
::subset_streamlines(const std::vector(new_idx)) = static_cast(len); - out_offsets(static_cast(new_idx + 1), 0) = - out_offsets(static_cast(new_idx), 0) + len; + out_offsets(static_cast(new_idx + 1), 0) = out_offsets(static_cast(new_idx), 0) + len; if (len > 0) { - out_positions.block(static_cast(cursor), 0, - static_cast(len), 3) = - this->streamlines->_data.block(static_cast(start), 0, - static_cast(len), 3); + out_positions.block(static_cast(cursor), 0, static_cast(len), 3) = + this->streamlines->_data.block(static_cast(start), 0, static_cast(len), 3); for (const auto &kv : this->data_per_vertex) { const std::string &name = kv.first; @@ -3186,10 +3203,8 @@ std::unique_ptr> TrxFile
::subset_streamlines(const std::vectorsecond->_data; auto &src_dpv = kv.second->_data; const Eigen::Index cols = src_dpv.cols(); - out_dpv.block(static_cast(cursor), 0, - static_cast(len), cols) = - src_dpv.block(static_cast(start), 0, - static_cast(len), cols); + out_dpv.block(static_cast(cursor), 0, static_cast(len), cols) = + src_dpv.block(static_cast(start), 0, static_cast(len), cols); } } @@ -3258,21 +3273,20 @@ std::unique_ptr> TrxFile
::subset_streamlines(const std::vector_matrix, dpg_filename); - std::tuple dpg_shape = std::make_tuple(field_kv.second->_matrix.rows(), - field_kv.second->_matrix.cols()); + std::tuple dpg_shape = + std::make_tuple(field_kv.second->_matrix.rows(), field_kv.second->_matrix.cols()); out->data_per_group[group_name][field_name] = std::make_unique>(); - out->data_per_group[group_name][field_name]->mmap = - _create_memmap(dpg_filename, dpg_shape, "w+", dpg_dtype); + out->data_per_group[group_name][field_name]->mmap = _create_memmap(dpg_filename, dpg_shape, "w+", dpg_dtype); trx::detail::remap(out->data_per_group[group_name][field_name]->_matrix, out->data_per_group[group_name][field_name]->mmap.data(), - std::get<0>(dpg_shape), std::get<1>(dpg_shape)); + std::get<0>(dpg_shape), + std::get<1>(dpg_shape)); for (int i = 0; i < out->data_per_group[group_name][field_name]->_matrix.rows(); ++i) { for (int j = 0; j < out->data_per_group[group_name][field_name]->_matrix.cols(); ++j) { - out->data_per_group[group_name][field_name]->_matrix(i, j) = - field_kv.second->_matrix(i, j); + out->data_per_group[group_name][field_name]->_matrix(i, j) = field_kv.second->_matrix(i, j); } } } 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..2b2773c --- /dev/null +++ b/src/legacy_io.cpp @@ -0,0 +1,1095 @@ +#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) { + try { + 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)); + if (!f) + return false; + 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)); + if (!f) + return false; + 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; + + const std::streampos offsets_start = f.tellg(); + f.seekg(0, std::ios::end); + const size_t off_bytes_avail = static_cast(f.tellg() - offsets_start); + f.seekg(offsets_start); + const size_t offset_elem_size = is_int64 ? sizeof(uint64_t) : sizeof(uint32_t); + if (num_offsets > off_bytes_avail / offset_elem_size) + return false; + + tr.offsets.resize(num_offsets); + for (size_t i = 0; i < num_offsets; ++i) { + if (is_int64) { + uint64_t val = 0; + f.read(reinterpret_cast(&val), 8); + if (!f) + return false; + 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 = 0; + f.read(reinterpret_cast(&val), 4); + if (!f) + return false; + val = swap_int32(val); + tr.offsets[i] = val; + } + } + return true; + } + f.seekg(pos_before_offsets); + + const std::streampos lines_start = f.tellg(); + f.seekg(0, std::ios::end); + const size_t lines_bytes_avail = static_cast(f.tellg() - lines_start); + f.seekg(lines_start); + + 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 = 0; + f.read(reinterpret_cast(&n_pts), sizeof(int32_t)); + if (!f) + return false; + n_pts = swap_int32(n_pts); + if (n_pts < 0) + return false; + if (n_pts == 0) + continue; + if (static_cast(n_pts) * sizeof(int32_t) > lines_bytes_avail) + return false; + 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)); + if (!f) + return false; + } + + if (tr.offsets.empty() || tr.offsets.size() != num_streamlines + 1) + return false; + + return true; + } catch (const std::exception &) { + return false; + } +} + +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]; + } + + // Preserve the vertex/streamline counts derived from the data. Assigning + // header_to_use directly would drop NB_VERTICES / NB_STREAMLINES for inputs + // whose header only carries DIMENSIONS / VOXEL_TO_RASMM (TRK/TCK/VTK), + // producing a TRX that fails to load. + auto header_obj = header_to_use.is_object() ? header_to_use.object_items() : std::map(); + header_obj["NB_VERTICES"] = static_cast(nb_vertices); + header_obj["NB_STREAMLINES"] = static_cast(nb_streamlines); + trx.header = header_obj; + + 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..59d8823 100644 --- a/src/trx.cpp +++ b/src/trx.cpp @@ -17,10 +17,10 @@ #include #include #include -#include #include #include #include +#include #include #include #include @@ -31,8 +31,8 @@ #endif #include -#include #include +#include // #define ZIP_DD_SIG 0x08074b50 // #define ZIP_CD_SIG 0x06054b50 @@ -52,8 +52,11 @@ namespace trx { // Forward declarations for functions defined later in this file. // These were previously declared in trx.h but are now internal. std::string extract_zip_to_directory(zip_t *zfolder); -void zip_from_folder(zip_t *zf, const std::string &root, const std::string &directory, - zip_uint32_t compression_standard, const std::unordered_set *skip); +void zip_from_folder(zip_t *zf, + const std::string &root, + const std::string &directory, + zip_uint32_t compression_standard, + const std::unordered_set *skip); json load_header(zip_t *zfolder); zip_uint32_t to_zip_compression(TrxCompression c) { @@ -104,7 +107,6 @@ std::string normalize_slashes(std::string path) { return path; } - bool parse_positions_dtype(const std::string &filename, std::string &out_dtype) { const std::string normalized = normalize_slashes(filename); try { @@ -188,6 +190,96 @@ 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,232 @@ 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") { + auto arr = read_entry_to_typed_array(static_cast(count_elems), 1); + arr.materialize_to_owned(); + const uint64_t nb_streamlines_u64 = static_cast(trx.header["NB_STREAMLINES"].number_value()); + const auto *vals = reinterpret_cast(arr.owned.data()); + for (size_t idx = 0; idx < count_elems; ++idx) { + if (static_cast(vals[idx]) >= nb_streamlines_u64) { + throw TrxFormatError("Group '" + base + "' contains a streamline index >= NB_STREAMLINES"); + } + } + trx.groups.emplace(base, std::move(arr)); + } 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 +825,89 @@ 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(); + const uint64_t nb_streamlines_u64 = static_cast(header["NB_STREAMLINES"].number_value()); + const auto *vals = reinterpret_cast(arr.owned.data()); + for (size_t idx = 0; idx < static_cast(size); ++idx) { + if (static_cast(vals[idx]) >= nb_streamlines_u64) { + throw TrxFormatError("Group '" + base + "' contains a streamline index >= NB_STREAMLINES"); + } + } + 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 +944,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,18 +1001,76 @@ 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) { + if (arr.empty()) { + return; + } + // If the array is purely memory-mapped and destination already exists, + // sync dirty pages to disk instead of truncating the file from under the active mmap + // (truncating an mmapped file leads to SIGBUS on subsequent reads). + std::error_code ec; + if (arr.owned.empty() && arr.mmap.is_open() && trx::fs::exists(path, ec)) { + const_cast(arr.mmap).sync(ec); + return; + } + + const auto bytes = arr.to_bytes(); + const std::string tmp_path = path + ".tmp"; + { + std::ofstream out(tmp_path, std::ios::binary | std::ios::out | std::ios::trunc); + if (!out.is_open()) { + throw TrxIOError("Failed to open output file: " + tmp_path); + } + if (bytes.data && bytes.size > 0) { + out.write(reinterpret_cast(bytes.data), static_cast(bytes.size)); + } + out.flush(); + } + trx::fs::rename(tmp_path, path, ec); + if (ec) { + trx::fs::copy_file(tmp_path, path, trx::fs::copy_options::overwrite_existing, ec); + trx::fs::remove(tmp_path, ec); + } +} +} // namespace + void AnyTrxFile::save(const std::string &filename, TrxCompression compression) { TrxSaveOptions options; options.compression = compression; save(filename, options); } -using trx::detail::TempFileGuard; using trx::detail::make_unique_temp_path; +using trx::detail::TempFileGuard; void AnyTrxFile::save(const std::string &filename, const TrxSaveOptions &options) { const std::string ext = get_ext(filename); @@ -636,128 +1079,162 @@ 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; + + 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()); - std::unordered_set skip = {"header.json"}; - // Guard deletes the temp positions file after commit (or on exception). - TempFileGuard tmp_pos_guard; + // 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); - zf.commit(filename); - } else { - std::error_code ec; - if (trx::fs::exists(filename, ec) && trx::fs::is_directory(filename, ec)) { - if (!options.overwrite_existing) { - throw TrxIOError("Output directory already exists: " + filename); + 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); } - if (rm_dir(filename) != 0) { - throw TrxIOError("Could not remove existing directory " + filename); + } + + // 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); } } - trx::fs::path dest_path(filename); - if (dest_path.has_parent_path()) { - mkdir_or_throw(dest_path.parent_path().string()); + + // 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); + } } - 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; - if (!same_directory) { - copy_dir(source_dir, filename); + // 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); + } + } } - 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); + zf.commit(filename); + } else { + // TrxSaveMode::Directory + std::error_code ec; + trx::fs::path dest_path(filename); + std::error_code source_ec, dest_ec; + const trx::fs::path source_path = _backing_directory.empty() + ? trx::fs::path() + : trx::fs::weakly_canonical(trx::fs::path(_backing_directory), source_ec); + const trx::fs::path normalized_dest = trx::fs::weakly_canonical(dest_path, dest_ec); + const bool same_directory = !_backing_directory.empty() && !source_ec && !dest_ec && source_path == normalized_dest; + + if (!same_directory) { + if (trx::fs::exists(filename, ec) && trx::fs::is_directory(filename, ec)) { + if (!options.overwrite_existing) { + throw TrxIOError("Output directory already exists: " + filename); + } + if (rm_dir(filename) != 0) { + throw TrxIOError("Could not remove existing directory " + filename); + } } + if (dest_path.has_parent_path()) { + mkdir_or_throw(dest_path.parent_path().string()); + } + mkdir_or_throw(filename); } const trx::fs::path final_header_path = dest_path / "header.json"; @@ -768,12 +1245,62 @@ void AnyTrxFile::save(const std::string &filename, const TrxSaveOptions &options 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 (options.positions_dtype.has_value() && !positions.empty()) { + const TrxScalarType target = *options.positions_dtype; + const std::string new_dtype_str = scalar_type_name(target); + 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); + } + + if (!offsets.empty()) { + const std::string off_path = filename + SEPARATOR + typed_array_filename("offsets", offsets); + write_typed_array_file(off_path, offsets); + } + + 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 (!trx::fs::exists(final_header_path)) { - throw TrxFormatError("Missing header.json in output directory: " + final_header_path.string()); + + if (!data_per_streamline.empty()) { + const std::string dps_dir = filename + SEPARATOR + "dps"; + trx::fs::create_directories(dps_dir, ec); + for (const auto &kv : data_per_streamline) { + write_typed_array_file(dps_dir + SEPARATOR + typed_array_filename(kv.first, kv.second), kv.second); + } + } + + if (!data_per_vertex.empty()) { + const std::string dpv_dir = filename + SEPARATOR + "dpv"; + trx::fs::create_directories(dpv_dir, ec); + for (const auto &kv : data_per_vertex) { + write_typed_array_file(dpv_dir + SEPARATOR + typed_array_filename(kv.first, kv.second), kv.second); + } + } + + if (!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 +1428,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'; } } @@ -916,7 +1443,7 @@ mio::shared_mmap_sink _create_memmap(std::string filename, static_cast(trx::detail::_sizeof_dtype(dtype)); // if file does not exist, create and allocate it - struct stat buffer {}; + struct stat buffer{}; if (stat(filename.c_str(), &buffer) != 0) { allocate_file(filename, filesize); } @@ -1068,8 +1595,7 @@ std::string make_temp_dir(const std::string &prefix) { static_cast(getpid()); #endif for (int attempt = 0; attempt < 100; ++attempt) { - const trx::fs::path candidate = - base_path / (prefix + "_" + std::to_string(pid) + "_" + std::to_string(dist(rng))); + const trx::fs::path candidate = base_path / (prefix + "_" + std::to_string(pid) + "_" + std::to_string(dist(rng))); ec.clear(); if (trx::fs::create_directory(candidate, ec)) { return candidate.string(); @@ -1199,11 +1725,11 @@ std::string extract_trx_archive(const std::string &zip_path) { } void write_trx_archive(const std::string &filename, - const std::string &source_dir, - TrxCompression compression, - const std::string &converted_positions_path, - const std::string &converted_positions_entry, - const std::unordered_set *skip) { + const std::string &source_dir, + TrxCompression compression, + const std::string &converted_positions_path, + const std::string &converted_positions_entry, + const std::unordered_set *skip) { const zip_uint32_t zip_comp = to_zip_compression(compression); int errorp; detail::ZipArchive zf(zip_open(filename.c_str(), ZIP_CREATE + ZIP_TRUNCATE, &errorp)); @@ -1215,8 +1741,8 @@ void write_trx_archive(const std::string &filename, zip_source_t *pos_src = zip_source_file(zf.get(), converted_positions_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(), converted_positions_entry.c_str(), pos_src, ZIP_FL_ENC_UTF_8 | ZIP_FL_OVERWRITE); + const zip_int64_t pos_idx = + zip_file_add(zf.get(), converted_positions_entry.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"); if (zip_set_file_compression(zf.get(), pos_idx, static_cast(zip_comp), 0) < 0) @@ -1251,26 +1777,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 { @@ -1296,7 +1802,7 @@ void AnyTrxFile::for_each_positions_chunk(size_t chunk_bytes, const PositionsChu const auto *base = bytes.data; const auto dtype = scalar_type_from_dtype(positions.dtype); for (size_t offset = 0; offset < total_points; offset += points_per_chunk) { - const size_t count = std::min(points_per_chunk, total_points - offset); + const size_t count = (std::min)(points_per_chunk, total_points - offset); const void *ptr = base + offset * bytes_per_point; fn(dtype, ptr, offset, count); } @@ -1325,7 +1831,7 @@ void AnyTrxFile::for_each_positions_chunk_mutable(size_t chunk_bytes, const Posi auto *base = bytes.data; const auto dtype = scalar_type_from_dtype(positions.dtype); for (size_t offset = 0; offset < total_points; offset += points_per_chunk) { - const size_t count = std::min(points_per_chunk, total_points - offset); + const size_t count = (std::min)(points_per_chunk, total_points - offset); void *ptr = base + offset * bytes_per_point; fn(dtype, ptr, offset, count); } @@ -1463,46 +1969,47 @@ void merge_trx_shards(const MergeTrxShardsOptions &options) { } }; - auto append_offsets_with_base = [](const std::string &dst, const std::string &src, uint64_t base_vertices, bool skip_first) { - std::ifstream in(src, std::ios::binary); - if (!in.is_open()) { - throw TrxIOError("Failed to open source offsets: " + src); - } - std::ofstream out(dst, std::ios::binary | std::ios::app); - if (!out.is_open()) { - throw TrxIOError("Failed to open destination offsets: " + dst); - } - constexpr size_t kChunkElems = (8 * 1024 * 1024) / sizeof(uint64_t); - std::vector buffer(kChunkElems); - bool first_value_pending = skip_first; - while (in) { - in.read(reinterpret_cast(buffer.data()), - static_cast(buffer.size() * sizeof(uint64_t))); - const std::streamsize bytes = in.gcount(); - if (bytes <= 0) { - break; - } - if (bytes % static_cast(sizeof(uint64_t)) != 0) { - throw TrxFormatError("Offsets file has invalid byte count: " + src); - } - const size_t count = static_cast(bytes) / sizeof(uint64_t); - size_t start_index = 0; - if (first_value_pending) { - if (count == 0) { - continue; + auto append_offsets_with_base = + [](const std::string &dst, const std::string &src, uint64_t base_vertices, bool skip_first) { + std::ifstream in(src, std::ios::binary); + if (!in.is_open()) { + throw TrxIOError("Failed to open source offsets: " + src); } - start_index = 1; - first_value_pending = false; - } - for (size_t i = start_index; i < count; ++i) { - buffer[i] += base_vertices; - } - if (count > start_index) { - out.write(reinterpret_cast(buffer.data() + start_index), - static_cast((count - start_index) * sizeof(uint64_t))); - } - } - }; + std::ofstream out(dst, std::ios::binary | std::ios::app); + if (!out.is_open()) { + throw TrxIOError("Failed to open destination offsets: " + dst); + } + constexpr size_t kChunkElems = (8 * 1024 * 1024) / sizeof(uint64_t); + std::vector buffer(kChunkElems); + bool first_value_pending = skip_first; + while (in) { + in.read(reinterpret_cast(buffer.data()), + static_cast(buffer.size() * sizeof(uint64_t))); + const std::streamsize bytes = in.gcount(); + if (bytes <= 0) { + break; + } + if (bytes % static_cast(sizeof(uint64_t)) != 0) { + throw TrxFormatError("Offsets file has invalid byte count: " + src); + } + const size_t count = static_cast(bytes) / sizeof(uint64_t); + size_t start_index = 0; + if (first_value_pending) { + if (count == 0) { + continue; + } + start_index = 1; + first_value_pending = false; + } + for (size_t i = start_index; i < count; ++i) { + buffer[i] += base_vertices; + } + if (count > start_index) { + out.write(reinterpret_cast(buffer.data() + start_index), + static_cast((count - start_index) * sizeof(uint64_t))); + } + } + }; auto append_group_indices_with_base = [](const std::string &dst, const std::string &src, uint32_t base_streamlines) { std::ifstream in(src, std::ios::binary); @@ -1516,8 +2023,7 @@ void merge_trx_shards(const MergeTrxShardsOptions &options) { constexpr size_t kChunkElems = (8 * 1024 * 1024) / sizeof(uint32_t); std::vector buffer(kChunkElems); while (in) { - in.read(reinterpret_cast(buffer.data()), - static_cast(buffer.size() * sizeof(uint32_t))); + in.read(reinterpret_cast(buffer.data()), static_cast(buffer.size() * sizeof(uint32_t))); const std::streamsize bytes = in.gcount(); if (bytes <= 0) { break; @@ -1556,12 +2062,13 @@ void merge_trx_shards(const MergeTrxShardsOptions &options) { return files; }; - auto ensure_schema_match = [&](const std::string &subdir, const std::vector &schema_files, const std::string &shard) { - const auto shard_files = list_subdir_files(shard, subdir); - if (shard_files != schema_files) { - throw TrxFormatError("Shard schema mismatch for subdir '" + subdir + "': " + shard); - } - }; + auto ensure_schema_match = + [&](const std::string &subdir, const std::vector &schema_files, const std::string &shard) { + const auto shard_files = list_subdir_files(shard, subdir); + if (shard_files != schema_files) { + throw TrxFormatError("Shard schema mismatch for subdir '" + subdir + "': " + shard); + } + }; std::error_code ec; for (const auto &dir : options.shard_directories) { @@ -1625,13 +2132,15 @@ void merge_trx_shards(const MergeTrxShardsOptions &options) { trx::fs::create_directories(output_dir + SEPARATOR + "groups", ec); } for (const auto &name : dps_schema) { - std::ofstream clear_file(output_dir + SEPARATOR + "dps" + SEPARATOR + name, std::ios::binary | std::ios::out | std::ios::trunc); + std::ofstream clear_file(output_dir + SEPARATOR + "dps" + SEPARATOR + name, + std::ios::binary | std::ios::out | std::ios::trunc); if (!clear_file.is_open()) { throw TrxIOError("Failed to create merged dps file: " + name); } } for (const auto &name : dpv_schema) { - std::ofstream clear_file(output_dir + SEPARATOR + "dpv" + SEPARATOR + name, std::ios::binary | std::ios::out | std::ios::trunc); + std::ofstream clear_file(output_dir + SEPARATOR + "dpv" + SEPARATOR + name, + std::ios::binary | std::ios::out | std::ios::trunc); if (!clear_file.is_open()) { throw TrxIOError("Failed to create merged dpv file: " + name); } @@ -1672,19 +2181,20 @@ void merge_trx_shards(const MergeTrxShardsOptions &options) { append_offsets_with_base(offsets_out, shard_offsets, total_vertices, i != 0); for (const auto &name : dps_schema) { - append_binary(output_dir + SEPARATOR + "dps" + SEPARATOR + name, shard_dir + SEPARATOR + "dps" + SEPARATOR + name); + append_binary(output_dir + SEPARATOR + "dps" + SEPARATOR + name, + shard_dir + SEPARATOR + "dps" + SEPARATOR + name); } for (const auto &name : dpv_schema) { - append_binary(output_dir + SEPARATOR + "dpv" + SEPARATOR + name, shard_dir + SEPARATOR + "dpv" + SEPARATOR + name); + append_binary(output_dir + SEPARATOR + "dpv" + SEPARATOR + name, + shard_dir + SEPARATOR + "dpv" + SEPARATOR + name); } for (const auto &name : groups_schema) { if (total_streamlines > static_cast(std::numeric_limits::max())) { throw TrxFormatError("Group index offset exceeds uint32 range during merge"); } - append_group_indices_with_base( - output_dir + SEPARATOR + "groups" + SEPARATOR + name, - shard_dir + SEPARATOR + "groups" + SEPARATOR + name, - static_cast(total_streamlines)); + append_group_indices_with_base(output_dir + SEPARATOR + "groups" + SEPARATOR + name, + shard_dir + SEPARATOR + "groups" + SEPARATOR + name, + static_cast(total_streamlines)); } total_vertices += shard_vertices; @@ -1732,8 +2242,12 @@ struct RawEntry { // Adds one entry to an already-open zip archive. The data is copied into a // malloc buffer that libzip takes ownership of (free=1). // If overwrite=false and the entry already exists the function is a no-op. -void zip_add_buffer_entry(zip_t *zf, const std::string &entry, const void *data, - std::size_t nbytes, zip_uint32_t compression, bool overwrite) { +void zip_add_buffer_entry(zip_t *zf, + const std::string &entry, + const void *data, + std::size_t nbytes, + zip_uint32_t compression, + bool overwrite) { std::string normalized_entry = entry; std::replace(normalized_entry.begin(), normalized_entry.end(), '\\', '/'); if (!overwrite) { @@ -1772,9 +2286,11 @@ void zip_add_buffer_entry(zip_t *zf, const std::string &entry, const void *data, // Opens the zip at `path` without truncating it, adds a directory entry for // `subdir` (harmless if already present), writes each RawEntry under that // subdir, then commits. If overwrite=false, existing entries are skipped. -void append_raw_entries_to_zip(const std::string &path, const std::string &subdir, - const std::vector &entries, zip_uint32_t compression, - bool overwrite) { +void append_raw_entries_to_zip(const std::string &path, + const std::string &subdir, + const std::vector &entries, + zip_uint32_t compression, + bool overwrite) { if (entries.empty()) { return; } @@ -1785,16 +2301,17 @@ void append_raw_entries_to_zip(const std::string &path, const std::string &subdi } zip_dir_add(zf.get(), subdir.c_str(), ZIP_FL_ENC_UTF_8); for (const auto &e : entries) { - zip_add_buffer_entry(zf.get(), subdir + "/" + e.filename, e.data, e.nbytes, compression, - overwrite); + zip_add_buffer_entry(zf.get(), subdir + "/" + e.filename, e.data, e.nbytes, compression, overwrite); } zf.commit(path); } // Creates `directory/subdir/` if absent, then writes each RawEntry as a // binary file. If overwrite=false, existing files are skipped. -void append_raw_entries_to_directory(const std::string &directory, const std::string &subdir, - const std::vector &entries, bool overwrite) { +void append_raw_entries_to_directory(const std::string &directory, + const std::string &subdir, + const std::vector &entries, + bool overwrite) { if (entries.empty()) { return; } @@ -1822,7 +2339,8 @@ void append_raw_entries_to_directory(const std::string &directory, const std::st void append_groups_to_zip(const std::string &path, const std::map> &groups, - TrxCompression compression, bool overwrite) { + TrxCompression compression, + bool overwrite) { std::vector entries; entries.reserve(groups.size()); for (const auto &kv : groups) { @@ -1842,8 +2360,10 @@ void append_groups_to_directory(const std::string &directory, append_raw_entries_to_directory(directory, "groups", entries, overwrite); } -void append_dps_to_zip(const std::string &path, const std::map &dps, - TrxCompression compression, bool overwrite) { +void append_dps_to_zip(const std::string &path, + const std::map &dps, + TrxCompression compression, + bool overwrite) { std::vector entries; entries.reserve(dps.size()); for (const auto &kv : dps) { @@ -1854,7 +2374,8 @@ void append_dps_to_zip(const std::string &path, const std::map &dps, bool overwrite) { + const std::map &dps, + bool overwrite) { std::vector entries; entries.reserve(dps.size()); for (const auto &kv : dps) { @@ -1864,8 +2385,10 @@ void append_dps_to_directory(const std::string &directory, append_raw_entries_to_directory(directory, "dps", entries, overwrite); } -void append_dpv_to_zip(const std::string &path, const std::map &dpv, - TrxCompression compression, bool overwrite) { +void append_dpv_to_zip(const std::string &path, + const std::map &dpv, + TrxCompression compression, + bool overwrite) { std::vector entries; entries.reserve(dpv.size()); for (const auto &kv : dpv) { @@ -1876,7 +2399,8 @@ void append_dpv_to_zip(const std::string &path, const std::map &dpv, bool overwrite) { + const std::map &dpv, + bool overwrite) { std::vector entries; entries.reserve(dpv.size()); for (const auto &kv : dpv) { @@ -1897,8 +2421,8 @@ static std::string make_prefix_key(const std::string &name, int depth) { return name.substr(0, pos - 1); } -std::string format_groups_summary(const std::map &groups, int prefix_depth, - const std::string &line_prefix) { +std::string +format_groups_summary(const std::map &groups, int prefix_depth, const std::string &line_prefix) { if (groups.empty()) return ""; std::ostringstream out; diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index f37f3db..89d5944 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,16 @@ 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) + +add_executable(test_legacy_io test_trx_legacy_io.cpp) +target_link_libraries(test_legacy_io PRIVATE trx GTest::gtest_main) +target_compile_features(test_legacy_io PRIVATE cxx_std_17) + include(GoogleTest) gtest_discover_tests(test_mmap PROPERTIES ENVIRONMENT "TRX_TEST_DATA_DIR=${TRX_TEST_DATA_DIR}" @@ -124,3 +140,11 @@ 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}" +) + +gtest_discover_tests(test_legacy_io 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_legacy_io.cpp b/tests/test_trx_legacy_io.cpp new file mode 100644 index 0000000..89aee0d --- /dev/null +++ b/tests/test_trx_legacy_io.cpp @@ -0,0 +1,258 @@ +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +namespace fs = std::filesystem; + +namespace { +std::string test_data_root() { + const auto *env = std::getenv("TRX_TEST_DATA_DIR"); // NOLINT(concurrency-mt-unsafe) + if (env != nullptr && !std::string(env).empty()) { + fs::path dir = fs::path(env) / "gs"; + if (fs::exists(dir / "gs.trx")) + return dir.string(); + if (fs::exists(fs::path(env) / "gs.trx")) + return std::string(env); + } + fs::path repo_data = fs::path(__FILE__).parent_path() / "test_data" / "gs"; + if (fs::exists(repo_data / "gs.trx")) + return repo_data.string(); + return {}; +} + +fs::path unique_temp_path(const std::string &stem, const std::string &ext) { + std::error_code ec; + const fs::path base = fs::temp_directory_path(ec); + if (ec) { + throw std::runtime_error("Failed to get temp directory: " + ec.message()); + } + return base / (stem + "_" + std::to_string(std::rand()) + ext); +} + +void expect_legacy_to_trx_round_trip(const fs::path &input, const std::string &ref_nifti, const std::string &stem) { + trx::legacy::Tractogram tr; + if (input.extension() == ".tck") { + ASSERT_TRUE(trx::legacy::load_tck(input.string(), tr)); + } else if (input.extension() == ".trk") { + ASSERT_TRUE(trx::legacy::load_trk(input.string(), tr)); + } else if (input.extension() == ".vtk") { + ASSERT_TRUE(trx::legacy::load_vtk(input.string(), tr)); + } else { + FAIL() << "unhandled input extension: " << input.string(); + } + ASSERT_FALSE(tr.offsets.empty()); + + const size_t expected_vertices = tr.pts.size() / 3; + const size_t expected_streamlines = tr.offsets.size() - 1; + + const fs::path out = unique_temp_path(stem, ".trx"); + std::error_code ec; + fs::remove(out, ec); + ASSERT_TRUE(trx::legacy::save_trx(tr, out.string(), ref_nifti)); + + auto loaded = trx::load_any(out.string()); + EXPECT_EQ(loaded.num_vertices(), expected_vertices); + EXPECT_EQ(loaded.num_streamlines(), expected_streamlines); + loaded.close(); + + fs::remove(out, ec); +} + +} // namespace + +TEST(LegacyIo, TckToTrxRoundTripPreservesHeaderCounts) { + const std::string root_str = test_data_root(); + if (root_str.empty()) + GTEST_SKIP() << "Test data not found"; + const fs::path root(root_str); + const fs::path tck = root / "gs.tck"; + const fs::path nii = root / "gs.nii"; + if (!fs::exists(tck) || !fs::exists(nii)) { + GTEST_SKIP() << "gs.tck / gs.nii not present in test data"; + } + expect_legacy_to_trx_round_trip(tck, nii.string(), "trx_legacy_tck_roundtrip"); +} + +TEST(LegacyIo, TrkToTrxRoundTripPreservesHeaderCounts) { + const std::string root_str = test_data_root(); + if (root_str.empty()) + GTEST_SKIP() << "Test data not found"; + const fs::path root(root_str); + const fs::path trk = root / "gs.trk"; + if (!fs::exists(trk)) { + GTEST_SKIP() << "gs.trk not present in test data"; + } + expect_legacy_to_trx_round_trip(trk, "", "trx_legacy_trk_roundtrip"); +} + +TEST(LegacyIo, VtkToTrxRoundTripPreservesHeaderCounts) { + const std::string root_str = test_data_root(); + if (root_str.empty()) + GTEST_SKIP() << "Test data not found"; + const fs::path root(root_str); + const fs::path vtk = root / "gs.vtk"; + const fs::path nii = root / "gs.nii"; + if (!fs::exists(vtk) || !fs::exists(nii)) { + GTEST_SKIP() << "gs.vtk / gs.nii not present in test data"; + } + expect_legacy_to_trx_round_trip(vtk, nii.string(), "trx_legacy_vtk_roundtrip"); +} + +TEST(LegacyIo, LoadNiftiHeaderValid) { + const std::string root_str = test_data_root(); + if (root_str.empty()) + GTEST_SKIP() << "Test data not found"; + const fs::path nii = fs::path(root_str) / "gs.nii"; + if (!fs::exists(nii)) { + GTEST_SKIP() << "gs.nii not present in test data"; + } + json11::Json header; + ASSERT_TRUE(trx::legacy::load_nifti_header(nii.string(), header)); + EXPECT_TRUE(header["VOXEL_TO_RASMM"].is_array()); + EXPECT_EQ(header["VOXEL_TO_RASMM"].array_items().size(), 4u); + EXPECT_TRUE(header["DIMENSIONS"].is_array()); + EXPECT_EQ(header["DIMENSIONS"].array_items().size(), 3u); +} + +TEST(LegacyIo, VtkMalformedInputFailsGracefully) { + trx::legacy::Tractogram tr; + + // Non-existent file + EXPECT_FALSE(trx::legacy::load_vtk("/nonexistent/path/file.vtk", tr)); + + // Empty file + const fs::path empty_file = unique_temp_path("empty_vtk", ".vtk"); + { + std::ofstream out(empty_file); + } + EXPECT_FALSE(trx::legacy::load_vtk(empty_file.string(), tr)); + fs::remove(empty_file); + + // File claiming huge points count with no data + const fs::path huge_pts_file = unique_temp_path("huge_pts", ".vtk"); + { + std::ofstream out(huge_pts_file); + out << "# vtk DataFile Version 4.2\nvtk output\nBINARY\nDATASET POLYDATA\nPOINTS 18446744073709551600 float\n"; + } + EXPECT_FALSE(trx::legacy::load_vtk(huge_pts_file.string(), tr)); + fs::remove(huge_pts_file); + + // File claiming lines but truncated + const fs::path truncated_file = unique_temp_path("trunc_vtk", ".vtk"); + { + std::ofstream out(truncated_file, std::ios::binary); + out << "# vtk DataFile Version 4.2\nvtk output\nBINARY\nDATASET POLYDATA\nPOINTS 3 float\n"; + std::array pts = {1.0f, 2.0f, 3.0f}; + out.write(reinterpret_cast(pts.data()), sizeof(float) * pts.size()); + out << "LINES 10 100\n"; + } + EXPECT_FALSE(trx::legacy::load_vtk(truncated_file.string(), tr)); + fs::remove(truncated_file); + + // File with negative cell count + const fs::path neg_cell_file = unique_temp_path("neg_cell_vtk", ".vtk"); + { + std::ofstream out(neg_cell_file, std::ios::binary); + out << "# vtk DataFile Version 4.2\nvtk output\nBINARY\nDATASET POLYDATA\nPOINTS 3 float\n"; + std::array pts = {1.0f, 2.0f, 3.0f}; + out.write(reinterpret_cast(pts.data()), sizeof(float) * pts.size()); + out << "LINES 1 10\n"; + int32_t neg_count = -5; + out.write(reinterpret_cast(&neg_count), sizeof(neg_count)); + } + EXPECT_FALSE(trx::legacy::load_vtk(neg_cell_file.string(), tr)); + fs::remove(neg_cell_file); +} + +TEST(LegacyIo, InPlaceDirectorySavePreservesData) { + const std::string root_str = test_data_root(); + if (root_str.empty()) + GTEST_SKIP() << "Test data not found"; + const fs::path root(root_str); + const fs::path src_trx = root / "gs.trx"; + if (!fs::exists(src_trx)) { + GTEST_SKIP() << "gs.trx not present in test data"; + } + + const fs::path temp_dir = unique_temp_path("inplace_dir_save", "_dir"); + std::error_code ec; + fs::remove_all(temp_dir, ec); + + // Load archive and save as uncompressed directory + { + auto trx = trx::load_any(src_trx.string()); + trx::TrxSaveOptions opts; + opts.mode = trx::TrxSaveMode::Directory; + trx.save(temp_dir.string(), opts); + trx.close(); + } + + ASSERT_TRUE(fs::exists(temp_dir / "header.json")); + + // Now load from the directory and save IN-PLACE with overwrite=true + { + auto dir_trx = trx::load_any(temp_dir.string()); + const size_t orig_v = dir_trx.num_vertices(); + const size_t orig_s = dir_trx.num_streamlines(); + + trx::TrxSaveOptions opts; + opts.mode = trx::TrxSaveMode::Directory; + opts.overwrite_existing = true; + EXPECT_NO_THROW(dir_trx.save(temp_dir.string(), opts)); + dir_trx.close(); + + // Verify directory still exists and loads correctly + auto reloaded = trx::load_any(temp_dir.string()); + EXPECT_EQ(reloaded.num_vertices(), orig_v); + EXPECT_EQ(reloaded.num_streamlines(), orig_s); + reloaded.close(); + } + + fs::remove_all(temp_dir, ec); +} + +TEST(LegacyIo, OutOfRangeUint32GroupThrows) { + const std::string root_str = test_data_root(); + if (root_str.empty()) + GTEST_SKIP() << "Test data not found"; + const fs::path root(root_str); + const fs::path src_trx = root / "gs.trx"; + if (!fs::exists(src_trx)) { + GTEST_SKIP() << "gs.trx not present in test data"; + } + + const fs::path temp_dir = unique_temp_path("bad_group_test", "_dir"); + std::error_code ec; + fs::remove_all(temp_dir, ec); + + // Save as directory first + { + auto trx = trx::load_any(src_trx.string()); + trx::TrxSaveOptions opts; + opts.mode = trx::TrxSaveMode::Directory; + trx.save(temp_dir.string(), opts); + trx.close(); + } + + // Inject a group with an invalid index (index >= NB_STREAMLINES) + const fs::path groups_dir = temp_dir / "groups"; + fs::create_directories(groups_dir, ec); + const fs::path bad_group_file = groups_dir / "bad_group.uint32"; + { + std::ofstream out(bad_group_file, std::ios::binary); + uint32_t invalid_idx = 999999; + out.write(reinterpret_cast(&invalid_idx), sizeof(invalid_idx)); + } + + // Loading from directory should now detect the out-of-range uint32 group index and throw + EXPECT_THROW(trx::load_any(temp_dir.string()), trx::TrxFormatError); + + fs::remove_all(temp_dir, ec); +} 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); +}