From 7bcb214491b0c5d16d187e4239d05480acfc0482 Mon Sep 17 00:00:00 2001 From: Sameeul Samee Date: Sat, 21 Feb 2026 09:36:41 -0500 Subject: [PATCH 1/4] feat: add NIfTI-1/2 reader support (.nii and .nii.gz) Add NiftiReaderCPP C++ class and NiftiReader Python wrapper for reading NIFTI-1 and NIFTI-2 format files. Uses direct file I/O + zlib for decompression instead of TensorStore, following the same interface convention as TsReaderCPP. - Supports .nii (plain) and .nii.gz (gzip-compressed) files - Auto-detects NIFTI-1 (348-byte header) and NIFTI-2 (540-byte header) - Handles byte-swapped files - Applies scl_slope/scl_inter scaling when present (returns float64) - Exposes physical voxel size via pixdim[1..3] - 14 unit tests covering shapes, pixel values, subregion reads, scaling, compressed files, physical size, and error handling Co-Authored-By: Claude Sonnet 4.6 --- CMakeLists.txt | 13 +- README.md | 9 + ci-utils/install_prereq_linux.sh | 16 +- src/cpp/interface/interface.cpp | 36 ++ src/cpp/reader/niftireader.cpp | 559 ++++++++++++++++++++++++++++++ src/cpp/reader/niftireader.h | 75 ++++ src/cpp/utilities/utilities.h | 2 +- src/python/bfiocpp/__init__.py | 1 + src/python/bfiocpp/niftireader.py | 79 +++++ tests/test_nifti_reader.py | 331 ++++++++++++++++++ 10 files changed, 1108 insertions(+), 13 deletions(-) create mode 100644 src/cpp/reader/niftireader.cpp create mode 100644 src/cpp/reader/niftireader.h create mode 100644 src/python/bfiocpp/niftireader.py create mode 100644 tests/test_nifti_reader.py diff --git a/CMakeLists.txt b/CMakeLists.txt index 1b910f8..ab0cc13 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -22,13 +22,14 @@ if(NOT CMAKE_BUILD_TYPE) endif() set(CMAKE_CXX_FLAGS_RELEASE "-O2") -set(SOURCE +set(SOURCE src/cpp/ts_driver/tiled_tiff/omexml.cc src/cpp/ts_driver/tiled_tiff/tiled_tiff_key_value_store.cc src/cpp/ts_driver/ometiff/metadata.cc src/cpp/ts_driver/ometiff/driver.cc src/cpp/interface/interface.cpp src/cpp/reader/tsreader.cpp + src/cpp/reader/niftireader.cpp src/cpp/utilities/utilities.cpp src/cpp/writer/tswriter.cpp ) @@ -61,6 +62,9 @@ else () message(STATUS "Unable to find threads. bfio_cpp must have a threading library i.e. pthreads.") endif () +#==== ZLIB (system, required for NIfTI .nii.gz decompression) +find_package(ZLIB REQUIRED) + #==== Pybind11 find_package(pybind11 CONFIG REQUIRED) @@ -76,7 +80,8 @@ if(CMAKE_CXX_COMPILER_ID STREQUAL "GNU") target_link_libraries(libbfiocpp PRIVATE stdc++fs) endif() -target_link_libraries(libbfiocpp PRIVATE - tensorstore::tensorstore - tensorstore::all_drivers) +target_link_libraries(libbfiocpp PRIVATE + tensorstore::tensorstore + tensorstore::all_drivers + ZLIB::ZLIB) target_link_libraries(libbfiocpp PRIVATE ${Build_LIBRARIES}) \ No newline at end of file diff --git a/README.md b/README.md index ba2f88f..1015b69 100644 --- a/README.md +++ b/README.md @@ -1,5 +1,14 @@ This is the new backend for `bfio`, using `Tensorstore` and other high-throughput IO library +## Supported Formats + +| Format | Read | Write | +|--------|------|-------| +| OME-TIFF | Yes | No | +| OME-Zarr v2 | Yes | Yes | +| OME-Zarr v3 | Yes | Yes | +| NIfTI (.nii, .nii.gz) | Yes | No | + ## Build Requirements `bfiocpp` uses `Tensorstore` for reading and writing OME Tiff and OME Zarr files. So `Tensorstore` build requirements are needed to be satisfied for `bfiocpp` also. diff --git a/ci-utils/install_prereq_linux.sh b/ci-utils/install_prereq_linux.sh index 9587920..4850ff6 100755 --- a/ci-utils/install_prereq_linux.sh +++ b/ci-utils/install_prereq_linux.sh @@ -14,9 +14,9 @@ fi mkdir -p $LOCAL_INSTALL_DIR curl -L https://github.com/pybind/pybind11/archive/refs/tags/v2.12.0.zip -o v2.12.0.zip -unzip v2.12.0.zip +unzip -o v2.12.0.zip cd pybind11-2.12.0 -mkdir build_man +mkdir -p build_man cd build_man cmake -DCMAKE_INSTALL_PREFIX=../../$LOCAL_INSTALL_DIR/ -DPYBIND11_TEST=OFF .. make install -j4 @@ -24,9 +24,9 @@ cd ../../ curl -L https://github.com/madler/zlib/releases/download/v1.3.1/zlib131.zip -o zlib131.zip -unzip zlib131.zip +unzip -o zlib131.zip cd zlib-1.3.1 -mkdir build_man +mkdir -p build_man cd build_man cmake -DCMAKE_POLICY_VERSION_MINIMUM=3.5 -DCMAKE_INSTALL_PREFIX=/usr/local .. cmake --build . @@ -34,9 +34,9 @@ cmake --build . --target install cd ../../ curl -L https://github.com/libjpeg-turbo/libjpeg-turbo/archive/refs/tags/3.1.0.zip -o 3.1.0.zip -unzip 3.1.0.zip +unzip -o 3.1.0.zip cd libjpeg-turbo-3.1.0 -mkdir build_man +mkdir -p build_man cd build_man cmake -DCMAKE_INSTALL_PREFIX=/usr/local -DENABLE_STATIC=FALSE -DCMAKE_BUILD_TYPE=Release .. if [[ "$OSTYPE" == "darwin"* ]]; then @@ -49,9 +49,9 @@ make install -j4 cd ../../ curl -L https://github.com/glennrp/libpng/archive/refs/tags/v1.6.53.zip -o v1.6.53.zip -unzip v1.6.53.zip +unzip -o v1.6.53.zip cd libpng-1.6.53 -mkdir build_man +mkdir -p build_man cd build_man cmake -DCMAKE_POLICY_VERSION_MINIMUM=3.5 -DCMAKE_INSTALL_PREFIX=/usr/local .. make install -j4 diff --git a/src/cpp/interface/interface.cpp b/src/cpp/interface/interface.cpp index 96374ef..ace4d86 100644 --- a/src/cpp/interface/interface.cpp +++ b/src/cpp/interface/interface.cpp @@ -4,6 +4,7 @@ #include #include #include "../reader/tsreader.h" +#include "../reader/niftireader.h" #include "../utilities/sequence.h" #include "../utilities/utilities.h" #include "../writer/tswriter.h" @@ -69,6 +70,19 @@ py::array get_iterator_requested_tile_data(bfiocpp::TsReaderCPP& tl, std::int64 return as_pyarray_shared_5d(tmp, ih, iw, id, nc, nt) ; } +py::array get_nifti_image_data(bfiocpp::NiftiReaderCPP& nr, + const Seq& rows, const Seq& cols, + const Seq& layers, const Seq& channels, + const Seq& tsteps) { + auto tmp = nr.GetImageData(rows, cols, layers, channels, tsteps); + auto ih = rows.Stop() - rows.Start() + 1; + auto iw = cols.Stop() - cols.Start() + 1; + auto id = layers.Stop() - layers.Start() + 1; + auto nc = channels.Stop()- channels.Start()+ 1; + auto nt = tsteps.Stop() - tsteps.Start() + 1; + return as_pyarray_shared_5d(tmp, ih, iw, id, nc, nt); +} + PYBIND11_MODULE(libbfiocpp, m) { py::class_>(m, "Seq") .def(py::init()); @@ -116,7 +130,29 @@ PYBIND11_MODULE(libbfiocpp, m) { .value("OmeTiff", bfiocpp::FileType::OmeTiff) .value("OmeZarrV2", bfiocpp::FileType::OmeZarrV2) .value("OmeZarrV3", bfiocpp::FileType::OmeZarrV3) + .value("Nifti", bfiocpp::FileType::Nifti) + .value("NiftiGz", bfiocpp::FileType::NiftiGz) .export_values(); + + py::class_>(m, "NiftiReaderCPP") + .def(py::init()) + .def("get_image_height", &bfiocpp::NiftiReaderCPP::GetImageHeight) + .def("get_image_width", &bfiocpp::NiftiReaderCPP::GetImageWidth) + .def("get_image_depth", &bfiocpp::NiftiReaderCPP::GetImageDepth) + .def("get_tile_height", &bfiocpp::NiftiReaderCPP::GetTileHeight) + .def("get_tile_width", &bfiocpp::NiftiReaderCPP::GetTileWidth) + .def("get_tile_depth", &bfiocpp::NiftiReaderCPP::GetTileDepth) + .def("get_channel_count", &bfiocpp::NiftiReaderCPP::GetChannelCount) + .def("get_tstep_count", &bfiocpp::NiftiReaderCPP::GetTstepCount) + .def("get_datatype", &bfiocpp::NiftiReaderCPP::GetDataType) + .def("get_physical_size_x", &bfiocpp::NiftiReaderCPP::GetPhysicalSizeX) + .def("get_physical_size_y", &bfiocpp::NiftiReaderCPP::GetPhysicalSizeY) + .def("get_physical_size_z", &bfiocpp::NiftiReaderCPP::GetPhysicalSizeZ) + .def("get_image_data", + [](bfiocpp::NiftiReaderCPP& nr, const Seq& rows, const Seq& cols, + const Seq& layers, const Seq& channels, const Seq& tsteps) { + return get_nifti_image_data(nr, rows, cols, layers, channels, tsteps); + }, py::return_value_policy::reference); m.def("get_ome_xml", &bfiocpp::GetOmeXml); diff --git a/src/cpp/reader/niftireader.cpp b/src/cpp/reader/niftireader.cpp new file mode 100644 index 0000000..f0676fe --- /dev/null +++ b/src/cpp/reader/niftireader.cpp @@ -0,0 +1,559 @@ +#include "niftireader.h" + +#include +#include +#include +#include +#include +#include + +#include + +namespace bfiocpp { + +// --------------------------------------------------------------------------- +// NIFTI-1 header (348 bytes, packed) +// --------------------------------------------------------------------------- +#pragma pack(push, 1) +struct Nifti1Header { + int32_t sizeof_hdr; // 4 + char data_type[10]; // 14 + char db_name[18]; // 32 + int32_t extents; // 36 + int16_t session_error; // 38 + char regular; // 39 + char dim_info; // 40 + int16_t dim[8]; // 56 + float intent_p1; // 60 + float intent_p2; // 64 + float intent_p3; // 68 + int16_t intent_code; // 70 + int16_t datatype; // 72 + int16_t bitpix; // 74 + int16_t slice_start; // 76 + float pixdim[8]; // 108 + float vox_offset; // 112 + float scl_slope; // 116 + float scl_inter; // 120 + int16_t slice_end; // 122 + char slice_code; // 123 + char xyzt_units; // 124 + float cal_max; // 128 + float cal_min; // 132 + float slice_duration;// 136 + float toffset; // 140 + int32_t glmax; // 144 + int32_t glmin; // 148 + char descrip[80]; // 228 + char aux_file[24]; // 252 + int16_t qform_code; // 254 + int16_t sform_code; // 256 + float quatern_b; // 260 + float quatern_c; // 264 + float quatern_d; // 268 + float qoffset_x; // 272 + float qoffset_y; // 276 + float qoffset_z; // 280 + float srow_x[4]; // 296 + float srow_y[4]; // 312 + float srow_z[4]; // 328 + char intent_name[16];// 344 + char magic[4]; // 348 +}; +#pragma pack(pop) +static_assert(sizeof(Nifti1Header) == 348, "NIFTI-1 header must be 348 bytes"); + +// --------------------------------------------------------------------------- +// NIFTI-2 header (540 bytes, packed) +// --------------------------------------------------------------------------- +#pragma pack(push, 1) +struct Nifti2Header { + int32_t sizeof_hdr; // 4 + char magic[8]; // 12 + int16_t datatype; // 14 + int16_t bitpix; // 16 + int64_t dim[8]; // 80 + double intent_p1; // 88 + double intent_p2; // 96 + double intent_p3; // 104 + double pixdim[8]; // 168 + int64_t vox_offset; // 176 + double scl_slope; // 184 + double scl_inter; // 192 + double cal_max; // 200 + double cal_min; // 208 + double slice_duration;// 216 + double toffset; // 224 + int64_t slice_start; // 232 + int64_t slice_end; // 240 + char descrip[80]; // 320 + char aux_file[24]; // 344 + int32_t qform_code; // 348 + int32_t sform_code; // 352 + double quatern_b; // 360 + double quatern_c; // 368 + double quatern_d; // 376 + double qoffset_x; // 384 + double qoffset_y; // 392 + double qoffset_z; // 400 + double srow_x[4]; // 432 + double srow_y[4]; // 464 + double srow_z[4]; // 496 + int32_t slice_code; // 500 + int32_t xyzt_units; // 504 + int32_t intent_code; // 508 + char intent_name[16];// 524 + char dim_info; // 525 + char unused_str[15];// 540 +}; +#pragma pack(pop) +static_assert(sizeof(Nifti2Header) == 540, "NIFTI-2 header must be 540 bytes"); + +// --------------------------------------------------------------------------- +// Byte-swap helpers +// --------------------------------------------------------------------------- +static inline void bswap2(void* p) { + auto* b = reinterpret_cast(p); + std::swap(b[0], b[1]); +} +static inline void bswap4(void* p) { + auto* b = reinterpret_cast(p); + std::swap(b[0], b[3]); + std::swap(b[1], b[2]); +} +static inline void bswap8(void* p) { + auto* b = reinterpret_cast(p); + std::swap(b[0], b[7]); + std::swap(b[1], b[6]); + std::swap(b[2], b[5]); + std::swap(b[3], b[4]); +} + +// --------------------------------------------------------------------------- +// Map NIFTI datatype to bfiocpp data_type string +// --------------------------------------------------------------------------- +static std::string nifti_dtype_to_string(int16_t dt) { + switch (dt) { + case 2: return "uint8"; + case 4: return "int16"; + case 8: return "int32"; + case 16: return "float"; + case 64: return "double"; + case 256: return "int8"; + case 512: return "uint16"; + case 768: return "uint32"; + case 1024: return "int64"; + case 1280: return "uint64"; + default: return "unknown"; + } +} + +// --------------------------------------------------------------------------- +// Constructor +// --------------------------------------------------------------------------- +NiftiReaderCPP::NiftiReaderCPP(const std::string& fname) + : _fname(fname) +{ + // Detect gzip by extension + if (fname.size() >= 7 && + fname.substr(fname.size() - 7) == ".nii.gz") { + _is_gz = true; + } else if (fname.size() >= 4 && + fname.substr(fname.size() - 4) == ".nii") { + _is_gz = false; + } else { + throw std::runtime_error("NiftiReaderCPP: unsupported extension (expected .nii or .nii.gz)"); + } + ReadHeader(); +} + +// --------------------------------------------------------------------------- +// ReadHeader — handles NIFTI-1 and NIFTI-2, with byte-swap detection +// --------------------------------------------------------------------------- +void NiftiReaderCPP::ReadHeader() { + // For .nii.gz we must inflate first to read the header + if (_is_gz) { + EnsureDataLoaded(); + if (_raw_data.size() < 4) { + throw std::runtime_error("NiftiReaderCPP: compressed file too small"); + } + + int32_t sh; + std::memcpy(&sh, _raw_data.data(), sizeof(sh)); + + bool swapped = false; + int32_t sh_swapped = sh; + bswap4(&sh_swapped); + + if (sh == 348 || sh == 540) { + swapped = false; + } else if (sh_swapped == 348 || sh_swapped == 540) { + swapped = true; + sh = sh_swapped; + } else { + throw std::runtime_error("NiftiReaderCPP: invalid sizeof_hdr in compressed file"); + } + _byte_swapped = swapped; + + if (sh == 348) { + if (_raw_data.size() < 348) { + throw std::runtime_error("NiftiReaderCPP: compressed file too small for NIFTI-1 header"); + } + Nifti1Header hdr; + std::memcpy(&hdr, _raw_data.data(), sizeof(hdr)); + + if (_byte_swapped) { + for (int i = 0; i < 8; ++i) bswap2(&hdr.dim[i]); + bswap2(&hdr.datatype); + for (int i = 0; i < 8; ++i) bswap4(&hdr.pixdim[i]); + bswap4(&hdr.vox_offset); + bswap4(&hdr.scl_slope); + bswap4(&hdr.scl_inter); + } + + int ndim = hdr.dim[0]; + _nx = (ndim >= 1) ? hdr.dim[1] : 1; + _ny = (ndim >= 2) ? hdr.dim[2] : 1; + _nz = (ndim >= 3) ? hdr.dim[3] : 1; + _nt = (ndim >= 4 && hdr.dim[4] > 0) ? hdr.dim[4] : 1; + _datatype = hdr.datatype; + _vox_offset = hdr.vox_offset; + _scl_slope = hdr.scl_slope; + _scl_inter = hdr.scl_inter; + for (int i = 0; i < 8; ++i) _pixdim[i] = hdr.pixdim[i]; + } else { + // NIFTI-2 + if (_raw_data.size() < 540) { + throw std::runtime_error("NiftiReaderCPP: compressed file too small for NIFTI-2 header"); + } + Nifti2Header hdr; + std::memcpy(&hdr, _raw_data.data(), sizeof(hdr)); + + if (_byte_swapped) { + for (int i = 0; i < 8; ++i) bswap8(&hdr.dim[i]); + bswap2(&hdr.datatype); + for (int i = 0; i < 8; ++i) bswap8(&hdr.pixdim[i]); + bswap8(&hdr.vox_offset); + bswap8(&hdr.scl_slope); + bswap8(&hdr.scl_inter); + } + + int64_t ndim = hdr.dim[0]; + _nx = (ndim >= 1) ? hdr.dim[1] : 1; + _ny = (ndim >= 2) ? hdr.dim[2] : 1; + _nz = (ndim >= 3) ? hdr.dim[3] : 1; + _nt = (ndim >= 4 && hdr.dim[4] > 0) ? hdr.dim[4] : 1; + _datatype = hdr.datatype; + _vox_offset = static_cast(hdr.vox_offset); + _scl_slope = static_cast(hdr.scl_slope); + _scl_inter = static_cast(hdr.scl_inter); + for (int i = 0; i < 8; ++i) _pixdim[i] = static_cast(hdr.pixdim[i]); + } + } else { + // Plain .nii — read directly from file + std::ifstream f(_fname, std::ios::binary); + if (!f.is_open()) { + throw std::runtime_error("NiftiReaderCPP: cannot open file: " + _fname); + } + + int32_t sh = 0; + f.read(reinterpret_cast(&sh), sizeof(sh)); + if (!f) throw std::runtime_error("NiftiReaderCPP: cannot read sizeof_hdr"); + + bool swapped = false; + int32_t sh_swapped = sh; + bswap4(&sh_swapped); + + if (sh == 348 || sh == 540) { + swapped = false; + } else if (sh_swapped == 348 || sh_swapped == 540) { + swapped = true; + sh = sh_swapped; + } else { + throw std::runtime_error("NiftiReaderCPP: invalid sizeof_hdr"); + } + _byte_swapped = swapped; + + f.seekg(0, std::ios::beg); + + if (sh == 348) { + Nifti1Header hdr; + f.read(reinterpret_cast(&hdr), sizeof(hdr)); + if (!f) throw std::runtime_error("NiftiReaderCPP: cannot read NIFTI-1 header"); + + if (_byte_swapped) { + for (int i = 0; i < 8; ++i) bswap2(&hdr.dim[i]); + bswap2(&hdr.datatype); + for (int i = 0; i < 8; ++i) bswap4(&hdr.pixdim[i]); + bswap4(&hdr.vox_offset); + bswap4(&hdr.scl_slope); + bswap4(&hdr.scl_inter); + } + + int ndim = hdr.dim[0]; + _nx = (ndim >= 1) ? hdr.dim[1] : 1; + _ny = (ndim >= 2) ? hdr.dim[2] : 1; + _nz = (ndim >= 3) ? hdr.dim[3] : 1; + _nt = (ndim >= 4 && hdr.dim[4] > 0) ? hdr.dim[4] : 1; + _datatype = hdr.datatype; + _vox_offset = hdr.vox_offset; + _scl_slope = hdr.scl_slope; + _scl_inter = hdr.scl_inter; + for (int i = 0; i < 8; ++i) _pixdim[i] = hdr.pixdim[i]; + } else { + // NIFTI-2 + Nifti2Header hdr; + f.read(reinterpret_cast(&hdr), sizeof(hdr)); + if (!f) throw std::runtime_error("NiftiReaderCPP: cannot read NIFTI-2 header"); + + if (_byte_swapped) { + for (int i = 0; i < 8; ++i) bswap8(&hdr.dim[i]); + bswap2(&hdr.datatype); + for (int i = 0; i < 8; ++i) bswap8(&hdr.pixdim[i]); + bswap8(&hdr.vox_offset); + bswap8(&hdr.scl_slope); + bswap8(&hdr.scl_inter); + } + + int64_t ndim = hdr.dim[0]; + _nx = (ndim >= 1) ? hdr.dim[1] : 1; + _ny = (ndim >= 2) ? hdr.dim[2] : 1; + _nz = (ndim >= 3) ? hdr.dim[3] : 1; + _nt = (ndim >= 4 && hdr.dim[4] > 0) ? hdr.dim[4] : 1; + _datatype = hdr.datatype; + _vox_offset = static_cast(hdr.vox_offset); + _scl_slope = static_cast(hdr.scl_slope); + _scl_inter = static_cast(hdr.scl_inter); + for (int i = 0; i < 8; ++i) _pixdim[i] = static_cast(hdr.pixdim[i]); + } + } + + // Validate datatype + _data_type = nifti_dtype_to_string(_datatype); + if (_data_type == "unknown") { + throw std::runtime_error("NiftiReaderCPP: unsupported NIFTI datatype code: " + + std::to_string(_datatype)); + } + + // Map to bfiocpp convention: X=width, Y=height, Z=depth, T=tsteps, C=channels + _image_width = _nx; + _image_height = _ny; + _image_depth = (_nz > 0) ? _nz : 1; + _num_tsteps = (_nt > 0) ? _nt : 1; + _num_channels = 1; + + // If scl_slope is non-zero, output will be float64 + if (_scl_slope != 0.0f) { + _data_type = "double"; + } +} + +// --------------------------------------------------------------------------- +// EnsureDataLoaded — inflate .nii.gz (called once via std::once_flag) +// --------------------------------------------------------------------------- +void NiftiReaderCPP::EnsureDataLoaded() { + std::call_once(_load_flag, [this]() { + std::ifstream f(_fname, std::ios::binary | std::ios::ate); + if (!f.is_open()) { + throw std::runtime_error("NiftiReaderCPP: cannot open .nii.gz file: " + _fname); + } + auto compressed_size = static_cast(f.tellg()); + f.seekg(0, std::ios::beg); + + std::vector compressed(compressed_size); + f.read(reinterpret_cast(compressed.data()), compressed_size); + f.close(); + + // Inflate using zlib with gzip support (windowBits = 15+16) + z_stream stream{}; + stream.next_in = compressed.data(); + stream.avail_in = static_cast(compressed_size); + + if (inflateInit2(&stream, 16 + 15) != Z_OK) { + throw std::runtime_error("NiftiReaderCPP: inflateInit2 failed"); + } + + const std::size_t CHUNK = 1024 * 1024; // 1 MB chunks + _raw_data.clear(); + + int ret; + do { + std::size_t old_size = _raw_data.size(); + _raw_data.resize(old_size + CHUNK); + stream.next_out = _raw_data.data() + old_size; + stream.avail_out = static_cast(CHUNK); + ret = inflate(&stream, Z_NO_FLUSH); + if (ret == Z_STREAM_ERROR || ret == Z_DATA_ERROR || ret == Z_MEM_ERROR) { + inflateEnd(&stream); + throw std::runtime_error("NiftiReaderCPP: inflate error code=" + std::to_string(ret)); + } + // Trim the unused portion of the last chunk + _raw_data.resize(_raw_data.size() - stream.avail_out); + } while (ret == Z_OK); // Z_STREAM_END exits; Z_BUF_ERROR is benign + + inflateEnd(&stream); + }); +} + +// --------------------------------------------------------------------------- +// Byte-swap a single value of type T +// --------------------------------------------------------------------------- +template +static T maybe_swap(T val, bool swap) { + if (!swap) return val; + uint8_t buf[sizeof(T)]; + std::memcpy(buf, &val, sizeof(T)); + std::reverse(buf, buf + sizeof(T)); + std::memcpy(&val, buf, sizeof(T)); + return val; +} + +// --------------------------------------------------------------------------- +// ReadRegion — reads a rectangular sub-region without scaling +// --------------------------------------------------------------------------- +template +std::shared_ptr NiftiReaderCPP::ReadRegion( + const Seq& rows, const Seq& cols, const Seq& layers, + const Seq& channels, const Seq& tsteps) +{ + // Clamp to valid ranges (0-indexed, inclusive) + auto y0 = std::max(rows.Start(), 0); + auto y1 = std::min(rows.Stop(), _image_height - 1); + auto x0 = std::max(cols.Start(), 0); + auto x1 = std::min(cols.Stop(), _image_width - 1); + auto z0 = std::max(layers.Start(), 0); + auto z1 = std::min(layers.Stop(), _image_depth - 1); + auto t0 = std::max(tsteps.Start(), 0); + auto t1 = std::min(tsteps.Stop(), _num_tsteps - 1); + + std::size_t ny = static_cast(y1 - y0 + 1); + std::size_t nx = static_cast(x1 - x0 + 1); + std::size_t nz = static_cast(z1 - z0 + 1); + std::size_t nt = static_cast(t1 - t0 + 1); + + std::vector out; + out.reserve(nt * 1 * nz * ny * nx); + + // vox_offset is where pixel data starts in the file + std::size_t vox_start = static_cast(_vox_offset); + std::size_t elem_size = sizeof(T); + + if (_is_gz) { + EnsureDataLoaded(); + + for (long it = t0; it <= t1; ++it) { + for (long iz = z0; iz <= z1; ++iz) { + for (long iy = y0; iy <= y1; ++iy) { + for (long ix = x0; ix <= x1; ++ix) { + // NIFTI flat index: x-fastest + std::size_t flat = static_cast( + ix + _nx * (iy + _ny * (iz + _nz * it))); + std::size_t byte_off = vox_start + flat * elem_size; + if (byte_off + elem_size > _raw_data.size()) { + throw std::runtime_error("NiftiReaderCPP: out-of-bounds read in compressed data"); + } + T val; + std::memcpy(&val, _raw_data.data() + byte_off, elem_size); + out.push_back(maybe_swap(val, _byte_swapped)); + } + } + } + } + } else { + std::ifstream f(_fname, std::ios::binary); + if (!f.is_open()) { + throw std::runtime_error("NiftiReaderCPP: cannot open file: " + _fname); + } + + for (long it = t0; it <= t1; ++it) { + for (long iz = z0; iz <= z1; ++iz) { + for (long iy = y0; iy <= y1; ++iy) { + for (long ix = x0; ix <= x1; ++ix) { + std::size_t flat = static_cast( + ix + _nx * (iy + _ny * (iz + _nz * it))); + std::size_t byte_off = vox_start + flat * elem_size; + f.seekg(static_cast(byte_off), std::ios::beg); + T val; + f.read(reinterpret_cast(&val), elem_size); + if (!f) { + throw std::runtime_error("NiftiReaderCPP: read error at offset " + + std::to_string(byte_off)); + } + out.push_back(maybe_swap(val, _byte_swapped)); + } + } + } + } + } + + return std::make_shared(std::move(out)); +} + +// --------------------------------------------------------------------------- +// ReadRegionScaled — reads and applies scl_slope / scl_inter, returns double +// --------------------------------------------------------------------------- +std::shared_ptr NiftiReaderCPP::ReadRegionScaled( + const Seq& rows, const Seq& cols, const Seq& layers, + const Seq& channels, const Seq& tsteps) +{ + // Read raw data into a temporary image_data, then convert to double + scale + std::shared_ptr raw; + switch (_datatype) { + case 2: raw = ReadRegion (rows, cols, layers, channels, tsteps); break; + case 4: raw = ReadRegion (rows, cols, layers, channels, tsteps); break; + case 8: raw = ReadRegion (rows, cols, layers, channels, tsteps); break; + case 16: raw = ReadRegion (rows, cols, layers, channels, tsteps); break; + case 64: raw = ReadRegion (rows, cols, layers, channels, tsteps); break; + case 256: raw = ReadRegion (rows, cols, layers, channels, tsteps); break; + case 512: raw = ReadRegion(rows, cols, layers, channels, tsteps); break; + case 768: raw = ReadRegion(rows, cols, layers, channels, tsteps); break; + case 1024: raw = ReadRegion (rows, cols, layers, channels, tsteps); break; + case 1280: raw = ReadRegion(rows, cols, layers, channels, tsteps); break; + default: + throw std::runtime_error("NiftiReaderCPP: unsupported datatype for scaling"); + } + + double slope = static_cast(_scl_slope); + double inter = static_cast(_scl_inter); + + // Convert raw variant → std::vector with scaling applied + std::vector out; + std::visit([&](auto&& vec) { + out.reserve(vec.size()); + for (auto v : vec) { + out.push_back(static_cast(v) * slope + inter); + } + }, *raw); + + return std::make_shared(std::move(out)); +} + +// --------------------------------------------------------------------------- +// GetImageData — public dispatch +// --------------------------------------------------------------------------- +std::shared_ptr NiftiReaderCPP::GetImageData( + const Seq& rows, const Seq& cols, + const Seq& layers, const Seq& channels, const Seq& tsteps) +{ + if (_scl_slope != 0.0f) { + return ReadRegionScaled(rows, cols, layers, channels, tsteps); + } + + switch (_datatype) { + case 2: return ReadRegion (rows, cols, layers, channels, tsteps); + case 4: return ReadRegion (rows, cols, layers, channels, tsteps); + case 8: return ReadRegion (rows, cols, layers, channels, tsteps); + case 16: return ReadRegion (rows, cols, layers, channels, tsteps); + case 64: return ReadRegion (rows, cols, layers, channels, tsteps); + case 256: return ReadRegion (rows, cols, layers, channels, tsteps); + case 512: return ReadRegion(rows, cols, layers, channels, tsteps); + case 768: return ReadRegion(rows, cols, layers, channels, tsteps); + case 1024: return ReadRegion (rows, cols, layers, channels, tsteps); + case 1280: return ReadRegion(rows, cols, layers, channels, tsteps); + default: + throw std::runtime_error("NiftiReaderCPP: unsupported datatype: " + + std::to_string(_datatype)); + } +} + +} // namespace bfiocpp diff --git a/src/cpp/reader/niftireader.h b/src/cpp/reader/niftireader.h new file mode 100644 index 0000000..7310087 --- /dev/null +++ b/src/cpp/reader/niftireader.h @@ -0,0 +1,75 @@ +#pragma once + +#include +#include +#include +#include +#include +#include "../utilities/sequence.h" +#include "../reader/tsreader.h" // for image_data + +namespace bfiocpp { + +class NiftiReaderCPP { +public: + explicit NiftiReaderCPP(const std::string& fname); + + std::int64_t GetImageHeight() const { return _image_height; } + std::int64_t GetImageWidth() const { return _image_width; } + std::int64_t GetImageDepth() const { return _image_depth; } + std::int64_t GetChannelCount() const { return _num_channels; } + std::int64_t GetTstepCount() const { return _num_tsteps; } + // Tile dims = full image dims (no tiling for NIFTI) + std::int64_t GetTileHeight() const { return _image_height; } + std::int64_t GetTileWidth() const { return _image_width; } + std::int64_t GetTileDepth() const { return _image_depth; } + + std::string GetDataType() const { return _data_type; } + + float GetPhysicalSizeX() const { return _pixdim[1]; } + float GetPhysicalSizeY() const { return _pixdim[2]; } + float GetPhysicalSizeZ() const { return _pixdim[3]; } + + std::shared_ptr GetImageData( + const Seq& rows, const Seq& cols, + const Seq& layers, const Seq& channels, const Seq& tsteps); + +private: + void ReadHeader(); + void EnsureDataLoaded(); // inflate .nii.gz → _raw_data (thread-safe) + + template + std::shared_ptr ReadRegion( + const Seq& rows, const Seq& cols, const Seq& layers, + const Seq& channels, const Seq& tsteps); + + std::shared_ptr ReadRegionScaled( + const Seq& rows, const Seq& cols, const Seq& layers, + const Seq& channels, const Seq& tsteps); + + std::string _fname; + bool _is_gz = false; + bool _byte_swapped = false; + + // NIFTI dimension (stored as int64 to cover both NIFTI-1 and NIFTI-2) + int64_t _nx = 1, _ny = 1, _nz = 1, _nt = 1; + int16_t _datatype = 0; // raw NIFTI datatype code + double _vox_offset = 0.0; + float _scl_slope = 0.0f; + float _scl_inter = 0.0f; + float _pixdim[8] = {}; + + // bfiocpp-style dimensions + std::int64_t _image_width = 1; + std::int64_t _image_height = 1; + std::int64_t _image_depth = 1; + std::int64_t _num_channels = 1; + std::int64_t _num_tsteps = 1; + std::string _data_type; + + // Compressed data cache + std::vector _raw_data; + std::once_flag _load_flag; +}; + +} // namespace bfiocpp diff --git a/src/cpp/utilities/utilities.h b/src/cpp/utilities/utilities.h index 8195b34..5f71f02 100644 --- a/src/cpp/utilities/utilities.h +++ b/src/cpp/utilities/utilities.h @@ -11,7 +11,7 @@ namespace bfiocpp { -enum class FileType {OmeTiff, OmeZarrV2, OmeZarrV3}; +enum class FileType {OmeTiff, OmeZarrV2, OmeZarrV3, Nifti, NiftiGz}; tensorstore::Spec GetOmeTiffSpecToRead(const std::string& filename); tensorstore::Spec GetZarrSpecToRead(const std::string& filename, FileType ft); diff --git a/src/python/bfiocpp/__init__.py b/src/python/bfiocpp/__init__.py index 5e89052..1bb64ee 100644 --- a/src/python/bfiocpp/__init__.py +++ b/src/python/bfiocpp/__init__.py @@ -1,5 +1,6 @@ from .tsreader import TSReader, Seq, FileType, get_ome_xml # NOQA: F401 from .tswriter import TSWriter # NOQA: F401 +from .niftireader import NiftiReader # NOQA: F401 from . import _version __version__ = _version.get_versions()["version"] diff --git a/src/python/bfiocpp/niftireader.py b/src/python/bfiocpp/niftireader.py new file mode 100644 index 0000000..836f9ca --- /dev/null +++ b/src/python/bfiocpp/niftireader.py @@ -0,0 +1,79 @@ +import numpy as np +from .libbfiocpp import NiftiReaderCPP, Seq # NOQA: F401 + + +class NiftiReader: + """Read NIfTI-1/2 files (.nii and .nii.gz).""" + + def __init__(self, file_name: str) -> None: + self._reader = NiftiReaderCPP(file_name) + self._X: int = self._reader.get_image_width() + self._Y: int = self._reader.get_image_height() + self._Z: int = self._reader.get_image_depth() + self._C: int = self._reader.get_channel_count() + self._T: int = self._reader.get_tstep_count() + self._datatype: str = self._reader.get_datatype() + + @property + def X(self) -> int: + return self._X + + @property + def Y(self) -> int: + return self._Y + + @property + def Z(self) -> int: + return self._Z + + @property + def C(self) -> int: + return self._C + + @property + def T(self) -> int: + return self._T + + @property + def datatype(self) -> str: + return self._datatype + + @property + def physical_size_x(self) -> float: + return self._reader.get_physical_size_x() + + @property + def physical_size_y(self) -> float: + return self._reader.get_physical_size_y() + + @property + def physical_size_z(self) -> float: + return self._reader.get_physical_size_z() + + def data( + self, + rows: Seq, + cols: Seq, + layers: Seq = None, + channels: Seq = None, + tsteps: Seq = None, + ) -> np.ndarray: + if layers is None: + layers = Seq(0, 0, 1) + if channels is None: + channels = Seq(0, 0, 1) + if tsteps is None: + tsteps = Seq(0, 0, 1) + return self._reader.get_image_data(rows, cols, layers, channels, tsteps) + + def close(self) -> None: + pass + + def __enter__(self) -> "NiftiReader": + return self + + def __del__(self) -> None: + self.close() + + def __exit__(self, type_class, value, traceback) -> None: + self.close() diff --git a/tests/test_nifti_reader.py b/tests/test_nifti_reader.py new file mode 100644 index 0000000..c689c27 --- /dev/null +++ b/tests/test_nifti_reader.py @@ -0,0 +1,331 @@ +"""Tests for NIfTI-1 reader support in bfiocpp.""" + +import gzip +import os +import struct +import tempfile +import unittest + +import numpy as np + +from bfiocpp import NiftiReader +from bfiocpp.libbfiocpp import Seq + + +# --------------------------------------------------------------------------- +# Helpers to build minimal NIFTI-1 files +# --------------------------------------------------------------------------- + +def _make_nifti1_header( + nx, ny, nz, nt=1, + datatype=512, # uint16 by default + vox_offset=352.0, + scl_slope=0.0, + scl_inter=0.0, + pixdim=None, +): + """Return 348-byte NIFTI-1 header as bytes.""" + if pixdim is None: + pixdim = [1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0] + + # Determine ndim and dim array + if nt > 1: + ndim = 4 + elif nz > 1: + ndim = 3 + elif ny > 1: + ndim = 2 + else: + ndim = 1 + + dims = [ndim, nx, ny, nz, nt, 1, 1, 1] + + # Map datatype to bitpix + bitpix_map = {2: 8, 4: 16, 8: 32, 16: 32, 64: 64, + 256: 8, 512: 16, 768: 32, 1024: 64, 1280: 64} + bitpix = bitpix_map.get(datatype, 16) + + # Pack header fields (all little-endian) + # Use 'h' for regular+dim_info pair: low byte = regular='r', high byte = dim_info=0 + header = struct.pack( + ' bytes: + """Return raw bytes of a numpy array in C order (x-fastest in NIFTI convention).""" + # NIFTI stores x-fastest, which maps to column-major in the XY plane. + # For a 3-D array indexed [z, y, x] this is just C-contiguous bytes. + return array.tobytes(order='C') + + +def _make_nifti1_file(array: np.ndarray, datatype: int, tmp_dir: str, + fname="test.nii", scl_slope=0.0, scl_inter=0.0, + pixdim=None): + """Write a minimal NIFTI-1 .nii file and return the path.""" + shape = array.shape # (nz, ny, nx) or (nt, nz, ny, nx) + if array.ndim == 3: + nz, ny, nx = shape + nt = 1 + elif array.ndim == 4: + nt, nz, ny, nx = shape + else: + raise ValueError("Only 3-D and 4-D arrays supported") + + vox_offset = 352.0 # 348-byte header + 4-byte extension block + header = _make_nifti1_header( + nx=nx, ny=ny, nz=nz, nt=nt, + datatype=datatype, + vox_offset=vox_offset, + scl_slope=scl_slope, + scl_inter=scl_inter, + pixdim=pixdim, + ) + pixel_data = _nifti_pixel_data(array) + + path = os.path.join(tmp_dir, fname) + with open(path, 'wb') as f: + f.write(header) + f.write(b'\x00' * 4) # 4-byte extension block (all zeros = no extension) + f.write(pixel_data) + return path + + +def _make_nifti1_gz_file(array: np.ndarray, datatype: int, tmp_dir: str, + fname="test.nii.gz", scl_slope=0.0, scl_inter=0.0, + pixdim=None): + """Write a NIFTI-1 .nii.gz file and return the path.""" + nii_path = _make_nifti1_file(array, datatype, tmp_dir, + fname=fname.replace('.gz', ''), + scl_slope=scl_slope, scl_inter=scl_inter, + pixdim=pixdim) + gz_path = os.path.join(tmp_dir, fname) + with open(nii_path, 'rb') as f_in: + with gzip.open(gz_path, 'wb') as f_out: + f_out.write(f_in.read()) + os.remove(nii_path) + return gz_path + + +# --------------------------------------------------------------------------- +# Tests +# --------------------------------------------------------------------------- + +class TestNiftiRead(unittest.TestCase): + + @classmethod + def setUpClass(cls): + cls._tmpdir = tempfile.mkdtemp(prefix="bfiocpp_nifti_test_") + + # 3-D uint16 image: shape (nz=4, ny=5, nx=6) + cls.arr_u16 = np.arange(4 * 5 * 6, dtype=np.uint16).reshape(4, 5, 6) + cls.nii_u16 = _make_nifti1_file( + cls.arr_u16, datatype=512, tmp_dir=cls._tmpdir, fname="u16.nii" + ) + + # 3-D float32 image + cls.arr_f32 = np.random.rand(3, 4, 5).astype(np.float32) + cls.nii_f32 = _make_nifti1_file( + cls.arr_f32, datatype=16, tmp_dir=cls._tmpdir, fname="f32.nii" + ) + + # 3-D uint16 image with scaling: output must be float64 + cls.arr_scale_raw = np.ones((2, 3, 4), dtype=np.uint16) * 10 + cls.nii_scaled = _make_nifti1_file( + cls.arr_scale_raw, datatype=512, tmp_dir=cls._tmpdir, + fname="scaled.nii", scl_slope=2.0, scl_inter=1.0 + ) + + # .nii.gz version of the uint16 image + cls.arr_gz = np.arange(2 * 3 * 4, dtype=np.uint16).reshape(2, 3, 4) + cls.nii_gz = _make_nifti1_gz_file( + cls.arr_gz, datatype=512, tmp_dir=cls._tmpdir, fname="gz.nii.gz" + ) + + # 3-D image with physical size info + cls.arr_phys = np.zeros((2, 3, 4), dtype=np.float32) + cls.nii_phys = _make_nifti1_file( + cls.arr_phys, datatype=16, tmp_dir=cls._tmpdir, fname="phys.nii", + pixdim=[1.0, 0.5, 0.25, 2.0, 0.0, 0.0, 0.0, 0.0], + ) + + @classmethod + def tearDownClass(cls): + import shutil + shutil.rmtree(cls._tmpdir, ignore_errors=True) + + # ------------------------------------------------------------------ + # Shape / dimension tests + # ------------------------------------------------------------------ + + def test_shape_3d_uint16(self): + with NiftiReader(self.nii_u16) as r: + self.assertEqual(r.X, 6) # nx + self.assertEqual(r.Y, 5) # ny + self.assertEqual(r.Z, 4) # nz + self.assertEqual(r.C, 1) + self.assertEqual(r.T, 1) + + def test_shape_3d_float32(self): + with NiftiReader(self.nii_f32) as r: + self.assertEqual(r.X, 5) + self.assertEqual(r.Y, 4) + self.assertEqual(r.Z, 3) + + def test_tile_dims_equal_image_dims(self): + with NiftiReader(self.nii_u16) as r: + self.assertEqual(r._reader.get_tile_height(), r.Y) + self.assertEqual(r._reader.get_tile_width(), r.X) + self.assertEqual(r._reader.get_tile_depth(), r.Z) + + # ------------------------------------------------------------------ + # Pixel value tests + # ------------------------------------------------------------------ + + def test_pixel_values_uint16(self): + """Read entire volume and compare to original array.""" + with NiftiReader(self.nii_u16) as r: + arr = r.data( + rows=Seq(0, r.Y - 1, 1), + cols=Seq(0, r.X - 1, 1), + layers=Seq(0, r.Z - 1, 1), + ) + # arr shape: (T, C, Z, Y, X) + self.assertEqual(arr.shape, (1, 1, 4, 5, 6)) + self.assertEqual(arr.dtype, np.uint16) + # Compare against original (z, y, x) → (T, C, Z, Y, X) + np.testing.assert_array_equal(arr[0, 0], self.arr_u16) + + def test_pixel_values_float32(self): + with NiftiReader(self.nii_f32) as r: + arr = r.data( + rows=Seq(0, r.Y - 1, 1), + cols=Seq(0, r.X - 1, 1), + layers=Seq(0, r.Z - 1, 1), + ) + self.assertEqual(arr.dtype, np.float32) + np.testing.assert_array_almost_equal(arr[0, 0], self.arr_f32, decimal=6) + + def test_subregion_read(self): + """Read a sub-region and check values match the slice.""" + with NiftiReader(self.nii_u16) as r: + arr = r.data( + rows=Seq(1, 3, 1), # y 1..3 + cols=Seq(2, 4, 1), # x 2..4 + layers=Seq(0, 1, 1), # z 0..1 + ) + self.assertEqual(arr.shape, (1, 1, 2, 3, 3)) + expected = self.arr_u16[0:2, 1:4, 2:5] + np.testing.assert_array_equal(arr[0, 0], expected) + + # ------------------------------------------------------------------ + # Scaling test + # ------------------------------------------------------------------ + + def test_scaling(self): + """scl_slope=2, scl_inter=1 → output float64 with val*2+1.""" + with NiftiReader(self.nii_scaled) as r: + self.assertEqual(r.datatype, "double") + arr = r.data( + rows=Seq(0, r.Y - 1, 1), + cols=Seq(0, r.X - 1, 1), + layers=Seq(0, r.Z - 1, 1), + ) + self.assertEqual(arr.dtype, np.float64) + expected = self.arr_scale_raw.astype(np.float64) * 2.0 + 1.0 + np.testing.assert_array_almost_equal(arr[0, 0], expected, decimal=10) + + # ------------------------------------------------------------------ + # Compressed (.nii.gz) test + # ------------------------------------------------------------------ + + def test_gz_shape(self): + with NiftiReader(self.nii_gz) as r: + self.assertEqual(r.X, 4) + self.assertEqual(r.Y, 3) + self.assertEqual(r.Z, 2) + + def test_gz_pixel_values(self): + with NiftiReader(self.nii_gz) as r: + arr = r.data( + rows=Seq(0, r.Y - 1, 1), + cols=Seq(0, r.X - 1, 1), + layers=Seq(0, r.Z - 1, 1), + ) + self.assertEqual(arr.dtype, np.uint16) + np.testing.assert_array_equal(arr[0, 0], self.arr_gz) + + # ------------------------------------------------------------------ + # Physical size test + # ------------------------------------------------------------------ + + def test_physical_size(self): + with NiftiReader(self.nii_phys) as r: + self.assertAlmostEqual(r.physical_size_x, 0.5, places=5) + self.assertAlmostEqual(r.physical_size_y, 0.25, places=5) + self.assertAlmostEqual(r.physical_size_z, 2.0, places=5) + + # ------------------------------------------------------------------ + # Datatype string tests + # ------------------------------------------------------------------ + + def test_datatype_uint16(self): + with NiftiReader(self.nii_u16) as r: + self.assertEqual(r.datatype, "uint16") + + def test_datatype_float32(self): + with NiftiReader(self.nii_f32) as r: + self.assertEqual(r.datatype, "float") + + # ------------------------------------------------------------------ + # Error handling + # ------------------------------------------------------------------ + + def test_invalid_extension(self): + with self.assertRaises(RuntimeError): + NiftiReader("/tmp/nofile.tiff") + + def test_nonexistent_file(self): + with self.assertRaises(Exception): + NiftiReader("/tmp/does_not_exist.nii") + + +if __name__ == "__main__": + unittest.main() From 186c2d54f10a93f7104fd1f0df87fa2bc396ea4a Mon Sep 17 00:00:00 2001 From: Sameeul Samee Date: Fri, 6 Mar 2026 07:06:19 -0500 Subject: [PATCH 2/4] Add nifti writer support --- CMakeLists.txt | 1 + README.md | 2 +- ci-utils/install_prereq_linux.sh | 22 +- src/cpp/interface/interface.cpp | 15 +- src/cpp/writer/niftiwriter.cpp | 413 ++++++++++++++++++++++++++++++ src/cpp/writer/niftiwriter.h | 57 +++++ src/python/bfiocpp/__init__.py | 1 + src/python/bfiocpp/niftiwriter.py | 66 +++++ tests/test_nifti_writer.py | 230 +++++++++++++++++ 9 files changed, 797 insertions(+), 10 deletions(-) create mode 100644 src/cpp/writer/niftiwriter.cpp create mode 100644 src/cpp/writer/niftiwriter.h create mode 100644 src/python/bfiocpp/niftiwriter.py create mode 100644 tests/test_nifti_writer.py diff --git a/CMakeLists.txt b/CMakeLists.txt index ab0cc13..5c04c7b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -32,6 +32,7 @@ set(SOURCE src/cpp/reader/niftireader.cpp src/cpp/utilities/utilities.cpp src/cpp/writer/tswriter.cpp + src/cpp/writer/niftiwriter.cpp ) include(FetchContent) diff --git a/README.md b/README.md index 1015b69..c46fc0f 100644 --- a/README.md +++ b/README.md @@ -7,7 +7,7 @@ This is the new backend for `bfio`, using `Tensorstore` and other high-throughpu | OME-TIFF | Yes | No | | OME-Zarr v2 | Yes | Yes | | OME-Zarr v3 | Yes | Yes | -| NIfTI (.nii, .nii.gz) | Yes | No | +| NIfTI (.nii, .nii.gz) | Yes | Yes | ## Build Requirements diff --git a/ci-utils/install_prereq_linux.sh b/ci-utils/install_prereq_linux.sh index 4850ff6..4d38a0e 100755 --- a/ci-utils/install_prereq_linux.sh +++ b/ci-utils/install_prereq_linux.sh @@ -13,12 +13,22 @@ fi mkdir -p $LOCAL_INSTALL_DIR +# Set CMAKE_INSTALL_PREFIX to /usr/local when running in GitHub Actions, +# otherwise install into LOCAL_INSTALL_DIR. +if [ -n "$SYS_INSTALL" ]; then + DEPS_INSTALL_PREFIX="/usr/local" +else + DEPS_INSTALL_PREFIX=$LOCAL_INSTALL_DIR +fi + + + curl -L https://github.com/pybind/pybind11/archive/refs/tags/v2.12.0.zip -o v2.12.0.zip unzip -o v2.12.0.zip cd pybind11-2.12.0 mkdir -p build_man cd build_man -cmake -DCMAKE_INSTALL_PREFIX=../../$LOCAL_INSTALL_DIR/ -DPYBIND11_TEST=OFF .. +cmake -DCMAKE_INSTALL_PREFIX="$DEPS_INSTALL_PREFIX" -DPYBIND11_TEST=OFF .. make install -j4 cd ../../ @@ -28,8 +38,8 @@ unzip -o zlib131.zip cd zlib-1.3.1 mkdir -p build_man cd build_man -cmake -DCMAKE_POLICY_VERSION_MINIMUM=3.5 -DCMAKE_INSTALL_PREFIX=/usr/local .. -cmake --build . +cmake -DCMAKE_POLICY_VERSION_MINIMUM=3.5 -DCMAKE_INSTALL_PREFIX="$DEPS_INSTALL_PREFIX" .. +cmake --build . cmake --build . --target install cd ../../ @@ -38,21 +48,19 @@ unzip -o 3.1.0.zip cd libjpeg-turbo-3.1.0 mkdir -p build_man cd build_man -cmake -DCMAKE_INSTALL_PREFIX=/usr/local -DENABLE_STATIC=FALSE -DCMAKE_BUILD_TYPE=Release .. +cmake -DCMAKE_INSTALL_PREFIX="$DEPS_INSTALL_PREFIX" -DENABLE_STATIC=FALSE -DCMAKE_BUILD_TYPE=Release .. if [[ "$OSTYPE" == "darwin"* ]]; then sudo make install -j4 else make install -j4 fi cd ../../ -make install -j4 -cd ../../ curl -L https://github.com/glennrp/libpng/archive/refs/tags/v1.6.53.zip -o v1.6.53.zip unzip -o v1.6.53.zip cd libpng-1.6.53 mkdir -p build_man cd build_man -cmake -DCMAKE_POLICY_VERSION_MINIMUM=3.5 -DCMAKE_INSTALL_PREFIX=/usr/local .. +cmake -DCMAKE_POLICY_VERSION_MINIMUM=3.5 -DCMAKE_INSTALL_PREFIX="$DEPS_INSTALL_PREFIX" .. make install -j4 cd ../../ diff --git a/src/cpp/interface/interface.cpp b/src/cpp/interface/interface.cpp index ace4d86..13da81f 100644 --- a/src/cpp/interface/interface.cpp +++ b/src/cpp/interface/interface.cpp @@ -8,6 +8,7 @@ #include "../utilities/sequence.h" #include "../utilities/utilities.h" #include "../writer/tswriter.h" +#include "../writer/niftiwriter.h" namespace py = pybind11; using bfiocpp::Seq; @@ -156,8 +157,18 @@ PYBIND11_MODULE(libbfiocpp, m) { m.def("get_ome_xml", &bfiocpp::GetOmeXml); - - // Writer class + // NIfTI writer + py::class_>(m, "NiftiWriterCPP") + .def(py::init&, + const std::string&, const std::string&>(), + py::arg("filename"), + py::arg("image_shape"), + py::arg("dtype"), + py::arg("dimension_order")) + .def("write_image_data", &bfiocpp::NiftiWriterCPP::WriteImageData) + .def("close", &bfiocpp::NiftiWriterCPP::Close); + + // Zarr/OME-TIFF writer py::class_>(m, "TsWriterCPP") .def(py::init&, const std::vector&, const std::string&, const std::string&, bfiocpp::FileType>(), py::arg("filename"), diff --git a/src/cpp/writer/niftiwriter.cpp b/src/cpp/writer/niftiwriter.cpp new file mode 100644 index 0000000..13c88c6 --- /dev/null +++ b/src/cpp/writer/niftiwriter.cpp @@ -0,0 +1,413 @@ +#include "niftiwriter.h" + +#include +#include +#include +#include +#include +#include +#include + +#include + +namespace bfiocpp { + +// --------------------------------------------------------------------------- +// NIFTI-1 header (348 bytes, packed) — same layout as in niftireader.cpp +// --------------------------------------------------------------------------- +#pragma pack(push, 1) +struct Nifti1Header_W { + int32_t sizeof_hdr; + char data_type[10]; + char db_name[18]; + int32_t extents; + int16_t session_error; + char regular; + char dim_info; + int16_t dim[8]; + float intent_p1; + float intent_p2; + float intent_p3; + int16_t intent_code; + int16_t datatype; + int16_t bitpix; + int16_t slice_start; + float pixdim[8]; + float vox_offset; + float scl_slope; + float scl_inter; + int16_t slice_end; + char slice_code; + char xyzt_units; + float cal_max; + float cal_min; + float slice_duration; + float toffset; + int32_t glmax; + int32_t glmin; + char descrip[80]; + char aux_file[24]; + int16_t qform_code; + int16_t sform_code; + float quatern_b; + float quatern_c; + float quatern_d; + float qoffset_x; + float qoffset_y; + float qoffset_z; + float srow_x[4]; + float srow_y[4]; + float srow_z[4]; + char intent_name[16]; + char magic[4]; +}; +#pragma pack(pop) +static_assert(sizeof(Nifti1Header_W) == 348, "NIFTI-1 writer header must be 348 bytes"); + +// --------------------------------------------------------------------------- +// Map dtype string → (NIFTI datatype code, element size in bytes) +// --------------------------------------------------------------------------- +static std::pair dtype_string_to_nifti(const std::string& s) { + if (s == "uint8") return {2, 1}; + if (s == "int16") return {4, 2}; + if (s == "int32") return {8, 4}; + if (s == "float32" || s == "float") return {16, 4}; + if (s == "float64" || s == "double") return {64, 8}; + if (s == "int8") return {256, 1}; + if (s == "uint16") return {512, 2}; + if (s == "uint32") return {768, 4}; + if (s == "int64") return {1024, 8}; + if (s == "uint64") return {1280, 8}; + throw std::runtime_error("NiftiWriterCPP: unsupported dtype \"" + s + "\""); +} + +// --------------------------------------------------------------------------- +// Constructor +// --------------------------------------------------------------------------- +NiftiWriterCPP::NiftiWriterCPP( + const std::string& fname, + const std::vector& image_shape, + const std::string& dtype_str, + const std::string& dimension_order) + : _fname(fname), _is_gz(false), _closed(false), + _nx(1), _ny(1), _nz(1), _nt(1), + _vox_offset(352.0) +{ + // --- Detect file format --- + if (fname.size() >= 7 && fname.substr(fname.size() - 7) == ".nii.gz") { + _is_gz = true; + } else if (fname.size() >= 4 && fname.substr(fname.size() - 4) == ".nii") { + _is_gz = false; + } else { + throw std::runtime_error( + "NiftiWriterCPP: unsupported extension (expected .nii or .nii.gz): " + fname); + } + + // --- Validate and parse dimension_order (same logic as TsWriterCPP) --- + if (dimension_order.size() < 2 || dimension_order.size() > 5) { + throw std::invalid_argument( + "NiftiWriterCPP: dimension_order must contain 2–5 characters, got: \"" + + dimension_order + "\""); + } + for (char ch : dimension_order) { + if (ch != 'X' && ch != 'Y' && ch != 'Z' && ch != 'C' && ch != 'T') { + throw std::invalid_argument( + "NiftiWriterCPP: dimension_order contains invalid character, only " + "T/C/Z/Y/X allowed: \"" + dimension_order + "\""); + } + } + + auto pos = dimension_order.find('X'); + if (pos == std::string::npos) + throw std::invalid_argument("NiftiWriterCPP: dimension_order must contain 'X'"); + _x_index = static_cast(pos); + + pos = dimension_order.find('Y'); + if (pos == std::string::npos) + throw std::invalid_argument("NiftiWriterCPP: dimension_order must contain 'Y'"); + _y_index = static_cast(pos); + + pos = dimension_order.find('Z'); + if (pos != std::string::npos) _z_index.emplace(static_cast(pos)); + + pos = dimension_order.find('C'); + if (pos != std::string::npos) _c_index.emplace(static_cast(pos)); + + pos = dimension_order.find('T'); + if (pos != std::string::npos) _t_index.emplace(static_cast(pos)); + + // --- Extract image dimensions --- + if (static_cast(image_shape.size()) != static_cast(dimension_order.size())) { + throw std::invalid_argument( + "NiftiWriterCPP: image_shape length must match dimension_order length"); + } + _nx = image_shape[_x_index]; + _ny = image_shape[_y_index]; + if (_z_index.has_value()) _nz = image_shape[_z_index.value()]; + if (_t_index.has_value()) _nt = image_shape[_t_index.value()]; + + // NIFTI-1 dim fields are int16_t — validate + auto check_dim = [](int64_t d, const char* name) { + if (d <= 0 || d > std::numeric_limits::max()) + throw std::invalid_argument( + std::string("NiftiWriterCPP: dimension ") + name + + " = " + std::to_string(d) + " is out of NIFTI-1 range [1, 32767]"); + }; + check_dim(_nx, "X"); + check_dim(_ny, "Y"); + check_dim(_nz, "Z"); + check_dim(_nt, "T"); + + // --- Map dtype --- + auto [nifti_code, elem_size] = dtype_string_to_nifti(dtype_str); + _nifti_datatype = nifti_code; + _elem_size = elem_size; + + // --- Initialise storage --- + int64_t pixel_bytes = _nx * _ny * _nz * _nt * static_cast(_elem_size); + + if (_is_gz) { + // Buffer the entire image in memory; flush (compressed) on Close() + _pixel_buf.assign(static_cast(pixel_bytes), 0); + } else { + // Create the .nii file, write header, then zero-fill the pixel region + // so that partial WriteImageData calls work without gaps. + std::ofstream f(_fname, std::ios::binary | std::ios::trunc); + if (!f.is_open()) + throw std::runtime_error("NiftiWriterCPP: cannot create file: " + _fname); + + WriteNifti1Header(f); // writes 348-byte header + 4-byte extension block + + // Zero-fill pixel data region + const int64_t CHUNK = 65536; + std::vector zeros(static_cast(std::min(pixel_bytes, CHUNK)), 0); + int64_t remaining = pixel_bytes; + while (remaining > 0) { + int64_t n = std::min(remaining, CHUNK); + f.write(zeros.data(), static_cast(n)); + remaining -= n; + } + if (!f) + throw std::runtime_error("NiftiWriterCPP: error initialising file: " + _fname); + } +} + +// --------------------------------------------------------------------------- +// Destructor +// --------------------------------------------------------------------------- +NiftiWriterCPP::~NiftiWriterCPP() { + try { Close(); } catch (...) {} +} + +// --------------------------------------------------------------------------- +// WriteNifti1Header — write 348-byte header + 4-byte extension block +// --------------------------------------------------------------------------- +void NiftiWriterCPP::WriteNifti1Header(std::ostream& out) { + Nifti1Header_W hdr; + std::memset(&hdr, 0, sizeof(hdr)); + + hdr.sizeof_hdr = 348; + hdr.regular = 'r'; + + // ndim: minimum 2, go up if Z or T > 1 + int16_t ndim = 2; + if (_nz > 1) ndim = 3; + if (_nt > 1) ndim = 4; + + hdr.dim[0] = ndim; + hdr.dim[1] = static_cast(_nx); + hdr.dim[2] = static_cast(_ny); + hdr.dim[3] = static_cast(_nz); + hdr.dim[4] = static_cast(_nt); + hdr.dim[5] = 1; + hdr.dim[6] = 1; + hdr.dim[7] = 1; + + hdr.datatype = _nifti_datatype; + hdr.bitpix = static_cast(_elem_size * 8); + hdr.vox_offset = static_cast(_vox_offset); + hdr.scl_slope = 0.0f; // no scaling + + // Unit pixel spacing + hdr.pixdim[0] = 1.0f; + hdr.pixdim[1] = 1.0f; + hdr.pixdim[2] = 1.0f; + hdr.pixdim[3] = 1.0f; + hdr.pixdim[4] = 1.0f; + + // Magic for single-file NIfTI ("n+1\0") + hdr.magic[0] = 'n'; + hdr.magic[1] = '+'; + hdr.magic[2] = '1'; + hdr.magic[3] = '\0'; + + out.write(reinterpret_cast(&hdr), sizeof(hdr)); + + // 4-byte extension block (all zeros = no extensions) + const char ext[4] = {0, 0, 0, 0}; + out.write(ext, 4); +} + +// --------------------------------------------------------------------------- +// WriteImageData +// --------------------------------------------------------------------------- +void NiftiWriterCPP::WriteImageData( + const py::array& py_image, + const Seq& rows, + const Seq& cols, + const std::optional& layers, + const std::optional& channels, + const std::optional& tsteps) +{ + if (_closed) + throw std::runtime_error("NiftiWriterCPP: cannot write after Close()"); + + // Determine ranges + long y0 = rows.Start(), y1 = rows.Stop(); + long x0 = cols.Start(), x1 = cols.Stop(); + long z0 = 0, z1 = 0; + long t0 = 0, t1 = 0; + + if (_z_index.has_value() && layers.has_value()) { + z0 = layers->Start(); + z1 = layers->Stop(); + } + if (_t_index.has_value() && tsteps.has_value()) { + t0 = tsteps->Start(); + t1 = tsteps->Stop(); + } + + // Clamp to valid image bounds + x0 = std::max(x0, 0L); x1 = std::min(x1, _nx - 1); + y0 = std::max(y0, 0L); y1 = std::min(y1, _ny - 1); + z0 = std::max(z0, 0L); z1 = std::min(z1, _nz - 1); + t0 = std::max(t0, 0L); t1 = std::min(t1, _nt - 1); + + if (x0 > x1 || y0 > y1) + throw std::runtime_error("NiftiWriterCPP: empty write region"); + + // Get raw byte pointer from the numpy array + py::buffer_info info = py_image.request(); + const uint8_t* src = static_cast(info.ptr); + + long row_elems = x1 - x0 + 1; + long row_bytes = row_elems * _elem_size; + + if (_is_gz) { + // Write rows into the pixel buffer (memory-only) + long input_i = 0; + for (long it = t0; it <= t1; ++it) { + for (long iz = z0; iz <= z1; ++iz) { + for (long iy = y0; iy <= y1; ++iy) { + int64_t flat0 = x0 + _nx * (iy + _ny * (iz + _nz * it)); + int64_t buf_off = flat0 * _elem_size; + std::memcpy(_pixel_buf.data() + buf_off, + src + input_i * _elem_size, + static_cast(row_bytes)); + input_i += row_elems; + } + } + } + } else { + // Write rows directly into the .nii file at computed byte offsets + std::fstream f(_fname, std::ios::binary | std::ios::in | std::ios::out); + if (!f.is_open()) + throw std::runtime_error("NiftiWriterCPP: cannot open for writing: " + _fname); + + long input_i = 0; + for (long it = t0; it <= t1; ++it) { + for (long iz = z0; iz <= z1; ++iz) { + for (long iy = y0; iy <= y1; ++iy) { + int64_t flat0 = x0 + _nx * (iy + _ny * (iz + _nz * it)); + int64_t file_off = static_cast(_vox_offset) + flat0 * _elem_size; + f.seekp(static_cast(file_off), std::ios::beg); + f.write(reinterpret_cast(src + input_i * _elem_size), + static_cast(row_bytes)); + input_i += row_elems; + } + } + } + if (!f) + throw std::runtime_error("NiftiWriterCPP: write error in file: " + _fname); + } +} + +// --------------------------------------------------------------------------- +// Close — flush .nii.gz to disk (no-op for .nii; idempotent) +// --------------------------------------------------------------------------- +void NiftiWriterCPP::Close() { + if (_closed) return; + _closed = true; + if (_is_gz) { + FlushGzip(); + } +} + +// --------------------------------------------------------------------------- +// FlushGzip — compress header + _pixel_buf and write to the .nii.gz file +// --------------------------------------------------------------------------- +void NiftiWriterCPP::FlushGzip() { + // Build the 352-byte uncompressed header + extension block + std::ostringstream hdr_oss; + WriteNifti1Header(hdr_oss); + std::string hdr_bytes = hdr_oss.str(); + assert(hdr_bytes.size() == 352); + + // Open output file + std::ofstream fout(_fname, std::ios::binary | std::ios::trunc); + if (!fout) + throw std::runtime_error("NiftiWriterCPP: cannot create .nii.gz file: " + _fname); + + // Initialise gzip deflate stream + z_stream zs{}; + if (deflateInit2(&zs, Z_DEFAULT_COMPRESSION, Z_DEFLATED, + 16 + 15, 8, Z_DEFAULT_STRATEGY) != Z_OK) { + throw std::runtime_error("NiftiWriterCPP: deflateInit2 failed"); + } + + const uInt OUT_SZ = 65536; + std::vector out_buf(OUT_SZ); + + // Helper: feed a block of bytes through deflate, writing output to fout + auto deflate_block = [&](const uint8_t* data, std::size_t len, int flush_mode) { + const uint8_t* ptr = data; + std::size_t rem = len; + + // Feed in uInt-sized chunks (handles >4 GB gracefully) + do { + uInt chunk = static_cast( + std::min(rem, static_cast(std::numeric_limits::max()))); + zs.next_in = const_cast(ptr); + zs.avail_in = chunk; + ptr += chunk; + rem -= chunk; + + int cur_flush = (rem == 0) ? flush_mode : Z_NO_FLUSH; + int ret; + do { + zs.next_out = out_buf.data(); + zs.avail_out = OUT_SZ; + ret = deflate(&zs, cur_flush); + if (ret == Z_STREAM_ERROR) { + deflateEnd(&zs); + throw std::runtime_error("NiftiWriterCPP: deflate stream error"); + } + uInt have = OUT_SZ - zs.avail_out; + fout.write(reinterpret_cast(out_buf.data()), + static_cast(have)); + } while (zs.avail_out == 0); + } while (rem > 0); + }; + + // Compress header, then pixel data + deflate_block(reinterpret_cast(hdr_bytes.data()), + hdr_bytes.size(), Z_NO_FLUSH); + deflate_block(_pixel_buf.data(), _pixel_buf.size(), Z_FINISH); + + deflateEnd(&zs); + + if (!fout) + throw std::runtime_error("NiftiWriterCPP: write error creating .nii.gz: " + _fname); +} + +} // namespace bfiocpp diff --git a/src/cpp/writer/niftiwriter.h b/src/cpp/writer/niftiwriter.h new file mode 100644 index 0000000..524cbc0 --- /dev/null +++ b/src/cpp/writer/niftiwriter.h @@ -0,0 +1,57 @@ +#pragma once + +#include +#include +#include +#include +#include +#include "../utilities/sequence.h" + +namespace py = pybind11; + +namespace bfiocpp { + +class NiftiWriterCPP { +public: + NiftiWriterCPP( + const std::string& fname, + const std::vector& image_shape, + const std::string& dtype_str, + const std::string& dimension_order + ); + ~NiftiWriterCPP(); + + void WriteImageData( + const py::array& py_image, + const Seq& rows, + const Seq& cols, + const std::optional& layers, + const std::optional& channels, + const std::optional& tsteps + ); + + void Close(); // explicit flush for .nii.gz; safe to call multiple times + +private: + void WriteNifti1Header(std::ostream& out); + void FlushGzip(); + + std::string _fname; + bool _is_gz; + bool _closed; + + // Image dimensions (NIFTI order: X=width, Y=height, Z=depth, T=tsteps) + int64_t _nx, _ny, _nz, _nt; + int16_t _nifti_datatype; + int _elem_size; // bytes per voxel + double _vox_offset; // 352.0 for NIFTI-1 single-file + + // Index of each dimension in the caller's image_shape / Seq arguments + int _x_index, _y_index; + std::optional _z_index, _c_index, _t_index; + + // Pixel buffer for .nii.gz (full image; flushed on Close()) + std::vector _pixel_buf; +}; + +} // namespace bfiocpp diff --git a/src/python/bfiocpp/__init__.py b/src/python/bfiocpp/__init__.py index 1bb64ee..0983332 100644 --- a/src/python/bfiocpp/__init__.py +++ b/src/python/bfiocpp/__init__.py @@ -1,6 +1,7 @@ from .tsreader import TSReader, Seq, FileType, get_ome_xml # NOQA: F401 from .tswriter import TSWriter # NOQA: F401 from .niftireader import NiftiReader # NOQA: F401 +from .niftiwriter import NiftiWriter # NOQA: F401 from . import _version __version__ = _version.get_versions()["version"] diff --git a/src/python/bfiocpp/niftiwriter.py b/src/python/bfiocpp/niftiwriter.py new file mode 100644 index 0000000..604b72b --- /dev/null +++ b/src/python/bfiocpp/niftiwriter.py @@ -0,0 +1,66 @@ +import numpy as np +from .libbfiocpp import NiftiWriterCPP, Seq # NOQA: F401 + + +class NiftiWriter: + """Write NIfTI-1 files (.nii and .nii.gz). + + Parameters + ---------- + file_name : str + Output path ending in ``.nii`` or ``.nii.gz``. + image_shape : list[int] + Full image shape matching *dimension_order*. + dtype : np.dtype + Element data type (e.g. ``np.dtype("uint16")``). + dimension_order : str + A combination of ``T``, ``C``, ``Z``, ``Y``, ``X`` + describing the axis order of *image_shape*. + Must contain at least ``X`` and ``Y``. + ``C`` is accepted but ignored by NIfTI. + """ + + def __init__( + self, + file_name: str, + image_shape: list, + dtype: np.dtype, + dimension_order: str, + ) -> None: + self._writer = NiftiWriterCPP( + file_name, + [int(d) for d in image_shape], + str(np.dtype(dtype)), + dimension_order, + ) + self._dtype = np.dtype(dtype) + + def write_image_data( + self, + image_data: np.ndarray, + rows: Seq, + cols: Seq, + layers: Seq = None, + channels: Seq = None, + tsteps: Seq = None, + ) -> None: + if not isinstance(image_data, np.ndarray): + raise ValueError("image_data must be a numpy ndarray") + # Ensure correct dtype and contiguous C layout before passing to C++ + flat = np.ascontiguousarray(image_data, dtype=self._dtype).flatten() + self._writer.write_image_data(flat, rows, cols, layers, channels, tsteps) + + def close(self) -> None: + self._writer.close() + + def __enter__(self) -> "NiftiWriter": + return self + + def __exit__(self, exc_type, exc_val, exc_tb) -> None: + self.close() + + def __del__(self) -> None: + try: + self.close() + except Exception: + pass diff --git a/tests/test_nifti_writer.py b/tests/test_nifti_writer.py new file mode 100644 index 0000000..583de65 --- /dev/null +++ b/tests/test_nifti_writer.py @@ -0,0 +1,230 @@ +"""Round-trip tests for NIfTI-1 writer support in bfiocpp.""" + +import os +import shutil +import tempfile +import unittest + +import numpy as np + +from bfiocpp import NiftiReader, NiftiWriter +from bfiocpp.libbfiocpp import Seq + + +class TestNiftiWrite(unittest.TestCase): + + @classmethod + def setUpClass(cls): + cls._tmpdir = tempfile.mkdtemp(prefix="bfiocpp_nifti_write_test_") + + @classmethod + def tearDownClass(cls): + shutil.rmtree(cls._tmpdir, ignore_errors=True) + + def _out(self, name): + return os.path.join(self._tmpdir, name) + + # ------------------------------------------------------------------ + # Helper: read the full volume back and return it as (Z, Y, X) or + # (T, Z, Y, X) shaped array. + # ------------------------------------------------------------------ + def _readback(self, path, shape_zyx): + """Read entire volume; returns array shaped (Z, Y, X).""" + with NiftiReader(path) as r: + nx, ny, nz = r.X, r.Y, r.Z + arr = r.data( + rows=Seq(0, ny - 1, 1), + cols=Seq(0, nx - 1, 1), + layers=Seq(0, nz - 1, 1), + ) + # arr shape: (T=1, C=1, Z, Y, X) + return arr[0, 0] + + # ------------------------------------------------------------------ + # test_write_read_uint16_3d + # ------------------------------------------------------------------ + def test_write_read_uint16_3d(self): + """Write a ZYX uint16 volume and read it back.""" + arr = np.arange(4 * 5 * 6, dtype=np.uint16).reshape(4, 5, 6) + path = self._out("u16_3d.nii") + + with NiftiWriter(path, [4, 5, 6], np.dtype("uint16"), "ZYX") as w: + w.write_image_data( + arr, + rows=Seq(0, 4, 1), + cols=Seq(0, 5, 1), + layers=Seq(0, 3, 1), + ) + + result = self._readback(path, arr.shape) + self.assertEqual(result.dtype, np.uint16) + np.testing.assert_array_equal(result, arr) + + # ------------------------------------------------------------------ + # test_write_read_float32_3d + # ------------------------------------------------------------------ + def test_write_read_float32_3d(self): + """Write a ZYX float32 volume and read it back.""" + rng = np.random.default_rng(42) + arr = rng.random((3, 4, 5)).astype(np.float32) + path = self._out("f32_3d.nii") + + with NiftiWriter(path, [3, 4, 5], np.dtype("float32"), "ZYX") as w: + w.write_image_data( + arr, + rows=Seq(0, 3, 1), + cols=Seq(0, 4, 1), + layers=Seq(0, 2, 1), + ) + + result = self._readback(path, arr.shape) + self.assertEqual(result.dtype, np.float32) + np.testing.assert_array_almost_equal(result, arr, decimal=6) + + # ------------------------------------------------------------------ + # test_write_read_gz + # ------------------------------------------------------------------ + def test_write_read_gz(self): + """Write a .nii.gz file and read it back.""" + arr = np.arange(2 * 3 * 4, dtype=np.uint16).reshape(2, 3, 4) + path = self._out("gz_test.nii.gz") + + with NiftiWriter(path, [2, 3, 4], np.dtype("uint16"), "ZYX") as w: + w.write_image_data( + arr, + rows=Seq(0, 2, 1), + cols=Seq(0, 3, 1), + layers=Seq(0, 1, 1), + ) + + result = self._readback(path, arr.shape) + self.assertEqual(result.dtype, np.uint16) + np.testing.assert_array_equal(result, arr) + + # ------------------------------------------------------------------ + # test_write_subregion + # ------------------------------------------------------------------ + def test_write_subregion(self): + """Two separate partial WriteImageData calls combine correctly.""" + # Full image: Z=4, Y=5, X=6, uint16 + full = np.arange(4 * 5 * 6, dtype=np.uint16).reshape(4, 5, 6) + path = self._out("subregion.nii") + + with NiftiWriter(path, [4, 5, 6], np.dtype("uint16"), "ZYX") as w: + # Write first two z-slices + w.write_image_data( + full[:2], + rows=Seq(0, 4, 1), + cols=Seq(0, 5, 1), + layers=Seq(0, 1, 1), + ) + # Write last two z-slices + w.write_image_data( + full[2:], + rows=Seq(0, 4, 1), + cols=Seq(0, 5, 1), + layers=Seq(2, 3, 1), + ) + + result = self._readback(path, full.shape) + np.testing.assert_array_equal(result, full) + + # ------------------------------------------------------------------ + # test_write_read_4d + # ------------------------------------------------------------------ + def test_write_read_4d(self): + """Write two time-steps of a TZYX volume and read back.""" + # shape: T=2, Z=3, Y=4, X=5 + arr = np.arange(2 * 3 * 4 * 5, dtype=np.int32).reshape(2, 3, 4, 5) + path = self._out("4d_tzyx.nii") + + with NiftiWriter(path, [2, 3, 4, 5], np.dtype("int32"), "TZYX") as w: + for t in range(2): + w.write_image_data( + arr[t], + rows=Seq(0, 3, 1), + cols=Seq(0, 4, 1), + layers=Seq(0, 2, 1), + tsteps=Seq(t, t, 1), + ) + + # Read back both time steps + with NiftiReader(path) as r: + nx, ny, nz, nt = r.X, r.Y, r.Z, r.T + self.assertEqual(nt, 2) + for t in range(2): + data = r.data( + rows=Seq(0, ny - 1, 1), + cols=Seq(0, nx - 1, 1), + layers=Seq(0, nz - 1, 1), + tsteps=Seq(t, t, 1), + ) + # data shape: (T=1, C=1, Z, Y, X) + np.testing.assert_array_equal(data[0, 0], arr[t]) + + # ------------------------------------------------------------------ + # test_write_subregion_xy + # ------------------------------------------------------------------ + def test_write_subregion_xy(self): + """Partial writes in X/Y produce the correct result.""" + full = np.zeros((4, 6), dtype=np.uint16) + full[1:3, 2:5] = np.array([[10, 11, 12], [20, 21, 22]], dtype=np.uint16) + path = self._out("subregion_xy.nii") + + with NiftiWriter(path, [4, 6], np.dtype("uint16"), "YX") as w: + # Write just the interior patch + patch = full[1:3, 2:5] + w.write_image_data( + patch, + rows=Seq(1, 2, 1), + cols=Seq(2, 4, 1), + ) + + result = self._readback(path, (1, 4, 6)) + # Z=1 slice; compare (Y, X) grid + np.testing.assert_array_equal(result[0], full) + + # ------------------------------------------------------------------ + # test_invalid_dtype + # ------------------------------------------------------------------ + def test_invalid_dtype(self): + """Unknown dtype string raises RuntimeError.""" + with self.assertRaises(Exception): + NiftiWriter(self._out("bad_dtype.nii"), [4, 5], np.dtype("complex64"), "YX") + + # ------------------------------------------------------------------ + # test_invalid_dimension_order_missing_x + # ------------------------------------------------------------------ + def test_invalid_dimension_order_missing_x(self): + """dimension_order without 'X' raises an error.""" + with self.assertRaises(Exception): + NiftiWriter(self._out("bad_order.nii"), [4, 5], np.dtype("uint16"), "YZ") + + # ------------------------------------------------------------------ + # test_invalid_dimension_order_missing_y + # ------------------------------------------------------------------ + def test_invalid_dimension_order_missing_y(self): + """dimension_order without 'Y' raises an error.""" + with self.assertRaises(Exception): + NiftiWriter(self._out("bad_order2.nii"), [4, 5], np.dtype("uint16"), "XZ") + + # ------------------------------------------------------------------ + # test_context_manager + # ------------------------------------------------------------------ + def test_context_manager_gz(self): + """Context manager correctly flushes a .nii.gz file.""" + arr = np.ones((2, 3, 4), dtype=np.float32) * 3.14 + path = self._out("ctx_mgr.nii.gz") + + with NiftiWriter(path, [2, 3, 4], np.dtype("float32"), "ZYX") as w: + w.write_image_data(arr, Seq(0, 2, 1), Seq(0, 3, 1), Seq(0, 1, 1)) + + self.assertTrue(os.path.exists(path)) + self.assertGreater(os.path.getsize(path), 0) + + result = self._readback(path, arr.shape) + np.testing.assert_array_almost_equal(result, arr, decimal=5) + + +if __name__ == "__main__": + unittest.main() From e36a90b7f72642ec84b75aa4a643cd3daf42fed5 Mon Sep 17 00:00:00 2001 From: Sameeul Samee Date: Fri, 6 Mar 2026 11:58:09 -0500 Subject: [PATCH 3/4] update path --- ci-utils/install_prereq_linux.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ci-utils/install_prereq_linux.sh b/ci-utils/install_prereq_linux.sh index 4d38a0e..90d291d 100755 --- a/ci-utils/install_prereq_linux.sh +++ b/ci-utils/install_prereq_linux.sh @@ -18,7 +18,7 @@ mkdir -p $LOCAL_INSTALL_DIR if [ -n "$SYS_INSTALL" ]; then DEPS_INSTALL_PREFIX="/usr/local" else - DEPS_INSTALL_PREFIX=$LOCAL_INSTALL_DIR + DEPS_INSTALL_PREFIX="$PWD/$LOCAL_INSTALL_DIR" fi From 2fa606471b0e38569fa31420a27016a0c5ca3ad5 Mon Sep 17 00:00:00 2001 From: Sameeul Samee Date: Fri, 6 Mar 2026 13:21:42 -0500 Subject: [PATCH 4/4] Fix std:min issue --- src/cpp/writer/niftiwriter.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/cpp/writer/niftiwriter.cpp b/src/cpp/writer/niftiwriter.cpp index 13c88c6..9994745 100644 --- a/src/cpp/writer/niftiwriter.cpp +++ b/src/cpp/writer/niftiwriter.cpp @@ -278,10 +278,10 @@ void NiftiWriterCPP::WriteImageData( } // Clamp to valid image bounds - x0 = std::max(x0, 0L); x1 = std::min(x1, _nx - 1); - y0 = std::max(y0, 0L); y1 = std::min(y1, _ny - 1); - z0 = std::max(z0, 0L); z1 = std::min(z1, _nz - 1); - t0 = std::max(t0, 0L); t1 = std::min(t1, _nt - 1); + x0 = std::max(x0, 0L); x1 = std::min(x1, static_cast(_nx - 1)); + y0 = std::max(y0, 0L); y1 = std::min(y1, static_cast(_ny - 1)); + z0 = std::max(z0, 0L); z1 = std::min(z1, static_cast(_nz - 1)); + t0 = std::max(t0, 0L); t1 = std::min(t1, static_cast(_nt - 1)); if (x0 > x1 || y0 > y1) throw std::runtime_error("NiftiWriterCPP: empty write region");