From 4177d722afb4d5d398708051c80a0c7142b71618 Mon Sep 17 00:00:00 2001 From: JesseMckinzie Date: Sun, 10 Nov 2024 21:12:02 -0700 Subject: [PATCH 01/24] Update PyramidCompositorCPP --- CMakeLists.txt | 2 + src/cpp/core/pyramid_compositor.cpp | 514 ++++++++++++++++++ src/cpp/core/pyramid_compositor.h | 100 ++++ .../interface/pyramid_python_interface.cpp | 12 + src/cpp/utilities/utilities.cpp | 94 ++++ src/cpp/utilities/utilities.h | 12 + src/python/argolid/pyramid_compositor.py | 100 +++- 7 files changed, 833 insertions(+), 1 deletion(-) create mode 100644 src/cpp/core/pyramid_compositor.cpp create mode 100644 src/cpp/core/pyramid_compositor.h diff --git a/CMakeLists.txt b/CMakeLists.txt index e8ee5c3..3376210 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -31,6 +31,8 @@ set(SOURCE src/cpp/core/chunked_pyramid_assembler.cpp src/cpp/core/chunked_base_to_pyr_gen.cpp src/cpp/core/ome_tiff_to_chunked_pyramid.cpp + src/cpp/core/pyramid_compositor.cpp + src/cpp/core/pyramid_view.cpp src/cpp/core/pyramid_view.cpp src/cpp/utilities/utilities.cpp ) diff --git a/src/cpp/core/pyramid_compositor.cpp b/src/cpp/core/pyramid_compositor.cpp new file mode 100644 index 0000000..7732fd1 --- /dev/null +++ b/src/cpp/core/pyramid_compositor.cpp @@ -0,0 +1,514 @@ +#include "pyramid_compositor.h" +#include "../utilities/utilities.h" + +#include +#include +#include +#include + +#include "tensorstore/context.h" +#include "tensorstore/array.h" +#include "tensorstore/driver/zarr/dtype.h" +#include "tensorstore/index_space/dim_expression.h" +#include "tensorstore/kvstore/kvstore.h" +#include "tensorstore/open.h" + +#include + +using ::tensorstore::internal_zarr::ChooseBaseDType; +using json = nlohmann::json; + +namespace fs = std::filesystem; + +namespace argolid { +PyramidCompositor::PyramidCompositor(const std::string& input_pyramids_loc, const std::string& out_dir, const std::string& output_pyramid_name): + _input_pyramids_loc(input_pyramids_loc), + _output_pyramid_name(out_dir + "/" + output_pyramid_name), + _ome_metadata_file(out_dir + "/" + output_pyramid_name + "/METADATA.ome.xml"), + _pyramid_levels(-1), + _num_channels(-1) { + + + tensorstore::TensorStore source; + + TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( + GetZarrSpecToRead(input_pyramids_loc), + tensorstore::OpenMode::open, + tensorstore::ReadWriteMode::read).result()); + + _image_dtype = source.dtype().name(); + _data_type_code = GetDataTypeCode(_image_dtype); + + // auto read_spec = GetZarrSpecToRead(input_pyramids_loc) + + // TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( + // read_spec, + // tensorstore::OpenMode::open, + // tensorstore::ReadWriteMode::read).result()); + + // auto image_shape = source.domain().shape(); + // const auto read_chunk_shape = source.chunk_layout().value().read_chunk_shape(); + + // if (image_shape.size() == 5){ + // _image_height = image_shape[3]; + // _image_width = image_shape[4]; + // _image_depth = image_shape[2]; + // _num_channels = image_shape[1]; + // _num_tsteps = image_shape[0]; + // _t_int.emplace(0); + // _c_int.emplace(1); + // _z_int.emplace(2); + // } else { + // assert(image_shape.size() >= 2); + // std::tie(_t_int, _c_int, _z_int) = ParseMultiscaleMetadata(axes_list, image_shape.size()); + // _image_height = image_shape[image_shape.size()-2]; + // _image_width = image_shape[image_shape.size()-1]; + // if (_t_int.has_value()) { + // _num_tsteps = image_shape[_t_int.value()]; + // } else { + // _num_tsteps = 1; + // } + + // if (_c_int.has_value()) { + // _num_channels = image_shape[_c_int.value()]; + // } else { + // _num_channels = 1; + // } + + // if (_z_int.has_value()) { + // _image_depth = image_shape[_z_int.value()]; + // } else { + // _image_depth = 1; + // } + // } + + // _data_type = source.dtype().name(); + // _data_type_code = GetDataTypeCode(_data_type); + + // TENSORSTORE_CHECK_OK_AND_ASSIGN(___mb_cur_max, tensorstore::Open( + // GetZarrSpecToWrite(_filename, _image_shape, _chunk_shape, GetEncodedType(_dtype_code)), + // tensorstore::OpenMode::create | + // tensorstore::OpenMode::delete_existing, + // tensorstore::ReadWriteMode::write).result() + // ); +} + +// Private Methods +void PyramidCompositor::create_xml() { + CreateXML( + this->_output_pyramid_name, + this->_ome_metadata_file, + this->_plate_image_shapes[0], + this->_image_dtype + ); +} + +void PyramidCompositor::create_zattr_file() { + WriteTSZattrFilePlateImage( + this->_output_pyramid_name, + this->_output_pyramid_name + "/data.zarr/0", + this->_plate_image_shapes + ); +} + +void PyramidCompositor::create_zgroup_file() { + WriteVivZgroupFiles(this->_output_pyramid_name); +} + +void PyramidCompositor::create_auxiliary_files() { + create_xml(); + create_zattr_file(); + create_zgroup_file(); +} + +void PyramidCompositor::write_zarr_chunk(int level, int channel, int y_int, int x_int) { + // Implementation of tile data writing logic + // Placeholder logic + + // std::tuple y_range = { + // y_int * CHUNK_SIZE, + // std::min((y_int+1) * CHUNK_SIZE, this->zarr_arrays_[level].size()[3]) + // }; + + // std::tuple x_range = { + // x_int * CHUNK_SIZE, + // min((x_int+1) * CHUNK_SIZE, this->zarr_arrays_[level].size()[4]) + // }; + + // int assembled_width = std::get<1>(x_range) - std::get<0>(x_range); + // int assembled_height = std::get<1>(y_range) - std::get<0>(y_range); + + // std::vector<> + + // Compute ranges in global coordinates + int y_start = y_int * CHUNK_SIZE; + int y_end = std::min((y_int + 1) * CHUNK_SIZE, std::get<1>(_zarr_arrays[level])[0]); + int x_start = x_int * CHUNK_SIZE; + int x_end = std::min((x_int + 1) * CHUNK_SIZE, std::get<1>(_zarr_arrays[level])[1]); + + int assembled_width = x_end - x_start; + int assembled_height = y_end - y_start; + + image_data assembled_image = GetAssembledImageVector(assembled_height*assembled_width); + + // Determine required input images + int unit_image_height = _unit_image_shapes[level].first; + int unit_image_width = _unit_image_shapes[level].second; + + int row_start_pos = y_start; + while (row_start_pos < y_end) { + int row = row_start_pos / unit_image_height; + int local_y_start = row_start_pos - y_start; + int tile_y_start = row_start_pos - row * unit_image_height; + int tile_y_dim = std::min((row + 1) * unit_image_height - row_start_pos, y_end - row_start_pos); + int tile_y_end = tile_y_start + tile_y_dim; + + int col_start_pos = x_start; + while (col_start_pos < x_end) { + int col = col_start_pos / unit_image_width; + int local_x_start = col_start_pos - x_start; + int tile_x_start = col_start_pos - col * unit_image_width; + int tile_x_dim = std::min((col + 1) * unit_image_width - col_start_pos, x_end - col_start_pos); + int tile_x_end = tile_x_start + tile_x_dim; + + // Fetch input file path from composition map + std::string input_file_name = _composition_map[{col, row, channel}]; + std::filesystem::path zarrArrayLoc = std::filesystem::path(input_file_name) / "data.zarr/0/" / std::to_string(level); + + auto read_data = std::get<0>(ReadImageData(zarrArrayLoc, Seq(tile_y_start, tile_y_end), Seq(tile_x_start, tile_x_end))).get(); + + for (int i = local_y_start; i < local_y_start + tile_y_end - tile_y_start; ++i) { + for (int j = local_x_start; i < local_x_start + tile_x_end - tile_x_start; ++j) { + assembled_image[i*(tile_x_end - tile_x_start) + j] = read_data[(i-local_y_start)*(tile_x_end - tile_x_start) + (j-local_x_start)]; + } + } + + // // Open the Zarr array for reading with TensorStore + // auto zarrSpec = tensorstore::Context::Default() + // | tensorstore::OpenMode::read + // | tensorstore::dtype_v + // | tensorstore::Shape({1, 1, 1, assembledHeight, assembledWidth}); + + // auto inputZarrArray = tensorstore::Open(zarrSpec).result(); + + // TENSORSTORE_CHECK_OK_AND_ASSIGN() + + // if (!inputZarrArray.ok()) { + // std::cerr << "Failed to open Zarr array: " << inputZarrArray.status() << "\n"; + // return; + // } + + // // Read data from input Zarr file into assembledImage + // tensorstore::Read(*inputZarrArray, + // tensorstore::Span(assembledImage.data(), assembledImage.size())).result(); + + col_start_pos += tile_x_end - tile_x_start; + } + row_start_pos += tile_y_end - tile_y_start; + } + + + WriteImageData( + _input_pyramids_loc, + std::get<0>(_zarr_arrays[level]).get(), + Seq(y_start,y_end), Seq(x_start, x_end), + Seq(0,0), Seq(channel, channel), Seq(0,0) + ); +} + +void PyramidCompositor::setComposition(const std::unordered_map, std::string>& comp_map) { + _composition_map = comp_map; + + for (const auto& coord : comp_map) { + std::string file_path = coord.second; + std::filesystem::path attr_file_loc = std::filesystem::path(file_path) / "data.zarr/0/.zattrs"; + + if (std::filesystem::exists(attr_file_loc)) { + std::ifstream attr_file(attr_file_loc); + json attrs; + attr_file >> attrs; + + const auto& multiscale_metadata = attrs["multiscales"][0]["datasets"]; + _pyramid_levels = multiscale_metadata.size(); + for (const auto& dic : multiscale_metadata) { + std::string res_key = dic["path"]; + std::filesystem::path zarr_array_loc = std::filesystem::path(file_path) / "data.zarr/0/" / res_key; + + auto read_spec = GetZarrSpecToRead(zarr_array_loc); + + tensorstore::TensorStore source; + + TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( + read_spec, + tensorstore::OpenMode::open, + tensorstore::ReadWriteMode::read).result()); + + auto image_shape = source.domain().shape(); + _image_dtype = source.dtype().name(); + _data_type_code = GetDataTypeCode(_image_dtype); + + _unit_image_shapes[std::stoi(res_key)] = { + image_shape[image_shape.size()-2], + image_shape[image_shape.size()-1] + }; + } + break; + } + } + + int num_rows = 0, num_cols = 0, num_channels = 0; + for (const auto& coord : comp_map) { + num_rows = std::max(num_rows, std::get<1>(coord.first)); + num_cols = std::max(num_cols, std::get<0>(coord.first)); + _num_channels = std::max(num_channels, std::get<2>(coord.first)); + } + num_cols += 1; + num_rows += 1; + num_channels += 1; + _num_channels = num_channels; + + for (const auto& [level, shape] : _unit_image_shapes) { + std::string path = _output_pyramid_name + "/data.zarr/0/" + std::to_string(level); + + auto zarr_array_data = ReadImageData( + path, + Seq(0, num_rows * shape.first), + Seq(0, num_cols * shape.second), + Seq(0,0), + Seq(0, num_channels), + Seq(0,0) + ); + + _zarr_arrays[level] = zarr_array_data; + } +} + +void PyramidCompositor::reset_composition() { + // Remove pyramid file and clear internal data structures + fs::remove_all(_output_pyramid_name); + _composition_map.clear(); + _plate_image_shapes.clear(); + _tile_cache.clear(); + _plate_image_shapes.clear(); + _zarr_arrays.clear(); +} + +image_data PyramidCompositor::GetAssembledImageVector(int size){ + switch (_image_dtype_code) + { + case (1): + return std::vector(size); + case (2): + return std::vector(size); + case (4): + return std::vector(size); + case (8): + return std::vector(size); + case (16): + return std::vector(size); + case (32): + return std::vector(size); + case (64): + return std::vector(size); + case (128): + return std::vector(size); + case (256): + return std::vector(size); + case (512): + return std::vector(size); + default: + throw std::runtime_error("Invalid data type."); + } +} + +void PyramidCompositor::WriteImageData( + const std::string& path, + image_data& image, + const Seq& rows, + const Seq& cols, + const std::optional& layers, + const std::optional& channels, + const std::optional& tsteps) { + + std::vector shape; + + tensorstore::TensorStore source; + + TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( + GetZarrSpecToRead(path), + tensorstore::OpenMode::open, + tensorstore::ReadWriteMode::write).result()); + + auto output_transform = tensorstore::IdentityTransform(source.domain()); + + if (_t_int.has_value() && tsteps.has_value()) { + output_transform = (std::move(output_transform) | tensorstore::Dims(_t_int.value()).ClosedInterval(tsteps.value().Start(), tsteps.value().Stop())).value(); + shape.emplace_back(tsteps.value().Stop() - tsteps.value().Start()+1); + } + + if (_c_int.has_value() && channels.has_value()) { + output_transform = (std::move(output_transform) | tensorstore::Dims(_c_int.value()).ClosedInterval(channels.value().Start(), channels.value().Stop())).value(); + shape.emplace_back(channels.value().Stop() - channels.value().Start()+1); + } + + if (_z_int.has_value() && layers.has_value()) { + output_transform = (std::move(output_transform) | tensorstore::Dims(_z_int.value()).ClosedInterval(layers.value().Start(), layers.value().Stop())).value(); + shape.emplace_back(layers.value().Stop() - layers.value().Start()+1); + } + + output_transform = (std::move(output_transform) | tensorstore::Dims(_y_int).ClosedInterval(rows.Start(), rows.Stop()) | + tensorstore::Dims(_x_int).ClosedInterval(cols.Start(), cols.Stop())).value(); + + shape.emplace_back(rows.Stop() - rows.Start()+1); + shape.emplace_back(cols.Stop() - cols.Start()+1); + + auto data_array = tensorstore::Array(image.data(), shape, tensorstore::c_order); + + // Write data array to TensorStore + auto write_result = tensorstore::Write(tensorstore::UnownedToShared(data_array), source | output_transform).result(); + + if (!write_result.ok()) { + std::cerr << "Error writing image: " << write_result.status() << std::endl; + } +} + +template +std::shared_ptr> PyramidCompositor::GetImageDataTemplated(const Seq& rows, const Seq& cols, const Seq& layers, const Seq& channels, const Seq& tsteps){ + + const auto data_height = rows.Stop() - rows.Start() + 1; + const auto data_width = cols.Stop() - cols.Start() + 1; + const auto data_depth = layers.Stop() - layers.Start() + 1; + const auto data_num_channels = channels.Stop() - channels.Start() + 1; + const auto data_tsteps = tsteps.Stop() - tsteps.Start() + 1; + + auto read_buffer = std::make_shared>(data_height*data_width*data_depth*data_num_channels*data_tsteps); + auto array = tensorstore::Array(read_buffer->data(), {data_tsteps, data_num_channels, data_depth, data_height, data_width}, tensorstore::c_order); + tensorstore::intTransform<> read_transform = tensorstore::IdentityTransform(source.domain()); + + int x_int=1, y_int=0; + if (_t_int.has_value()){ + read_transform = (std::move(read_transform) | tensorstore::Dims(_t_int.value()).ClosedInterval(tsteps.Start(), tsteps.Stop())).value(); + x_int++; + y_int++; + } + if (_c_int.has_value()){ + read_transform = (std::move(read_transform) | tensorstore::Dims(_c_int.value()).ClosedInterval(channels.Start(), channels.Stop())).value(); + x_int++; + y_int++; + } + if (_z_int.has_value()){ + read_transform = (std::move(read_transform) | tensorstore::Dims(_z_int.value()).ClosedInterval(layers.Start(), layers.Stop())).value(); + x_int++; + y_int++; + } + read_transform = (std::move(read_transform) | tensorstore::Dims(y_int).ClosedInterval(rows.Start(), rows.Stop()) | + tensorstore::Dims(x_int).ClosedInterval(cols.Start(), cols.Stop())).value(); + + + tensorstore::Read(source | read_transform, tensorstore::UnownedToShared(array)).value(); + + return read_buffer; +} + +std::tuple, std::vector> PyramidCompositor::ReadImageData( + const std::string& fname, + const Seq& rows, const Seq& cols, + const Seq& layers = Seq(0,0), + const Seq& channels = Seq(0,0), + const Seq& tsteps = Seq(0,0)) { + + auto read_spec = GetZarrSpecToRead(fname); + + tensorstore::TensorStore source; + + TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( + read_spec, + tensorstore::OpenMode::open, + tensorstore::ReadWriteMode::read).result()); + + auto image_shape = source.domain().shape(); + const auto read_chunk_shape = source.chunk_layout().value().read_chunk_shape(); + + // if (image_shape.size() == 5){ + // _image_height = image_shape[3]; + // _image_width = image_shape[4]; + // _image_depth = image_shape[2]; + // _num_channels = image_shape[1]; + // _num_tsteps = image_shape[0]; + // _t_index.emplace(0); + // _c_index.emplace(1); + // _z_index.emplace(2); + // } else { + // assert(image_shape.size() >= 2); + // std::tie(_t_index, _c_index, _z_index) = ParseMultiscaleMetadata(axes_list, image_shape.size()); + // _image_height = image_shape[image_shape.size()-2]; + // _image_width = image_shape[image_shape.size()-1]; + // if (_t_index.has_value()) { + // _num_tsteps = image_shape[_t_index.value()]; + // } else { + // _num_tsteps = 1; + // } + + // if (_c_index.has_value()) { + // _num_channels = image_shape[_c_index.value()]; + // } else { + // _num_channels = 1; + // } + + // if (_z_index.has_value()) { + // _image_depth = image_shape[_z_index.value()]; + // } else { + // _image_depth = 1; + // } + // } + + _image_dtype = source.dtype().name(); + _data_type_code = GetDataTypeCode(_image_dtype); + + switch (_data_type_code) + { + case (1): + return std::make_tuple( + std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), + std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); + case (2): + return std::make_tuple( + std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), + std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); + case (4): + return std::make_tuple( + std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), + std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); + case (8): + return std::make_tuple( + std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), + std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); + case (16): + return std::make_tuple( + std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), + std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); + case (32): + return std::make_tuple( + std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), + std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); + case (64): + return std::make_tuple( + std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), + std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); + case (128): + return std::make_tuple( + std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), + std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); + case (256): + return std::make_tuple( + std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), + std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); + case (512): + return std::make_tuple( + std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), + std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); + default: + break; + } +} +} // ns argolid \ No newline at end of file diff --git a/src/cpp/core/pyramid_compositor.h b/src/cpp/core/pyramid_compositor.h new file mode 100644 index 0000000..10a7e4a --- /dev/null +++ b/src/cpp/core/pyramid_compositor.h @@ -0,0 +1,100 @@ +#pragma once +#include +#include +#include +#include +#include +#include + +#include "tensorstore/tensorstore.h" +#include "tensorstore/spec.h" + +namespace fs = std::filesystem; +namespace argolid { +static constexpr int CHUNK_SIZE = 1024; + +using image_data = std::variant, + std::vector, + std::vector, + std::vector, + std::vector, + std::vector, + std::vector, + std::vector, + std::vector, + std::vector>; + +class Seq +{ + private: + long start_index_, stop_index_, step_; + public: + inline Seq(const long start, const long stop, const long step=1):start_index_(start), stop_index_(stop), step_(step){} + inline long Start() const {return start_index_;} + inline long Stop() const {return stop_index_;} + inline long Step() const {return step_;} +}; + +class PyramidCompositor { +public: + + PyramidCompositor(const std::string& well_pyramid_loc, const std::string& out_dir, const std::string& pyramid_file_name); + + void reset_composition(); + void get_tile_data(int level, int channel, int y_index, int x_index); + + void create_xml(); + void create_zattr_file(); + void create_zgroup_file(); + void create_auxiliary_files(); + void write_zarr_chunk(int level, int channel, int y_index, int x_index); +private: + + template + void WriteImageData( + const std::string& path, + image_data& image, + const Seq& rows, + const Seq& cols, + const std::optional& layers, + const std::optional& channels, + const std::optional& tsteps); + + template + std::shared_ptr> GetImageDataTemplated( + const Seq& rows, + const Seq& cols, + const Seq& layers, + const Seq& channels, + const Seq& tsteps); + + std::tuple, std::vector> ReadImageData( + const std::string& fname, + const Seq& rows, const Seq& cols, + const Seq& layers = Seq(0,0), + const Seq& channels = Seq(0,0), + const Seq& tsteps = Seq(0,0)); + + image_data GetAssembledImageVector(int size); + + void setComposition(const std::unordered_map, std::string>& comp_map); + + // members for pyramid composition + std::string _input_pyramids_loc; + std::string _output_pyramid_name; + std::string _ome_metadata_file; + std::unordered_map> _plate_image_shapes; + // std::unordered_map> _zarr_arrays; + std::unordered_map, std::vector>> _zarr_arrays; // tuple at 0 is the image and at 1 is the size + + std::unordered_map> _unit_image_shapes; + std::unordered_map, std::string> _composition_map; + + std::set> _tile_cache; + int _pyramid_levels; + int _num_channels; + + std::string _image_dtype; + std::uint16_t _data_type_code; +}; +} // ns argolid diff --git a/src/cpp/interface/pyramid_python_interface.cpp b/src/cpp/interface/pyramid_python_interface.cpp index adefd7f..8a4031b 100644 --- a/src/cpp/interface/pyramid_python_interface.cpp +++ b/src/cpp/interface/pyramid_python_interface.cpp @@ -1,5 +1,6 @@ #include "../core/ome_tiff_to_chunked_pyramid.h" #include "../core/pyramid_view.h" +#include "../core/pyramid_compositor.h" #include "../utilities/utilities.h" #include #include @@ -17,6 +18,17 @@ PYBIND11_MODULE(libargolid, m) { .def("GeneratePyramid", &argolid::PyramidView::GeneratePyramid) \ .def("AssembleBaseLevel", &argolid::PyramidView::AssembleBaseLevel) ; + + py::class_(m, "PyramidCompositorCPP") \ + .def(py::init) \ + .def("set_well_map", &argolid::PyramidCompositor::set_well_map) \ + .def("reset_composition", &argolid::PyramidCompositor::reset_composition) \ + .def("get_tile_data", &argolid::PyramidCompositor::get_tile_data) \ + .def("create_xml", &argolid::PyramidCompositor::create_xml) \ + .def("create_zattr_file", &argolid::PyramidCompositor::create_zattr_file) \ + .def("create_auxiliary_files", &argolid::PyramidCompositor::create_auxiliary_files) \ + .def("write_zarr_chunk", &argolid::PyramidCompositor::write_zarr_chunk) ; + py::enum_(m, "VisType") .value("NG_Zarr", argolid::VisType::NG_Zarr) .value("PCNG", argolid::VisType::PCNG) diff --git a/src/cpp/utilities/utilities.cpp b/src/cpp/utilities/utilities.cpp index 45dd9cd..3f52fee 100644 --- a/src/cpp/utilities/utilities.cpp +++ b/src/cpp/utilities/utilities.cpp @@ -202,6 +202,51 @@ void WriteTSZattrFile(const std::string& tiff_file_name, const std::string& zarr } } +void WriteTSZattrFilePlateImage( + const std::string& tiff_file_name, + const std::string& zarr_root_dir, + const std::unordered_map>& plate_image_shape){ + + json zarr_multiscale_axes; + zarr_multiscale_axes = json::parse(R"([ + {"name": "c", "type": "channel"}, + {"name": "z", "type": "space", "unit": "micrometer"}, + {"name": "y", "type": "space", "unit": "micrometer"}, + {"name": "x", "type": "space", "unit": "micrometer"} + ])"); + + + float level = 1.0; + json scale_metadata_list = json::array(); + for(const auto& key_val: plate_image_shape){ + json scale_metadata; + scale_metadata["path"] = std::to_string(std::get<0>(key_val)); + scale_metadata["coordinateTransformations"] = {{{"type", "scale"}, {"scale", {1.0, 1.0, level, level}}}}; + scale_metadata_list.push_back(scale_metadata); + level = level*2; + } + + json combined_metadata; + combined_metadata["datasets"] = scale_metadata_list; + combined_metadata["version"] = "0.4"; + combined_metadata["axes"] = zarr_multiscale_axes; + combined_metadata["name"] = tiff_file_name; + combined_metadata["metadata"] = {{"method", "mean"}}; + json final_formated_metadata; +#if defined(__clang__) || defined(_MSC_VER) +// more details here: https://github.com/nlohmann/json/issues/2311 + final_formated_metadata["multiscales"][0] = {combined_metadata}; +#else + final_formated_metadata["multiscales"] = {combined_metadata}; +#endif + std::ofstream f(zarr_root_dir + "/.zattrs",std::ios_base::trunc |std::ios_base::out); + if (f.is_open()){ + f << final_formated_metadata; + } else { + PLOG_INFO <<"Unable to write .zattr file at "<< zarr_root_dir << "."; + } +} + void WriteVivZattrFile(const std::string& tiff_file_name, const std::string& zattr_file_loc, int min_level, int max_level){ json scale_metadata_list = json::array(); @@ -318,6 +363,55 @@ void GenerateOmeXML(const std::string& image_name, const std::string& output_fil doc.save_file(output_file.c_str()); } +void CreateXML( + const std::string& image_name, + const std::string& output_file, + const std::tuple& image_shape, + const std::string& image_dtype) { + + pugi::xml_document doc; + + // Create the root element + pugi::xml_node omeNode = doc.append_child("OME"); + + // Add the namespaces and attributes to the root element + omeNode.append_attribute("xmlns") = "http://www.openmicroscopy.org/Schemas/OME/2016-06"; + omeNode.append_attribute("xmlns:xsi") = "http://www.w3.org/2001/XMLSchema-instance"; + auto creator = std::string{"Argolid "} + std::string{VERSION_INFO}; + omeNode.append_attribute("Creator") = creator.c_str(); + omeNode.append_attribute("UUID") = "urn:uuid:ce3367ae-0512-4e87-a045-20d87db14001"; + omeNode.append_attribute("xsi:schemaLocation") = "http://www.openmicroscopy.org/Schemas/OME/2016-06 http://www.openmicroscopy.org/Schemas/OME/2016-06/ome.xsd"; + + // Create the element + pugi::xml_node imageNode = omeNode.append_child("Image"); + imageNode.append_attribute("ID") = "Image:0"; + imageNode.append_attribute("Name") = image_name.c_str(); + + // Create the element + pugi::xml_node pixelsNode = imageNode.append_child("Pixels"); + pixelsNode.append_attribute("BigEndian") = "false"; + pixelsNode.append_attribute("DimensionOrder") = "XYZCT"; + pixelsNode.append_attribute("ID") = "Pixels:0"; + pixelsNode.append_attribute("Interleaved") = "false"; + pixelsNode.append_attribute("SizeC") = std::to_string(std::get<1>(image_shape)).c_str(); + pixelsNode.append_attribute("SizeZ") = std::to_string(std::get<2>(image_shape)).c_str(); + pixelsNode.append_attribute("SizeT") = std::to_string(std::get<0>(image_shape)).c_str(); + pixelsNode.append_attribute("SizeX") = std::to_string(std::get<4>(image_shape)).c_str(); + pixelsNode.append_attribute("SizeY") = std::to_string(std::get<3>(image_shape)).c_str(); + pixelsNode.append_attribute("Type") = image_dtype.c_str(); + + // Create the elements + for(std::int64_t i=0; i(image_shape); ++i){ + pugi::xml_node channelNode = pixelsNode.append_child("Channel"); + channelNode.append_attribute("ID") = ("Channel:0:" + std::to_string(i)).c_str(); + channelNode.append_attribute("SamplesPerPixel") = "1"; + // Create the elements + channelNode.append_child("LightPath"); + } + + doc.save_file(output_file.c_str()); +} + void WriteMultiscaleMetadataForSingleFile( const std::string& input_file , const std::string& output_dir, int min_level, int max_level, VisType v) { diff --git a/src/cpp/utilities/utilities.h b/src/cpp/utilities/utilities.h index fb6fd36..1a65203 100644 --- a/src/cpp/utilities/utilities.h +++ b/src/cpp/utilities/utilities.h @@ -60,4 +60,16 @@ inline std::tuple GetZarrParams(VisType v){ } std::optional> GetTiffDims (const std::string filename); + +void CreateXML( + const std::string& image_name, + const std::string& output_file, + const std::tuple& image_shape, + const std::string& image_dtype); + +void WriteTSZattrFilePlateImage( + const std::string& tiff_file_name, + const std::string& zarr_root_dir, + const std::unordered_map>& plate_image_shape); + } // ns argolid \ No newline at end of file diff --git a/src/python/argolid/pyramid_compositor.py b/src/python/argolid/pyramid_compositor.py index 8d0c15d..37b31e7 100644 --- a/src/python/argolid/pyramid_compositor.py +++ b/src/python/argolid/pyramid_compositor.py @@ -8,6 +8,8 @@ import ome_types import tensorstore as ts +from .libargolid import PyramidCompositorCPP + CHUNK_SIZE: int = 1024 OME_DTYPE = { @@ -93,7 +95,7 @@ def get_zarr_write_spec( } -class PyramidCompositor: +class PyramidCompositorPY: """ A class for composing a group of pyramid images into an assembled pyramid structure. """ @@ -399,3 +401,99 @@ def get_zarr_chunk( self._write_zarr_chunk(level, channel, y_index, x_index) self._chunk_cache.add((level, channel, y_index, x_index)) return + +class PyramidCompositor2: + """ + A class for composing a group of pyramid images into an assembled pyramid structure. + """ + + def __init__( + self, input_pyramids_loc: str, out_dir: str, output_pyramid_name: str + ) -> None: + """ + Initializes the PyramidCompositor object. + + Args: + input_pyramids_loc (str): The location of the input pyramid images. + out_dir (str): The output directory for the composed zarr pyramid file. + output_pyramid_name (str): The name of the zarr pyramid file. + """ + self._pyramid_compositor_cpp = PyramidCompositorCPP(input_pyramids_loc, out_dir, output_pyramid_name) + + + def _create_xml(self) -> None: + """ + Writes an OME-XML metadata file for the pyramid. + """ + self._pyramid_compositor_cpp.create_xml() + + + def _create_zattr_file(self) -> None: + """ + Creates a .zattrs file for the zarr pyramid. + """ + self._pyramid_compositor_cpp.create_zattr_file() + + def _create_zgroup_file(self) -> None: + """ + Creates .zgroup files for the zarr pyramid. + """ + self._pyramid_compositor_cpp.create_zgroup_file() + + def _create_auxilary_files(self) -> None: + """ + Creates auxiliary files (OME-XML, .zattrs, .zgroup) for the pyramid. + """ + self._pyramid_compositor_cpp.create_auxilary_files() + + + def _write_zarr_chunk( + self, level: int, channel: int, y_index: int, x_index: int + ) -> None: + """ + Writes the chunk file at the specified level, channel, y_index, and x_index. + + Args: + level (int): The level of the pyramid. + channel (int): The channel of the pyramid. + y_index (int): The y-index of the tile. + x_index (int): The x-index of the tile. + """ + self._pyramid_compositor_cpp.write_zarr_chunk(level, channel, y_index, x_index) + + def set_composition(self, composition_map: dict) -> None: + """ + Sets the composition for the pyramid. + + Args: + composition_map (dict): A dictionary mapping composition images to file paths. + """ + self._pyramid_compositor_cpp.set_composition(composition_map) + + def reset_composition(self) -> None: + """ + Resets the pyramid composition by removing the pyramid file and clearing internal data structures. + """ + self._pyramid_compositor_cpp.reset_composition() + + def get_zarr_chunk( + self, level: int, channel: int, y_index: int, x_index: int + ) -> None: + """ + Retrieves zarr chunk data from the pyramid at the specified level, channel, y_index, and x_index. + + This method checks if the requested zarr chunk is already in the cache. If not, it calls + the `_write_zarr_chunk` method to generate the chunk and add it to the cache. + + Args: + level (int): The level of the pyramid. + channel (int): The channel of the pyramid. + y_index (int): The y-index of the tile. + x_index (int): The x-index of the tile. + + Raises: + ValueError: If the composition map is not set, the requested level does not exist, + the requested channel does not exist, or the requested y_index or x_index + is out of bounds. + """ + self._pyramid_compositor_cpp.get_zarr_chunk(level, channel, y_index, x_index) From c2e2e87bdd8eaf0747fef17bd89f588b3ed65933 Mon Sep 17 00:00:00 2001 From: JesseMckinzie Date: Tue, 14 Jan 2025 08:26:23 -0700 Subject: [PATCH 02/24] Update c++ pyramid compositor --- src/cpp/core/pyramid_compositor.cpp | 213 ++++---- src/cpp/core/pyramid_compositor.h | 32 +- .../core/pyramid_compositor_from_trash.cpp | 514 ++++++++++++++++++ src/cpp/core/pyramid_compositor_old.cpp | 514 ++++++++++++++++++ src/cpp/utilities/utilities.cpp | 25 + src/cpp/utilities/utilities.h | 2 + 6 files changed, 1193 insertions(+), 107 deletions(-) create mode 100644 src/cpp/core/pyramid_compositor_from_trash.cpp create mode 100644 src/cpp/core/pyramid_compositor_old.cpp diff --git a/src/cpp/core/pyramid_compositor.cpp b/src/cpp/core/pyramid_compositor.cpp index 7732fd1..08ae8a1 100644 --- a/src/cpp/core/pyramid_compositor.cpp +++ b/src/cpp/core/pyramid_compositor.cpp @@ -37,50 +37,26 @@ PyramidCompositor::PyramidCompositor(const std::string& input_pyramids_loc, cons tensorstore::ReadWriteMode::read).result()); _image_dtype = source.dtype().name(); - _data_type_code = GetDataTypeCode(_image_dtype); + _image_dtype_code = GetDataTypeCode(_image_dtype); - // auto read_spec = GetZarrSpecToRead(input_pyramids_loc) + auto read_spec = GetZarrSpecToRead(input_pyramids_loc); - // TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( - // read_spec, - // tensorstore::OpenMode::open, - // tensorstore::ReadWriteMode::read).result()); + TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( + read_spec, + tensorstore::OpenMode::open, + tensorstore::ReadWriteMode::read).result()); - // auto image_shape = source.domain().shape(); - // const auto read_chunk_shape = source.chunk_layout().value().read_chunk_shape(); - - // if (image_shape.size() == 5){ - // _image_height = image_shape[3]; - // _image_width = image_shape[4]; - // _image_depth = image_shape[2]; - // _num_channels = image_shape[1]; - // _num_tsteps = image_shape[0]; - // _t_int.emplace(0); - // _c_int.emplace(1); - // _z_int.emplace(2); - // } else { - // assert(image_shape.size() >= 2); - // std::tie(_t_int, _c_int, _z_int) = ParseMultiscaleMetadata(axes_list, image_shape.size()); - // _image_height = image_shape[image_shape.size()-2]; - // _image_width = image_shape[image_shape.size()-1]; - // if (_t_int.has_value()) { - // _num_tsteps = image_shape[_t_int.value()]; - // } else { - // _num_tsteps = 1; - // } - - // if (_c_int.has_value()) { - // _num_channels = image_shape[_c_int.value()]; - // } else { - // _num_channels = 1; - // } + auto image_shape = source.domain().shape(); + const auto read_chunk_shape = source.chunk_layout().value().read_chunk_shape(); - // if (_z_int.has_value()) { - // _image_depth = image_shape[_z_int.value()]; - // } else { - // _image_depth = 1; - // } - // } + if (image_shape.size() == 5){ + _t_index.emplace(0); + _c_index.emplace(1); + _z_index.emplace(2); + } else { + assert(image_shape.size() >= 2); + std::tie(_t_index, _c_index, _z_index) = ParseMultiscaleMetadata("XYCTZ", image_shape.size()); + } // _data_type = source.dtype().name(); // _data_type_code = GetDataTypeCode(_data_type); @@ -142,9 +118,9 @@ void PyramidCompositor::write_zarr_chunk(int level, int channel, int y_int, int // Compute ranges in global coordinates int y_start = y_int * CHUNK_SIZE; - int y_end = std::min((y_int + 1) * CHUNK_SIZE, std::get<1>(_zarr_arrays[level])[0]); + int y_end = std::min((y_int + 1) * CHUNK_SIZE, (int)std::get<1>(_zarr_arrays[level])[0]); int x_start = x_int * CHUNK_SIZE; - int x_end = std::min((x_int + 1) * CHUNK_SIZE, std::get<1>(_zarr_arrays[level])[1]); + int x_end = std::min((x_int + 1) * CHUNK_SIZE, (int)std::get<1>(_zarr_arrays[level])[1]); int assembled_width = x_end - x_start; int assembled_height = y_end - y_start; @@ -175,13 +151,28 @@ void PyramidCompositor::write_zarr_chunk(int level, int channel, int y_int, int std::string input_file_name = _composition_map[{col, row, channel}]; std::filesystem::path zarrArrayLoc = std::filesystem::path(input_file_name) / "data.zarr/0/" / std::to_string(level); - auto read_data = std::get<0>(ReadImageData(zarrArrayLoc, Seq(tile_y_start, tile_y_end), Seq(tile_x_start, tile_x_end))).get(); + auto read_data = *std::get<0>(ReadImageData(zarrArrayLoc, Seq(tile_y_start, tile_y_end), Seq(tile_x_start, tile_x_end))).get(); - for (int i = local_y_start; i < local_y_start + tile_y_end - tile_y_start; ++i) { - for (int j = local_x_start; i < local_x_start + tile_x_end - tile_x_start; ++j) { - assembled_image[i*(tile_x_end - tile_x_start) + j] = read_data[(i-local_y_start)*(tile_x_end - tile_x_start) + (j-local_x_start)]; + // Use std::visit to handle the different types in the variant + std::visit([&](auto& assembled_vec) { + using T = typename std::decay_t::value_type; + + for (int i = local_y_start; i < local_y_start + tile_y_end - tile_y_start; ++i) { + for (int j = local_x_start; j < local_x_start + tile_x_end - tile_x_start; ++j) { + // Use std::visit to extract the value from read_data + std::visit([&](auto& value) { + assembled_vec[i * assembled_width + j] = static_cast(value[(i - local_y_start) * (tile_x_end - tile_x_start) + (j - local_x_start)]); + }, read_data); + } } - } + }, assembled_image); + + + // for (int i = local_y_start; i < local_y_start + tile_y_end - tile_y_start; ++i) { + // for (int j = local_x_start; i < local_x_start + tile_x_end - tile_x_start; ++j) { + // assembled_image[i*(tile_x_end - tile_x_start) + j] = read_data[(i-local_y_start)*(tile_x_end - tile_x_start) + (j-local_x_start)]; + // } + // } // // Open the Zarr array for reading with TensorStore // auto zarrSpec = tensorstore::Context::Default() @@ -207,16 +198,34 @@ void PyramidCompositor::write_zarr_chunk(int level, int channel, int y_int, int row_start_pos += tile_y_end - tile_y_start; } - - WriteImageData( - _input_pyramids_loc, - std::get<0>(_zarr_arrays[level]).get(), - Seq(y_start,y_end), Seq(x_start, x_end), - Seq(0,0), Seq(channel, channel), Seq(0,0) - ); + auto data = *(std::get<0>(_zarr_arrays[level]).get()); + // WriteImageData( + // _input_pyramids_loc, + // data, + // Seq(y_start,y_end), Seq(x_start, x_end), + // Seq(0,0), Seq(channel, channel), Seq(0,0) + // ); + + // Assuming _zarr_arrays[level] is of type std::tuple>, std::vector> + //auto& [variant_ptr, some_vector] = _zarr_arrays[level]; // Structured binding (C++17) + + std::visit([&](const auto& array) { + using T = std::decay_t; + + // Use external variables and call WriteImageData + WriteImageData( + _input_pyramids_loc, + array, // Access the data from the std::shared_ptr + Seq(y_start, y_end), + Seq(x_start, x_end), + Seq(0, 0), + Seq(channel, channel), + Seq(0, 0) + ); + }, data); // Dereference the shared_ptr to get the variant } -void PyramidCompositor::setComposition(const std::unordered_map, std::string>& comp_map) { +void PyramidCompositor::setComposition(const std::unordered_map, std::string, TupleHash>& comp_map) { _composition_map = comp_map; for (const auto& coord : comp_map) { @@ -245,7 +254,7 @@ void PyramidCompositor::setComposition(const std::unordered_map void PyramidCompositor::WriteImageData( const std::string& path, - image_data& image, + std::vector& image, const Seq& rows, const Seq& cols, const std::optional& layers, @@ -341,23 +352,23 @@ void PyramidCompositor::WriteImageData( auto output_transform = tensorstore::IdentityTransform(source.domain()); - if (_t_int.has_value() && tsteps.has_value()) { - output_transform = (std::move(output_transform) | tensorstore::Dims(_t_int.value()).ClosedInterval(tsteps.value().Start(), tsteps.value().Stop())).value(); + if (_t_index.has_value() && tsteps.has_value()) { + output_transform = (std::move(output_transform) | tensorstore::Dims(_t_index.value()).ClosedInterval(tsteps.value().Start(), tsteps.value().Stop())).value(); shape.emplace_back(tsteps.value().Stop() - tsteps.value().Start()+1); } - if (_c_int.has_value() && channels.has_value()) { - output_transform = (std::move(output_transform) | tensorstore::Dims(_c_int.value()).ClosedInterval(channels.value().Start(), channels.value().Stop())).value(); + if (_c_index.has_value() && channels.has_value()) { + output_transform = (std::move(output_transform) | tensorstore::Dims(_c_index.value()).ClosedInterval(channels.value().Start(), channels.value().Stop())).value(); shape.emplace_back(channels.value().Stop() - channels.value().Start()+1); } - if (_z_int.has_value() && layers.has_value()) { - output_transform = (std::move(output_transform) | tensorstore::Dims(_z_int.value()).ClosedInterval(layers.value().Start(), layers.value().Stop())).value(); + if (_z_index.has_value() && layers.has_value()) { + output_transform = (std::move(output_transform) | tensorstore::Dims(_z_index.value()).ClosedInterval(layers.value().Start(), layers.value().Stop())).value(); shape.emplace_back(layers.value().Stop() - layers.value().Start()+1); } - output_transform = (std::move(output_transform) | tensorstore::Dims(_y_int).ClosedInterval(rows.Start(), rows.Stop()) | - tensorstore::Dims(_x_int).ClosedInterval(cols.Start(), cols.Stop())).value(); + output_transform = (std::move(output_transform) | tensorstore::Dims(_y_index).ClosedInterval(rows.Start(), rows.Stop()) | + tensorstore::Dims(_x_index).ClosedInterval(cols.Start(), cols.Stop())).value(); shape.emplace_back(rows.Stop() - rows.Start()+1); shape.emplace_back(cols.Stop() - cols.Start()+1); @@ -373,7 +384,13 @@ void PyramidCompositor::WriteImageData( } template -std::shared_ptr> PyramidCompositor::GetImageDataTemplated(const Seq& rows, const Seq& cols, const Seq& layers, const Seq& channels, const Seq& tsteps){ +std::shared_ptr> PyramidCompositor::GetImageDataTemplated( + tensorstore::TensorStore& source, + const Seq& rows, + const Seq& cols, + const Seq& layers, + const Seq& channels, + const Seq& tsteps) { const auto data_height = rows.Stop() - rows.Start() + 1; const auto data_width = cols.Stop() - cols.Start() + 1; @@ -383,21 +400,21 @@ std::shared_ptr> PyramidCompositor::GetImageDataTemplated(const S auto read_buffer = std::make_shared>(data_height*data_width*data_depth*data_num_channels*data_tsteps); auto array = tensorstore::Array(read_buffer->data(), {data_tsteps, data_num_channels, data_depth, data_height, data_width}, tensorstore::c_order); - tensorstore::intTransform<> read_transform = tensorstore::IdentityTransform(source.domain()); + tensorstore::IndexTransform<> read_transform = tensorstore::IdentityTransform(source.domain()); int x_int=1, y_int=0; - if (_t_int.has_value()){ - read_transform = (std::move(read_transform) | tensorstore::Dims(_t_int.value()).ClosedInterval(tsteps.Start(), tsteps.Stop())).value(); + if (_t_index.has_value()){ + read_transform = (std::move(read_transform) | tensorstore::Dims(_t_index.value()).ClosedInterval(tsteps.Start(), tsteps.Stop())).value(); x_int++; y_int++; } - if (_c_int.has_value()){ - read_transform = (std::move(read_transform) | tensorstore::Dims(_c_int.value()).ClosedInterval(channels.Start(), channels.Stop())).value(); + if (_c_index.has_value()){ + read_transform = (std::move(read_transform) | tensorstore::Dims(_c_index.value()).ClosedInterval(channels.Start(), channels.Stop())).value(); x_int++; y_int++; } - if (_z_int.has_value()){ - read_transform = (std::move(read_transform) | tensorstore::Dims(_z_int.value()).ClosedInterval(layers.Start(), layers.Stop())).value(); + if (_z_index.has_value()){ + read_transform = (std::move(read_transform) | tensorstore::Dims(_z_index.value()).ClosedInterval(layers.Start(), layers.Stop())).value(); x_int++; y_int++; } @@ -410,12 +427,12 @@ std::shared_ptr> PyramidCompositor::GetImageDataTemplated(const S return read_buffer; } -std::tuple, std::vector> PyramidCompositor::ReadImageData( +std::tuple, std::vector> PyramidCompositor::ReadImageData( const std::string& fname, const Seq& rows, const Seq& cols, - const Seq& layers = Seq(0,0), - const Seq& channels = Seq(0,0), - const Seq& tsteps = Seq(0,0)) { + const Seq& layers, + const Seq& channels, + const Seq& tsteps) { auto read_spec = GetZarrSpecToRead(fname); @@ -463,50 +480,50 @@ std::tuple, std::vector> PyramidCompositor::Rea // } _image_dtype = source.dtype().name(); - _data_type_code = GetDataTypeCode(_image_dtype); + _image_dtype_code = GetDataTypeCode(_image_dtype); switch (_data_type_code) { case (1): return std::make_tuple( - std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), - std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); + std::make_shared(std::move(*(GetImageDataTemplated(source, rows, cols, layers, channels, tsteps)))), + std::vector{image_shape[image_shape.size()-2], image_shape[image_shape.size()-1]}); case (2): return std::make_tuple( - std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), - std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); + std::make_shared(std::move(*(GetImageDataTemplated(source, rows, cols, layers, channels, tsteps)))), + std::vector{image_shape[image_shape.size()-2], image_shape[image_shape.size()-1]}); case (4): return std::make_tuple( - std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), - std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); + std::make_shared(std::move(*(GetImageDataTemplated(source, rows, cols, layers, channels, tsteps)))), + std::vector{image_shape[image_shape.size()-2], image_shape[image_shape.size()-1]}); case (8): return std::make_tuple( - std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), - std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); + std::make_shared(std::move(*(GetImageDataTemplated(source, rows, cols, layers, channels, tsteps)))), + std::vector{image_shape[image_shape.size()-2], image_shape[image_shape.size()-1]}); case (16): return std::make_tuple( - std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), - std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); + std::make_shared(std::move(*(GetImageDataTemplated(source, rows, cols, layers, channels, tsteps)))), + std::vector{image_shape[image_shape.size()-2], image_shape[image_shape.size()-1]}); case (32): return std::make_tuple( - std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), - std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); + std::make_shared(std::move(*(GetImageDataTemplated(source, rows, cols, layers, channels, tsteps)))), + std::vector{image_shape[image_shape.size()-2], image_shape[image_shape.size()-1]}); case (64): return std::make_tuple( - std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), - std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); + std::make_shared(std::move(*(GetImageDataTemplated(source, rows, cols, layers, channels, tsteps)))), + std::vector{image_shape[image_shape.size()-2], image_shape[image_shape.size()-1]}); case (128): return std::make_tuple( - std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), - std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); + std::make_shared(std::move(*(GetImageDataTemplated(source, rows, cols, layers, channels, tsteps)))), + std::vector{image_shape[image_shape.size()-2], image_shape[image_shape.size()-1]}); case (256): return std::make_tuple( - std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), - std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); + std::make_shared(std::move(*(GetImageDataTemplated(source, rows, cols, layers, channels, tsteps)))), + std::vector{image_shape[image_shape.size()-2], image_shape[image_shape.size()-1]}); case (512): return std::make_tuple( - std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), - std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); + std::make_shared(std::move(*(GetImageDataTemplated(source, rows, cols, layers, channels, tsteps)))), + std::vector{image_shape[image_shape.size()-2], image_shape[image_shape.size()-1]}); default: break; } diff --git a/src/cpp/core/pyramid_compositor.h b/src/cpp/core/pyramid_compositor.h index 10a7e4a..860e748 100644 --- a/src/cpp/core/pyramid_compositor.h +++ b/src/cpp/core/pyramid_compositor.h @@ -13,6 +13,16 @@ namespace fs = std::filesystem; namespace argolid { static constexpr int CHUNK_SIZE = 1024; +// custom hashing for using std::tuple as key +struct TupleHash { + template + std::size_t operator()(const std::tuple& t) const { + return std::apply([](const auto&... args) { + return (std::hash>{}(args) ^ ...); + }, t); + } +}; + using image_data = std::variant, std::vector, std::vector, @@ -41,7 +51,6 @@ class PyramidCompositor { PyramidCompositor(const std::string& well_pyramid_loc, const std::string& out_dir, const std::string& pyramid_file_name); void reset_composition(); - void get_tile_data(int level, int channel, int y_index, int x_index); void create_xml(); void create_zattr_file(); @@ -49,11 +58,11 @@ class PyramidCompositor { void create_auxiliary_files(); void write_zarr_chunk(int level, int channel, int y_index, int x_index); private: - + template void WriteImageData( const std::string& path, - image_data& image, + std::vector& image, const Seq& rows, const Seq& cols, const std::optional& layers, @@ -62,22 +71,24 @@ class PyramidCompositor { template std::shared_ptr> GetImageDataTemplated( + tensorstore::TensorStore& source, const Seq& rows, const Seq& cols, const Seq& layers, const Seq& channels, const Seq& tsteps); - std::tuple, std::vector> ReadImageData( + std::tuple, std::vector> ReadImageData( const std::string& fname, const Seq& rows, const Seq& cols, const Seq& layers = Seq(0,0), const Seq& channels = Seq(0,0), - const Seq& tsteps = Seq(0,0)); + const Seq& tsteps = Seq(0,0) + ); image_data GetAssembledImageVector(int size); - void setComposition(const std::unordered_map, std::string>& comp_map); + void setComposition(const std::unordered_map, std::string, TupleHash>& comp_map); // members for pyramid composition std::string _input_pyramids_loc; @@ -85,16 +96,19 @@ class PyramidCompositor { std::string _ome_metadata_file; std::unordered_map> _plate_image_shapes; // std::unordered_map> _zarr_arrays; - std::unordered_map, std::vector>> _zarr_arrays; // tuple at 0 is the image and at 1 is the size + std::unordered_map, std::vector>> _zarr_arrays; // tuple at 0 is the image and at 1 is the size std::unordered_map> _unit_image_shapes; - std::unordered_map, std::string> _composition_map; + std::unordered_map, std::string, TupleHash> _composition_map; std::set> _tile_cache; int _pyramid_levels; int _num_channels; std::string _image_dtype; - std::uint16_t _data_type_code; + std::uint16_t _image_dtype_code; + + std::optional_z_index, _c_index, _t_index; + int _x_index, _y_index; }; } // ns argolid diff --git a/src/cpp/core/pyramid_compositor_from_trash.cpp b/src/cpp/core/pyramid_compositor_from_trash.cpp new file mode 100644 index 0000000..444db5a --- /dev/null +++ b/src/cpp/core/pyramid_compositor_from_trash.cpp @@ -0,0 +1,514 @@ +#include "pyramid_compositor.h" +#include "../utilities/utilities.h" + +#include +#include +#include +#include + +#include "tensorstore/context.h" +#include "tensorstore/array.h" +#include "tensorstore/driver/zarr/dtype.h" +#include "tensorstore/index_space/dim_expression.h" +#include "tensorstore/kvstore/kvstore.h" +#include "tensorstore/open.h" + +#include + +using ::tensorstore::internal_zarr::ChooseBaseDType; +using json = nlohmann::json; + +namespace fs = std::filesystem; + +namespace argolid { +PyramidCompositor::PyramidCompositor(const std::string& input_pyramids_loc, const std::string& out_dir, const std::string& output_pyramid_name): + _input_pyramids_loc(input_pyramids_loc), + _output_pyramid_name(out_dir + "/" + output_pyramid_name), + _ome_metadata_file(out_dir + "/" + output_pyramid_name + "/METADATA.ome.xml"), + _pyramid_levels(-1), + _num_channels(-1) { + + + tensorstore::TensorStore source; + + TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( + GetZarrSpecToRead(input_pyramids_loc), + tensorstore::OpenMode::open, + tensorstore::ReadWriteMode::read).result()); + + _image_dtype = source.dtype().name(); + _data_type_code = GetDataTypeCode(_image_dtype); + + // auto read_spec = GetZarrSpecToRead(input_pyramids_loc) + + // TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( + // read_spec, + // tensorstore::OpenMode::open, + // tensorstore::ReadWriteMode::read).result()); + + // auto image_shape = source.domain().shape(); + // const auto read_chunk_shape = source.chunk_layout().value().read_chunk_shape(); + + // if (image_shape.size() == 5){ + // _image_height = image_shape[3]; + // _image_width = image_shape[4]; + // _image_depth = image_shape[2]; + // _num_channels = image_shape[1]; + // _num_tsteps = image_shape[0]; + // _t_int.emplace(0); + // _c_int.emplace(1); + // _z_int.emplace(2); + // } else { + // assert(image_shape.size() >= 2); + // std::tie(_t_int, _c_int, _z_int) = ParseMultiscaleMetadata(axes_list, image_shape.size()); + // _image_height = image_shape[image_shape.size()-2]; + // _image_width = image_shape[image_shape.size()-1]; + // if (_t_int.has_value()) { + // _num_tsteps = image_shape[_t_int.value()]; + // } else { + // _num_tsteps = 1; + // } + + // if (_c_int.has_value()) { + // _num_channels = image_shape[_c_int.value()]; + // } else { + // _num_channels = 1; + // } + + // if (_z_int.has_value()) { + // _image_depth = image_shape[_z_int.value()]; + // } else { + // _image_depth = 1; + // } + // } + + // _data_type = source.dtype().name(); + // _data_type_code = GetDataTypeCode(_data_type); + + // TENSORSTORE_CHECK_OK_AND_ASSIGN(___mb_cur_max, tensorstore::Open( + // GetZarrSpecToWrite(_filename, _image_shape, _chunk_shape, GetEncodedType(_dtype_code)), + // tensorstore::OpenMode::create | + // tensorstore::OpenMode::delete_existing, + // tensorstore::ReadWriteMode::write).result() + // ); +} + +// Private Methods +void PyramidCompositor::create_xml() { + CreateXML( + this->_output_pyramid_name, + this->_ome_metadata_file, + this->_plate_image_shapes[0], + this->_image_dtype + ); +} + +void PyramidCompositor::create_zattr_file() { + WriteTSZattrFilePlateImage( + this->_output_pyramid_name, + this->_output_pyramid_name + "/data.zarr/0", + this->_plate_image_shapes + ); +} + +void PyramidCompositor::create_zgroup_file() { + WriteVivZgroupFiles(this->_output_pyramid_name); +} + +void PyramidCompositor::create_auxiliary_files() { + create_xml(); + create_zattr_file(); + create_zgroup_file(); +} + +void PyramidCompositor::write_zarr_chunk(int level, int channel, int y_int, int x_int) { + // Implementation of tile data writing logic + // Placeholder logic + + // std::tuple y_range = { + // y_int * CHUNK_SIZE, + // std::min((y_int+1) * CHUNK_SIZE, this->zarr_arrays_[level].size()[3]) + // }; + + // std::tuple x_range = { + // x_int * CHUNK_SIZE, + // min((x_int+1) * CHUNK_SIZE, this->zarr_arrays_[level].size()[4]) + // }; + + // int assembled_width = std::get<1>(x_range) - std::get<0>(x_range); + // int assembled_height = std::get<1>(y_range) - std::get<0>(y_range); + + // std::vector<> + + // Compute ranges in global coordinates + int y_start = y_int * CHUNK_SIZE; + int y_end = std::min((y_int + 1) * CHUNK_SIZE, std::get<1>(_zarr_arrays[level])[0]); + int x_start = x_int * CHUNK_SIZE; + int x_end = std::min((x_int + 1) * CHUNK_SIZE, std::get<1>(_zarr_arrays[level])[1]); + + int assembled_width = x_end - x_start; + int assembled_height = y_end - y_start; + + image_data assembled_image = GetAssembledImageVector(assembled_height*assembled_width); + + // Determine required input images + int unit_image_height = _unit_image_shapes[level].first; + int unit_image_width = _unit_image_shapes[level].second; + + int row_start_pos = y_start; + while (row_start_pos < y_end) { + int row = row_start_pos / unit_image_height; + int local_y_start = row_start_pos - y_start; + int tile_y_start = row_start_pos - row * unit_image_height; + int tile_y_dim = std::min((row + 1) * unit_image_height - row_start_pos, y_end - row_start_pos); + int tile_y_end = tile_y_start + tile_y_dim; + + int col_start_pos = x_start; + while (col_start_pos < x_end) { + int col = col_start_pos / unit_image_width; + int local_x_start = col_start_pos - x_start; + int tile_x_start = col_start_pos - col * unit_image_width; + int tile_x_dim = std::min((col + 1) * unit_image_width - col_start_pos, x_end - col_start_pos); + int tile_x_end = tile_x_start + tile_x_dim; + + // Fetch input file path from composition map + std::string input_file_name = _composition_map[{col, row, channel}]; + std::filesystem::path zarrArrayLoc = std::filesystem::path(input_file_name) / "data.zarr/0/" / std::to_string(level); + + auto read_data = std::get<0>(ReadImageData(zarrArrayLoc, Seq(tile_y_start, tile_y_end), Seq(tile_x_start, tile_x_end))).get(); + + for (int i = local_y_start; i < local_y_start + tile_y_end - tile_y_start; ++i) { + for (int j = local_x_start; i < local_x_start + tile_x_end - tile_x_start; ++j) { + assembled_image[i*(tile_x_end - tile_x_start) + j] = read_data[(i-local_y_start)*(tile_x_end - tile_x_start) + (j-local_x_start)]; + } + } + + // // Open the Zarr array for reading with TensorStore + // auto zarrSpec = tensorstore::Context::Default() + // | tensorstore::OpenMode::read + // | tensorstore::dtype_v + // | tensorstore::Shape({1, 1, 1, assembledHeight, assembledWidth}); + + // auto inputZarrArray = tensorstore::Open(zarrSpec).result(); + + // TENSORSTORE_CHECK_OK_AND_ASSIGN() + + // if (!inputZarrArray.ok()) { + // std::cerr << "Failed to open Zarr array: " << inputZarrArray.status() << "\n"; + // return; + // } + + // // Read data from input Zarr file into assembledImage + // tensorstore::Read(*inputZarrArray, + // tensorstore::Span(assembledImage.data(), assembledImage.size())).result(); + + col_start_pos += tile_x_end - tile_x_start; + } + row_start_pos += tile_y_end - tile_y_start; + } + + + WriteImageData( + _input_pyramids_loc, + std::get<0>(_zarr_arrays[level]).get(), + Seq(y_start,y_end), Seq(x_start, x_end), + Seq(0,0), Seq(channel, channel), Seq(0,0) + ); +} + +void PyramidCompositor::setComposition(const std::unordered_map, std::string>& comp_map) { + _composition_map = comp_map; + + for (const auto& coord : comp_map) { + std::string file_path = coord.second; + std::filesystem::path attr_file_loc = std::filesystem::path(file_path) / "data.zarr/0/.zattrs"; + + if (std::filesystem::exists(attr_file_loc)) { + std::ifstream attr_file(attr_file_loc); + json attrs; + attr_file >> attrs; + + const auto& multiscale_metadata = attrs["multiscales"][0]["datasets"]; + _pyramid_levels = multiscale_metadata.size(); + for (const auto& dic : multiscale_metadata) { + std::string res_key = dic["path"]; + std::filesystem::path zarr_array_loc = std::filesystem::path(file_path) / "data.zarr/0/" / res_key; + + auto read_spec = GetZarrSpecToRead(zarr_array_loc); + + tensorstore::TensorStore source; + + TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( + read_spec, + tensorstore::OpenMode::open, + tensorstore::ReadWriteMode::read).result()); + + auto image_shape = source.domain().shape(); + _image_dtype = source.dtype().name(); + _data_type_code = GetDataTypeCode(_image_dtype); + + _unit_image_shapes[std::stoi(res_key)] = { + image_shape[image_shape.size()-2], + image_shape[image_shape.size()-1] + }; + } + break; + } + } + + int num_rows = 0, num_cols = 0, num_channels = 0; + for (const auto& coord : comp_map) { + num_rows = std::max(num_rows, std::get<1>(coord.first)); + num_cols = std::max(num_cols, std::get<0>(coord.first)); + _num_channels = std::max(num_channels, std::get<2>(coord.first)); + } + num_cols += 1; + num_rows += 1; + num_channels += 1; + _num_channels = num_channels; + + for (const auto& [level, shape] : _unit_image_shapes) { + std::string path = _output_pyramid_name + "/data.zarr/0/" + std::to_string(level); + + auto zarr_array_data = ReadImageData( + path, + Seq(0, num_rows * shape.first), + Seq(0, num_cols * shape.second), + Seq(0,0), + Seq(0, num_channels), + Seq(0,0) + ); + + _zarr_arrays[level] = zarr_array_data; + } +} + +void PyramidCompositor::reset_composition() { + // Remove pyramid file and clear internal data structures + fs::remove_all(_output_pyramid_name); + _composition_map.clear(); + _plate_image_shapes.clear(); + _tile_cache.clear(); + _plate_image_shapes.clear(); + _zarr_arrays.clear(); +} + +image_data PyramidCompositor::GetAssembledImageVector(int size){ + switch (_image_dtype_code) + { + case (1): + return std::vector(size); + case (2): + return std::vector(size); + case (4): + return std::vector(size); + case (8): + return std::vector(size); + case (16): + return std::vector(size); + case (32): + return std::vector(size); + case (64): + return std::vector(size); + case (128): + return std::vector(size); + case (256): + return std::vector(size); + case (512): + return std::vector(size); + default: + throw std::runtime_error("Invalid data type."); + } +} + +void PyramidCompositor::WriteImageData( + const std::string& path, + image_data& image, + const Seq& rows, + const Seq& cols, + const std::optional& layers, + const std::optional& channels, + const std::optional& tsteps) { + + std::vector shape; + + tensorstore::TensorStore source; + + TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( + GetZarrSpecToRead(path), + tensorstore::OpenMode::open, + tensorstore::ReadWriteMode::write).result()); + + auto output_transform = tensorstore::IdentityTransform(source.domain()); + + if (_t_index.has_value() && tsteps.has_value()) { + output_transform = (std::move(output_transform) | tensorstore::Dims(_t_int.value()).ClosedInterval(tsteps.value().Start(), tsteps.value().Stop())).value(); + shape.emplace_back(tsteps.value().Stop() - tsteps.value().Start()+1); + } + + if (_c_int.has_value() && channels.has_value()) { + output_transform = (std::move(output_transform) | tensorstore::Dims(_c_int.value()).ClosedInterval(channels.value().Start(), channels.value().Stop())).value(); + shape.emplace_back(channels.value().Stop() - channels.value().Start()+1); + } + + if (_z_int.has_value() && layers.has_value()) { + output_transform = (std::move(output_transform) | tensorstore::Dims(_z_int.value()).ClosedInterval(layers.value().Start(), layers.value().Stop())).value(); + shape.emplace_back(layers.value().Stop() - layers.value().Start()+1); + } + + output_transform = (std::move(output_transform) | tensorstore::Dims(_y_int).ClosedInterval(rows.Start(), rows.Stop()) | + tensorstore::Dims(_x_int).ClosedInterval(cols.Start(), cols.Stop())).value(); + + shape.emplace_back(rows.Stop() - rows.Start()+1); + shape.emplace_back(cols.Stop() - cols.Start()+1); + + auto data_array = tensorstore::Array(image.data(), shape, tensorstore::c_order); + + // Write data array to TensorStore + auto write_result = tensorstore::Write(tensorstore::UnownedToShared(data_array), source | output_transform).result(); + + if (!write_result.ok()) { + std::cerr << "Error writing image: " << write_result.status() << std::endl; + } +} + +template +std::shared_ptr> PyramidCompositor::GetImageDataTemplated(const Seq& rows, const Seq& cols, const Seq& layers, const Seq& channels, const Seq& tsteps){ + + const auto data_height = rows.Stop() - rows.Start() + 1; + const auto data_width = cols.Stop() - cols.Start() + 1; + const auto data_depth = layers.Stop() - layers.Start() + 1; + const auto data_num_channels = channels.Stop() - channels.Start() + 1; + const auto data_tsteps = tsteps.Stop() - tsteps.Start() + 1; + + auto read_buffer = std::make_shared>(data_height*data_width*data_depth*data_num_channels*data_tsteps); + auto array = tensorstore::Array(read_buffer->data(), {data_tsteps, data_num_channels, data_depth, data_height, data_width}, tensorstore::c_order); + tensorstore::intTransform<> read_transform = tensorstore::IdentityTransform(source.domain()); + + int x_int=1, y_int=0; + if (_t_int.has_value()){ + read_transform = (std::move(read_transform) | tensorstore::Dims(_t_int.value()).ClosedInterval(tsteps.Start(), tsteps.Stop())).value(); + x_int++; + y_int++; + } + if (_c_int.has_value()){ + read_transform = (std::move(read_transform) | tensorstore::Dims(_c_int.value()).ClosedInterval(channels.Start(), channels.Stop())).value(); + x_int++; + y_int++; + } + if (_z_int.has_value()){ + read_transform = (std::move(read_transform) | tensorstore::Dims(_z_int.value()).ClosedInterval(layers.Start(), layers.Stop())).value(); + x_int++; + y_int++; + } + read_transform = (std::move(read_transform) | tensorstore::Dims(y_int).ClosedInterval(rows.Start(), rows.Stop()) | + tensorstore::Dims(x_int).ClosedInterval(cols.Start(), cols.Stop())).value(); + + + tensorstore::Read(source | read_transform, tensorstore::UnownedToShared(array)).value(); + + return read_buffer; +} + +std::tuple, std::vector> PyramidCompositor::ReadImageData( + const std::string& fname, + const Seq& rows, const Seq& cols, + const Seq& layers = Seq(0,0), + const Seq& channels = Seq(0,0), + const Seq& tsteps = Seq(0,0)) { + + auto read_spec = GetZarrSpecToRead(fname); + + tensorstore::TensorStore source; + + TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( + read_spec, + tensorstore::OpenMode::open, + tensorstore::ReadWriteMode::read).result()); + + auto image_shape = source.domain().shape(); + const auto read_chunk_shape = source.chunk_layout().value().read_chunk_shape(); + + // if (image_shape.size() == 5){ + // _image_height = image_shape[3]; + // _image_width = image_shape[4]; + // _image_depth = image_shape[2]; + // _num_channels = image_shape[1]; + // _num_tsteps = image_shape[0]; + // _t_index.emplace(0); + // _c_index.emplace(1); + // _z_index.emplace(2); + // } else { + // assert(image_shape.size() >= 2); + // std::tie(_t_index, _c_index, _z_index) = ParseMultiscaleMetadata(axes_list, image_shape.size()); + // _image_height = image_shape[image_shape.size()-2]; + // _image_width = image_shape[image_shape.size()-1]; + // if (_t_index.has_value()) { + // _num_tsteps = image_shape[_t_index.value()]; + // } else { + // _num_tsteps = 1; + // } + + // if (_c_index.has_value()) { + // _num_channels = image_shape[_c_index.value()]; + // } else { + // _num_channels = 1; + // } + + // if (_z_index.has_value()) { + // _image_depth = image_shape[_z_index.value()]; + // } else { + // _image_depth = 1; + // } + // } + + _image_dtype = source.dtype().name(); + _data_type_code = GetDataTypeCode(_image_dtype); + + switch (_data_type_code) + { + case (1): + return std::make_tuple( + std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), + std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); + case (2): + return std::make_tuple( + std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), + std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); + case (4): + return std::make_tuple( + std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), + std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); + case (8): + return std::make_tuple( + std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), + std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); + case (16): + return std::make_tuple( + std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), + std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); + case (32): + return std::make_tuple( + std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), + std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); + case (64): + return std::make_tuple( + std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), + std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); + case (128): + return std::make_tuple( + std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), + std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); + case (256): + return std::make_tuple( + std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), + std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); + case (512): + return std::make_tuple( + std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), + std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); + default: + break; + } +} +} // ns argolid \ No newline at end of file diff --git a/src/cpp/core/pyramid_compositor_old.cpp b/src/cpp/core/pyramid_compositor_old.cpp new file mode 100644 index 0000000..40fa5ff --- /dev/null +++ b/src/cpp/core/pyramid_compositor_old.cpp @@ -0,0 +1,514 @@ +#include "pyramid_compositor.h" +#include "../utilities/utilities.h" + +#include +#include +#include +#include + +#include "tensorstore/context.h" +#include "tensorstore/array.h" +#include "tensorstore/driver/zarr/dtype.h" +#include "tensorstore/index_space/dim_expression.h" +#include "tensorstore/kvstore/kvstore.h" +#include "tensorstore/open.h" + +#include + +using ::tensorstore::internal_zarr::ChooseBaseDType; +using json = nlohmann::json; + +namespace fs = std::filesystem; + +namespace argolid { +PyramidCompositor::PyramidCompositor(const std::string& input_pyramids_loc, const std::string& out_dir, const std::string& output_pyramid_name): + _input_pyramids_loc(input_pyramids_loc), + _output_pyramid_name(out_dir + "/" + output_pyramid_name), + _ome_metadata_file(out_dir + "/" + output_pyramid_name + "/METADATA.ome.xml"), + _pyramid_levels(-1), + _num_channels(-1) { + + + tensorstore::TensorStore source; + + TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( + GetZarrSpecToRead(input_pyramids_loc), + tensorstore::OpenMode::open, + tensorstore::ReadWriteMode::read).result()); + + _image_dtype = source.dtype().name(); + _data_type_code = GetDataTypeCode(_image_dtype); + + // auto read_spec = GetZarrSpecToRead(input_pyramids_loc) + + // TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( + // read_spec, + // tensorstore::OpenMode::open, + // tensorstore::ReadWriteMode::read).result()); + + // auto image_shape = source.domain().shape(); + // const auto read_chunk_shape = source.chunk_layout().value().read_chunk_shape(); + + // if (image_shape.size() == 5){ + // _image_height = image_shape[3]; + // _image_width = image_shape[4]; + // _image_depth = image_shape[2]; + // _num_channels = image_shape[1]; + // _num_tsteps = image_shape[0]; + // _t_index.emplace(0); + // _c_index.emplace(1); + // _z_index.emplace(2); + // } else { + // assert(image_shape.size() >= 2); + // std::tie(_t_index, _c_index, _z_index) = ParseMultiscaleMetadata(axes_list, image_shape.size()); + // _image_height = image_shape[image_shape.size()-2]; + // _image_width = image_shape[image_shape.size()-1]; + // if (_t_index.has_value()) { + // _num_tsteps = image_shape[_t_index.value()]; + // } else { + // _num_tsteps = 1; + // } + + // if (_c_index.has_value()) { + // _num_channels = image_shape[_c_index.value()]; + // } else { + // _num_channels = 1; + // } + + // if (_z_index.has_value()) { + // _image_depth = image_shape[_z_index.value()]; + // } else { + // _image_depth = 1; + // } + // } + + // _data_type = source.dtype().name(); + // _data_type_code = GetDataTypeCode(_data_type); + + // TENSORSTORE_CHECK_OK_AND_ASSIGN(___mb_cur_max, tensorstore::Open( + // GetZarrSpecToWrite(_filename, _image_shape, _chunk_shape, GetEncodedType(_dtype_code)), + // tensorstore::OpenMode::create | + // tensorstore::OpenMode::delete_existing, + // tensorstore::ReadWriteMode::write).result() + // ); +} + +// Private Methods +void PyramidCompositor::create_xml() { + CreateXML( + this->_output_pyramid_name, + this->_ome_metadata_file, + this->_plate_image_shapes[0], + this->_image_dtype + ); +} + +void PyramidCompositor::create_zattr_file() { + WriteTSZattrFilePlateImage( + this->_output_pyramid_name, + this->_output_pyramid_name + "/data.zarr/0", + this->_plate_image_shapes + ); +} + +void PyramidCompositor::create_zgroup_file() { + WriteVivZgroupFiles(this->_output_pyramid_name); +} + +void PyramidCompositor::create_auxiliary_files() { + create_xml(); + create_zattr_file(); + create_zgroup_file(); +} + +void PyramidCompositor::write_zarr_chunk(int level, int channel, int y_int, int x_int) { + // Implementation of tile data writing logic + // Placeholder logic + + // std::tuple y_range = { + // y_int * CHUNK_SIZE, + // std::min((y_int+1) * CHUNK_SIZE, this->zarr_arrays_[level].size()[3]) + // }; + + // std::tuple x_range = { + // x_int * CHUNK_SIZE, + // min((x_int+1) * CHUNK_SIZE, this->zarr_arrays_[level].size()[4]) + // }; + + // int assembled_width = std::get<1>(x_range) - std::get<0>(x_range); + // int assembled_height = std::get<1>(y_range) - std::get<0>(y_range); + + // std::vector<> + + // Compute ranges in global coordinates + int y_start = y_int * CHUNK_SIZE; + int y_end = std::min((y_int + 1) * CHUNK_SIZE, std::get<1>(_zarr_arrays[level])[0]); + int x_start = x_int * CHUNK_SIZE; + int x_end = std::min((x_int + 1) * CHUNK_SIZE, std::get<1>(_zarr_arrays[level])[1]); + + int assembled_width = x_end - x_start; + int assembled_height = y_end - y_start; + + image_data assembled_image = GetAssembledImageVector(assembled_height*assembled_width); + + // Determine required input images + int unit_image_height = _unit_image_shapes[level].first; + int unit_image_width = _unit_image_shapes[level].second; + + int row_start_pos = y_start; + while (row_start_pos < y_end) { + int row = row_start_pos / unit_image_height; + int local_y_start = row_start_pos - y_start; + int tile_y_start = row_start_pos - row * unit_image_height; + int tile_y_dim = std::min((row + 1) * unit_image_height - row_start_pos, y_end - row_start_pos); + int tile_y_end = tile_y_start + tile_y_dim; + + int col_start_pos = x_start; + while (col_start_pos < x_end) { + int col = col_start_pos / unit_image_width; + int local_x_start = col_start_pos - x_start; + int tile_x_start = col_start_pos - col * unit_image_width; + int tile_x_dim = std::min((col + 1) * unit_image_width - col_start_pos, x_end - col_start_pos); + int tile_x_end = tile_x_start + tile_x_dim; + + // Fetch input file path from composition map + std::string input_file_name = _composition_map[{col, row, channel}]; + std::filesystem::path zarrArrayLoc = std::filesystem::path(input_file_name) / "data.zarr/0/" / std::to_string(level); + + auto read_data = std::get<0>(ReadImageData(zarrArrayLoc, Seq(tile_y_start, tile_y_end), Seq(tile_x_start, tile_x_end))).get(); + + for (int i = local_y_start; i < local_y_start + tile_y_end - tile_y_start; ++i) { + for (int j = local_x_start; i < local_x_start + tile_x_end - tile_x_start; ++j) { + assembled_image[i*(tile_x_end - tile_x_start) + j] = read_data[(i-local_y_start)*(tile_x_end - tile_x_start) + (j-local_x_start)]; + } + } + + // // Open the Zarr array for reading with TensorStore + // auto zarrSpec = tensorstore::Context::Default() + // | tensorstore::OpenMode::read + // | tensorstore::dtype_v + // | tensorstore::Shape({1, 1, 1, assembledHeight, assembledWidth}); + + // auto inputZarrArray = tensorstore::Open(zarrSpec).result(); + + // TENSORSTORE_CHECK_OK_AND_ASSIGN() + + // if (!inputZarrArray.ok()) { + // std::cerr << "Failed to open Zarr array: " << inputZarrArray.status() << "\n"; + // return; + // } + + // // Read data from input Zarr file into assembledImage + // tensorstore::Read(*inputZarrArray, + // tensorstore::Span(assembledImage.data(), assembledImage.size())).result(); + + col_start_pos += tile_x_end - tile_x_start; + } + row_start_pos += tile_y_end - tile_y_start; + } + + + WriteImageData( + _input_pyramids_loc, + std::get<0>(_zarr_arrays[level]).get(), + Seq(y_start,y_end), Seq(x_start, x_end), + Seq(0,0), Seq(channel, channel), Seq(0,0) + ); +} + +void PyramidCompositor::setComposition(const std::unordered_map, std::string>& comp_map) { + _composition_map = comp_map; + + for (const auto& coord : comp_map) { + std::string file_path = coord.second; + std::filesystem::path attr_file_loc = std::filesystem::path(file_path) / "data.zarr/0/.zattrs"; + + if (std::filesystem::exists(attr_file_loc)) { + std::ifstream attr_file(attr_file_loc); + json attrs; + attr_file >> attrs; + + const auto& multiscale_metadata = attrs["multiscales"][0]["datasets"]; + _pyramid_levels = multiscale_metadata.size(); + for (const auto& dic : multiscale_metadata) { + std::string res_key = dic["path"]; + std::filesystem::path zarr_array_loc = std::filesystem::path(file_path) / "data.zarr/0/" / res_key; + + auto read_spec = GetZarrSpecToRead(zarr_array_loc); + + tensorstore::TensorStore source; + + TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( + read_spec, + tensorstore::OpenMode::open, + tensorstore::ReadWriteMode::read).result()); + + auto image_shape = source.domain().shape(); + _image_dtype = source.dtype().name(); + _data_type_code = GetDataTypeCode(_image_dtype); + + _unit_image_shapes[std::stoi(res_key)] = { + image_shape[image_shape.size()-2], + image_shape[image_shape.size()-1] + }; + } + break; + } + } + + int num_rows = 0, num_cols = 0, num_channels = 0; + for (const auto& coord : comp_map) { + num_rows = std::max(num_rows, std::get<1>(coord.first)); + num_cols = std::max(num_cols, std::get<0>(coord.first)); + _num_channels = std::max(num_channels, std::get<2>(coord.first)); + } + num_cols += 1; + num_rows += 1; + num_channels += 1; + _num_channels = num_channels; + + for (const auto& [level, shape] : _unit_image_shapes) { + std::string path = _output_pyramid_name + "/data.zarr/0/" + std::to_string(level); + + auto zarr_array_data = ReadImageData( + path, + Seq(0, num_rows * shape.first), + Seq(0, num_cols * shape.second), + Seq(0,0), + Seq(0, num_channels), + Seq(0,0) + ); + + _zarr_arrays[level] = zarr_array_data; + } +} + +void PyramidCompositor::reset_composition() { + // Remove pyramid file and clear internal data structures + fs::remove_all(_output_pyramid_name); + _composition_map.clear(); + _plate_image_shapes.clear(); + _tile_cache.clear(); + _plate_image_shapes.clear(); + _zarr_arrays.clear(); +} + +image_data PyramidCompositor::GetAssembledImageVector(int size){ + switch (_image_dtype_code) + { + case (1): + return std::vector(size); + case (2): + return std::vector(size); + case (4): + return std::vector(size); + case (8): + return std::vector(size); + case (16): + return std::vector(size); + case (32): + return std::vector(size); + case (64): + return std::vector(size); + case (128): + return std::vector(size); + case (256): + return std::vector(size); + case (512): + return std::vector(size); + default: + throw std::runtime_error("Invalid data type."); + } +} + +void PyramidCompositor::WriteImageData( + const std::string& path, + image_data& image, + const Seq& rows, + const Seq& cols, + const std::optional& layers, + const std::optional& channels, + const std::optional& tsteps) { + + std::vector shape; + + tensorstore::TensorStore source; + + TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( + GetZarrSpecToRead(path), + tensorstore::OpenMode::open, + tensorstore::ReadWriteMode::write).result()); + + auto output_transform = tensorstore::IdentityTransform(source.domain()); + + if (_t_index.has_value() && tsteps.has_value()) { + output_transform = (std::move(output_transform) | tensorstore::Dims(_t_index.value()).ClosedInterval(tsteps.value().Start(), tsteps.value().Stop())).value(); + shape.emplace_back(tsteps.value().Stop() - tsteps.value().Start()+1); + } + + if (_c_index.has_value() && channels.has_value()) { + output_transform = (std::move(output_transform) | tensorstore::Dims(_c_index.value()).ClosedInterval(channels.value().Start(), channels.value().Stop())).value(); + shape.emplace_back(channels.value().Stop() - channels.value().Start()+1); + } + + if (_z_index.has_value() && layers.has_value()) { + output_transform = (std::move(output_transform) | tensorstore::Dims(_z_index.value()).ClosedInterval(layers.value().Start(), layers.value().Stop())).value(); + shape.emplace_back(layers.value().Stop() - layers.value().Start()+1); + } + + output_transform = (std::move(output_transform) | tensorstore::Dims(_y_index).ClosedInterval(rows.Start(), rows.Stop()) | + tensorstore::Dims(_x_index).ClosedInterval(cols.Start(), cols.Stop())).value(); + + shape.emplace_back(rows.Stop() - rows.Start()+1); + shape.emplace_back(cols.Stop() - cols.Start()+1); + + auto data_array = tensorstore::Array(image.data(), shape, tensorstore::c_order); + + // Write data array to TensorStore + auto write_result = tensorstore::Write(tensorstore::UnownedToShared(data_array), source | output_transform).result(); + + if (!write_result.ok()) { + std::cerr << "Error writing image: " << write_result.status() << std::endl; + } +} + +template +std::shared_ptr> PyramidCompositor::GetImageDataTemplated(const Seq& rows, const Seq& cols, const Seq& layers, const Seq& channels, const Seq& tsteps){ + + const auto data_height = rows.Stop() - rows.Start() + 1; + const auto data_width = cols.Stop() - cols.Start() + 1; + const auto data_depth = layers.Stop() - layers.Start() + 1; + const auto data_num_channels = channels.Stop() - channels.Start() + 1; + const auto data_tsteps = tsteps.Stop() - tsteps.Start() + 1; + + auto read_buffer = std::make_shared>(data_height*data_width*data_depth*data_num_channels*data_tsteps); + auto array = tensorstore::Array(read_buffer->data(), {data_tsteps, data_num_channels, data_depth, data_height, data_width}, tensorstore::c_order); + tensorstore::intTransform<> read_transform = tensorstore::IdentityTransform(source.domain()); + + int x_int=1, y_int=0; + if (_t_index.has_value()){ + read_transform = (std::move(read_transform) | tensorstore::Dims(_t_index.value()).ClosedInterval(tsteps.Start(), tsteps.Stop())).value(); + x_int++; + y_int++; + } + if (_c_index.has_value()){ + read_transform = (std::move(read_transform) | tensorstore::Dims(_c_index.value()).ClosedInterval(channels.Start(), channels.Stop())).value(); + x_int++; + y_int++; + } + if (_z_index.has_value()){ + read_transform = (std::move(read_transform) | tensorstore::Dims(_z_index.value()).ClosedInterval(layers.Start(), layers.Stop())).value(); + x_int++; + y_int++; + } + read_transform = (std::move(read_transform) | tensorstore::Dims(y_int).ClosedInterval(rows.Start(), rows.Stop()) | + tensorstore::Dims(x_int).ClosedInterval(cols.Start(), cols.Stop())).value(); + + + tensorstore::Read(source | read_transform, tensorstore::UnownedToShared(array)).value(); + + return read_buffer; +} + +std::tuple, std::vector> PyramidCompositor::ReadImageData( + const std::string& fname, + const Seq& rows, const Seq& cols, + const Seq& layers = Seq(0,0), + const Seq& channels = Seq(0,0), + const Seq& tsteps = Seq(0,0)) { + + auto read_spec = GetZarrSpecToRead(fname); + + tensorstore::TensorStore source; + + TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( + read_spec, + tensorstore::OpenMode::open, + tensorstore::ReadWriteMode::read).result()); + + auto image_shape = source.domain().shape(); + const auto read_chunk_shape = source.chunk_layout().value().read_chunk_shape(); + + // if (image_shape.size() == 5){ + // _image_height = image_shape[3]; + // _image_width = image_shape[4]; + // _image_depth = image_shape[2]; + // _num_channels = image_shape[1]; + // _num_tsteps = image_shape[0]; + // _t_index.emplace(0); + // _c_index.emplace(1); + // _z_index.emplace(2); + // } else { + // assert(image_shape.size() >= 2); + // std::tie(_t_index, _c_index, _z_index) = ParseMultiscaleMetadata(axes_list, image_shape.size()); + // _image_height = image_shape[image_shape.size()-2]; + // _image_width = image_shape[image_shape.size()-1]; + // if (_t_index.has_value()) { + // _num_tsteps = image_shape[_t_index.value()]; + // } else { + // _num_tsteps = 1; + // } + + // if (_c_index.has_value()) { + // _num_channels = image_shape[_c_index.value()]; + // } else { + // _num_channels = 1; + // } + + // if (_z_index.has_value()) { + // _image_depth = image_shape[_z_index.value()]; + // } else { + // _image_depth = 1; + // } + // } + + _image_dtype = source.dtype().name(); + _data_type_code = GetDataTypeCode(_image_dtype); + + switch (_data_type_code) + { + case (1): + return std::make_tuple( + std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), + std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); + case (2): + return std::make_tuple( + std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), + std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); + case (4): + return std::make_tuple( + std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), + std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); + case (8): + return std::make_tuple( + std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), + std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); + case (16): + return std::make_tuple( + std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), + std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); + case (32): + return std::make_tuple( + std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), + std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); + case (64): + return std::make_tuple( + std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), + std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); + case (128): + return std::make_tuple( + std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), + std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); + case (256): + return std::make_tuple( + std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), + std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); + case (512): + return std::make_tuple( + std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), + std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); + default: + break; + } +} +} // ns argolid \ No newline at end of file diff --git a/src/cpp/utilities/utilities.cpp b/src/cpp/utilities/utilities.cpp index 3f52fee..e51187a 100644 --- a/src/cpp/utilities/utilities.cpp +++ b/src/cpp/utilities/utilities.cpp @@ -441,4 +441,29 @@ std::optional> GetTiffDims (const std:: } +std::tuple, std::optional, std::optional>ParseMultiscaleMetadata(const std::string& axes_list, int len){ + + std::optional t_index{std::nullopt}, c_index{std::nullopt}, z_index{std::nullopt}; + + assert(axes_list.length() <= 5); + + if (axes_list.length() == len){ + // no speculation + for (int i=0; i>& plate_image_shape); +std::tuple, std::optional, std::optional>ParseMultiscaleMetadata(const std::string& axes_list, int len); + } // ns argolid \ No newline at end of file From c9173a06234a18aec1be90c7e7976010c03477cf Mon Sep 17 00:00:00 2001 From: JesseMckinzie Date: Thu, 16 Jan 2025 08:31:19 -0700 Subject: [PATCH 03/24] Update cpp PyramidCompositor --- src/cpp/core/pyramid_compositor.cpp | 435 +++++++++++------- src/cpp/core/pyramid_compositor.h | 98 +++- .../interface/pyramid_python_interface.cpp | 16 +- src/python/argolid/__init__.py | 2 +- src/python/argolid/pyramid_compositor.py | 2 +- 5 files changed, 351 insertions(+), 202 deletions(-) diff --git a/src/cpp/core/pyramid_compositor.cpp b/src/cpp/core/pyramid_compositor.cpp index 08ae8a1..3969e7a 100644 --- a/src/cpp/core/pyramid_compositor.cpp +++ b/src/cpp/core/pyramid_compositor.cpp @@ -36,6 +36,7 @@ PyramidCompositor::PyramidCompositor(const std::string& input_pyramids_loc, cons tensorstore::OpenMode::open, tensorstore::ReadWriteMode::read).result()); + tensorstore::DataType _image_ts_dtype = source.dtype(); _image_dtype = source.dtype().name(); _image_dtype_code = GetDataTypeCode(_image_dtype); @@ -98,6 +99,11 @@ void PyramidCompositor::create_auxiliary_files() { } void PyramidCompositor::write_zarr_chunk(int level, int channel, int y_int, int x_int) { + this->write_zarr_chunk(level, channel, y_int, x_int); +} + +template +void PyramidCompositor::_write_zarr_chunk(int level, int channel, int y_int, int x_int) { // Implementation of tile data writing logic // Placeholder logic @@ -125,7 +131,14 @@ void PyramidCompositor::write_zarr_chunk(int level, int channel, int y_int, int int assembled_width = x_end - x_start; int assembled_height = y_end - y_start; - image_data assembled_image = GetAssembledImageVector(assembled_height*assembled_width); + std::vector assembled_image_vec(assembled_width*assembled_height); + + // auto assembled_image = tensorstore::AllocateArray({ + // assembled_width, + // assembled_height + // }, tensorstore::c_order, + // tensorstore::value_init, + // _image_ts_dtype); // Determine required input images int unit_image_height = _unit_image_shapes[level].first; @@ -151,21 +164,100 @@ void PyramidCompositor::write_zarr_chunk(int level, int channel, int y_int, int std::string input_file_name = _composition_map[{col, row, channel}]; std::filesystem::path zarrArrayLoc = std::filesystem::path(input_file_name) / "data.zarr/0/" / std::to_string(level); - auto read_data = *std::get<0>(ReadImageData(zarrArrayLoc, Seq(tile_y_start, tile_y_end), Seq(tile_x_start, tile_x_end))).get(); + // Read inage + //auto read_data = ReadImageData(zarrArrayLoc, Seq(tile_y_start, tile_y_end), Seq(tile_x_start, tile_x_end)); + + auto read_spec = GetZarrSpecToRead(zarrArrayLoc.u8string()); + + tensorstore::TensorStore source; + + TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( + read_spec, + tensorstore::OpenMode::open, + tensorstore::ReadWriteMode::read).result()); + + + std::vector read_buffer((tile_x_end-tile_x_start)*(tile_y_end-tile_y_start)); + auto array = tensorstore::Array(read_buffer.data(), {tile_y_end-tile_y_start, tile_x_end-tile_x_start}, tensorstore::c_order); + + tensorstore::IndexTransform<> read_transform = tensorstore::IdentityTransform(source.domain()); + + Seq rows = Seq(y_start, y_end); + Seq cols = Seq(x_start, x_end); + Seq tsteps = Seq(0, 0); + Seq channels = Seq(channel, channel); + Seq layers = Seq(0, 0); + + int x_int=1, y_int=0; + if (_t_index.has_value()){ + read_transform = (std::move(read_transform) | tensorstore::Dims(_t_index.value()).ClosedInterval(tsteps.Start(), tsteps.Stop())).value(); + x_int++; + y_int++; + } + if (_c_index.has_value()){ + read_transform = (std::move(read_transform) | tensorstore::Dims(_c_index.value()).ClosedInterval(channels.Start(), channels.Stop())).value(); + x_int++; + y_int++; + } + if (_z_index.has_value()){ + read_transform = (std::move(read_transform) | tensorstore::Dims(_z_index.value()).ClosedInterval(layers.Start(), layers.Stop())).value(); + x_int++; + y_int++; + } + read_transform = (std::move(read_transform) | tensorstore::Dims(y_int).ClosedInterval(rows.Start(), rows.Stop()) | + tensorstore::Dims(x_int).ClosedInterval(cols.Start(), cols.Stop())).value(); + + + tensorstore::Read(source | read_transform, tensorstore::UnownedToShared(array)).value(); + + for (int i = local_y_start; i < local_y_start + tile_y_end - tile_y_start; ++i) { + for (int j = local_x_start; j < local_x_start + tile_x_end - tile_x_start; ++j) { + assembled_image_vec[i * assembled_width + j] = read_buffer[(i - local_y_start) * (tile_x_end - tile_x_start) + (j - local_x_start)]; + } + } + + // // auto read_spec = GetZarrSpecToRead(zarrArrayLoc.u8string()); + + // // tensorstore::TensorStore source; + + // // TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( + // // read_spec, + // // tensorstore::OpenMode::open, + // // tensorstore::ReadWriteMode::read).result()); + // // auto image_shape = source.domain().shape(); + // // const auto read_chunk_shape = source.chunk_layout().value().read_chunk_shape(); + + // // auto image_width = image_shape[4]; + // // auto image_height = image_shape[3]; + + // // auto array = tensorstore::AllocateArray({ + // // image_height, + // // image_width + // // }, tensorstore::c_order, + // // tensorstore::value_init, source.dtype()); + + // // // initiate a read + // // tensorstore::Read(source | + // // tensorstore::Dims(3).ClosedInterval(0, image_height - 1) | + // // tensorstore::Dims(4).ClosedInterval(0, image_width - 1), + // // array).value(); + + //auto read_data = *std::get<0>(ReadImageData(zarrArrayLoc, Seq(tile_y_start, tile_y_end), Seq(tile_x_start, tile_x_end))).get(); + // Use std::visit to handle the different types in the variant - std::visit([&](auto& assembled_vec) { - using T = typename std::decay_t::value_type; + // std::visit([&](auto& assembled_vec) { + // using T = typename std::decay_t::value_type; - for (int i = local_y_start; i < local_y_start + tile_y_end - tile_y_start; ++i) { - for (int j = local_x_start; j < local_x_start + tile_x_end - tile_x_start; ++j) { - // Use std::visit to extract the value from read_data - std::visit([&](auto& value) { - assembled_vec[i * assembled_width + j] = static_cast(value[(i - local_y_start) * (tile_x_end - tile_x_start) + (j - local_x_start)]); - }, read_data); - } - } - }, assembled_image); + // for (int i = local_y_start; i < local_y_start + tile_y_end - tile_y_start; ++i) { + // for (int j = local_x_start; j < local_x_start + tile_x_end - tile_x_start; ++j) { + // // Use std::visit to extract the value from read_data + // std::visit([&](auto& value) { + // assembled_vec[i * assembled_width + j] = static_cast(value[(i - local_y_start) * (tile_x_end - tile_x_start) + (j - local_x_start)]); + // }, read_data); + // } + // } + // }, assembled_image); // for (int i = local_y_start; i < local_y_start + tile_y_end - tile_y_start; ++i) { @@ -198,7 +290,17 @@ void PyramidCompositor::write_zarr_chunk(int level, int channel, int y_int, int row_start_pos += tile_y_end - tile_y_start; } - auto data = *(std::get<0>(_zarr_arrays[level]).get()); + WriteImageData( + _input_pyramids_loc, + assembled_image_vec, // Access the data from the std::shared_ptr + Seq(y_start, y_end), + Seq(x_start, x_end), + Seq(0, 0), + Seq(channel, channel), + Seq(0, 0) + ); + + //auto data = *(std::get<0>(_zarr_arrays[level]).get()); // WriteImageData( // _input_pyramids_loc, // data, @@ -209,20 +311,20 @@ void PyramidCompositor::write_zarr_chunk(int level, int channel, int y_int, int // Assuming _zarr_arrays[level] is of type std::tuple>, std::vector> //auto& [variant_ptr, some_vector] = _zarr_arrays[level]; // Structured binding (C++17) - std::visit([&](const auto& array) { - using T = std::decay_t; - - // Use external variables and call WriteImageData - WriteImageData( - _input_pyramids_loc, - array, // Access the data from the std::shared_ptr - Seq(y_start, y_end), - Seq(x_start, x_end), - Seq(0, 0), - Seq(channel, channel), - Seq(0, 0) - ); - }, data); // Dereference the shared_ptr to get the variant + // std::visit([&](const auto& array) { + // using T = std::decay_t; + + // // Use external variables and call WriteImageData + // WriteImageData( + // _input_pyramids_loc, + // array, // Access the data from the std::shared_ptr + // Seq(y_start, y_end), + // Seq(x_start, x_end), + // Seq(0, 0), + // Seq(channel, channel), + // Seq(0, 0) + // ); + //}, data); // Dereference the shared_ptr to get the variant } void PyramidCompositor::setComposition(const std::unordered_map, std::string, TupleHash>& comp_map) { @@ -276,20 +378,20 @@ void PyramidCompositor::setComposition(const std::unordered_map(zarr_array_data); + // } } void PyramidCompositor::reset_composition() { @@ -302,35 +404,6 @@ void PyramidCompositor::reset_composition() { _zarr_arrays.clear(); } -image_data PyramidCompositor::GetAssembledImageVector(int size){ - switch (_image_dtype_code) - { - case (1): - return std::vector(size); - case (2): - return std::vector(size); - case (4): - return std::vector(size); - case (8): - return std::vector(size); - case (16): - return std::vector(size); - case (32): - return std::vector(size); - case (64): - return std::vector(size); - case (128): - return std::vector(size); - case (256): - return std::vector(size); - case (512): - return std::vector(size); - default: - throw std::runtime_error("Invalid data type."); - } -} - - template void PyramidCompositor::WriteImageData( const std::string& path, @@ -345,6 +418,8 @@ void PyramidCompositor::WriteImageData( tensorstore::TensorStore source; + auto result_array = tensorstore::Array(image->data(), {rows.Stop()-rows.Start(), cols.Stop()-cols.Start()}, tensorstore::c_order); + TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( GetZarrSpecToRead(path), tensorstore::OpenMode::open, @@ -373,33 +448,49 @@ void PyramidCompositor::WriteImageData( shape.emplace_back(rows.Stop() - rows.Start()+1); shape.emplace_back(cols.Stop() - cols.Start()+1); - auto data_array = tensorstore::Array(image.data(), shape, tensorstore::c_order); // Write data array to TensorStore - auto write_result = tensorstore::Write(tensorstore::UnownedToShared(data_array), source | output_transform).result(); + auto write_result = tensorstore::Write(tensorstore::UnownedToShared(result_array), source | output_transform).result(); if (!write_result.ok()) { std::cerr << "Error writing image: " << write_result.status() << std::endl; } } -template -std::shared_ptr> PyramidCompositor::GetImageDataTemplated( - tensorstore::TensorStore& source, + +auto PyramidCompositor::ReadImageData( + const std::string& fname, const Seq& rows, const Seq& cols, const Seq& layers, const Seq& channels, const Seq& tsteps) { + auto read_spec = GetZarrSpecToRead(fname); + + tensorstore::TensorStore source; + + TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( + read_spec, + tensorstore::OpenMode::open, + tensorstore::ReadWriteMode::read).result()); + + auto image_shape = source.domain().shape(); + const auto read_chunk_shape = source.chunk_layout().value().read_chunk_shape(); + const auto data_height = rows.Stop() - rows.Start() + 1; const auto data_width = cols.Stop() - cols.Start() + 1; const auto data_depth = layers.Stop() - layers.Start() + 1; const auto data_num_channels = channels.Stop() - channels.Start() + 1; const auto data_tsteps = tsteps.Stop() - tsteps.Start() + 1; - auto read_buffer = std::make_shared>(data_height*data_width*data_depth*data_num_channels*data_tsteps); - auto array = tensorstore::Array(read_buffer->data(), {data_tsteps, data_num_channels, data_depth, data_height, data_width}, tensorstore::c_order); + auto array = tensorstore::AllocateArray( + {data_tsteps, data_num_channels, data_depth, data_height, data_width}, + tensorstore::c_order, + tensorstore::value_init, + source.dtype() + ); + tensorstore::IndexTransform<> read_transform = tensorstore::IdentityTransform(source.domain()); int x_int=1, y_int=0; @@ -424,108 +515,108 @@ std::shared_ptr> PyramidCompositor::GetImageDataTemplated( tensorstore::Read(source | read_transform, tensorstore::UnownedToShared(array)).value(); - return read_buffer; + return array; } -std::tuple, std::vector> PyramidCompositor::ReadImageData( - const std::string& fname, - const Seq& rows, const Seq& cols, - const Seq& layers, - const Seq& channels, - const Seq& tsteps) { +// std::tuple, std::vector> PyramidCompositor::ReadImageData( +// const std::string& fname, +// const Seq& rows, const Seq& cols, +// const Seq& layers, +// const Seq& channels, +// const Seq& tsteps) { - auto read_spec = GetZarrSpecToRead(fname); +// auto read_spec = GetZarrSpecToRead(fname); - tensorstore::TensorStore source; +// tensorstore::TensorStore source; - TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( - read_spec, - tensorstore::OpenMode::open, - tensorstore::ReadWriteMode::read).result()); +// TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( +// read_spec, +// tensorstore::OpenMode::open, +// tensorstore::ReadWriteMode::read).result()); - auto image_shape = source.domain().shape(); - const auto read_chunk_shape = source.chunk_layout().value().read_chunk_shape(); +// auto image_shape = source.domain().shape(); +// const auto read_chunk_shape = source.chunk_layout().value().read_chunk_shape(); - // if (image_shape.size() == 5){ - // _image_height = image_shape[3]; - // _image_width = image_shape[4]; - // _image_depth = image_shape[2]; - // _num_channels = image_shape[1]; - // _num_tsteps = image_shape[0]; - // _t_index.emplace(0); - // _c_index.emplace(1); - // _z_index.emplace(2); - // } else { - // assert(image_shape.size() >= 2); - // std::tie(_t_index, _c_index, _z_index) = ParseMultiscaleMetadata(axes_list, image_shape.size()); - // _image_height = image_shape[image_shape.size()-2]; - // _image_width = image_shape[image_shape.size()-1]; - // if (_t_index.has_value()) { - // _num_tsteps = image_shape[_t_index.value()]; - // } else { - // _num_tsteps = 1; - // } - - // if (_c_index.has_value()) { - // _num_channels = image_shape[_c_index.value()]; - // } else { - // _num_channels = 1; - // } - - // if (_z_index.has_value()) { - // _image_depth = image_shape[_z_index.value()]; - // } else { - // _image_depth = 1; - // } - // } +// // if (image_shape.size() == 5){ +// // _image_height = image_shape[3]; +// // _image_width = image_shape[4]; +// // _image_depth = image_shape[2]; +// // _num_channels = image_shape[1]; +// // _num_tsteps = image_shape[0]; +// // _t_index.emplace(0); +// // _c_index.emplace(1); +// // _z_index.emplace(2); +// // } else { +// // assert(image_shape.size() >= 2); +// // std::tie(_t_index, _c_index, _z_index) = ParseMultiscaleMetadata(axes_list, image_shape.size()); +// // _image_height = image_shape[image_shape.size()-2]; +// // _image_width = image_shape[image_shape.size()-1]; +// // if (_t_index.has_value()) { +// // _num_tsteps = image_shape[_t_index.value()]; +// // } else { +// // _num_tsteps = 1; +// // } + +// // if (_c_index.has_value()) { +// // _num_channels = image_shape[_c_index.value()]; +// // } else { +// // _num_channels = 1; +// // } + +// // if (_z_index.has_value()) { +// // _image_depth = image_shape[_z_index.value()]; +// // } else { +// // _image_depth = 1; +// // } +// // } - _image_dtype = source.dtype().name(); - _image_dtype_code = GetDataTypeCode(_image_dtype); +// _image_dtype = source.dtype().name(); +// _image_dtype_code = GetDataTypeCode(_image_dtype); - switch (_data_type_code) - { - case (1): - return std::make_tuple( - std::make_shared(std::move(*(GetImageDataTemplated(source, rows, cols, layers, channels, tsteps)))), - std::vector{image_shape[image_shape.size()-2], image_shape[image_shape.size()-1]}); - case (2): - return std::make_tuple( - std::make_shared(std::move(*(GetImageDataTemplated(source, rows, cols, layers, channels, tsteps)))), - std::vector{image_shape[image_shape.size()-2], image_shape[image_shape.size()-1]}); - case (4): - return std::make_tuple( - std::make_shared(std::move(*(GetImageDataTemplated(source, rows, cols, layers, channels, tsteps)))), - std::vector{image_shape[image_shape.size()-2], image_shape[image_shape.size()-1]}); - case (8): - return std::make_tuple( - std::make_shared(std::move(*(GetImageDataTemplated(source, rows, cols, layers, channels, tsteps)))), - std::vector{image_shape[image_shape.size()-2], image_shape[image_shape.size()-1]}); - case (16): - return std::make_tuple( - std::make_shared(std::move(*(GetImageDataTemplated(source, rows, cols, layers, channels, tsteps)))), - std::vector{image_shape[image_shape.size()-2], image_shape[image_shape.size()-1]}); - case (32): - return std::make_tuple( - std::make_shared(std::move(*(GetImageDataTemplated(source, rows, cols, layers, channels, tsteps)))), - std::vector{image_shape[image_shape.size()-2], image_shape[image_shape.size()-1]}); - case (64): - return std::make_tuple( - std::make_shared(std::move(*(GetImageDataTemplated(source, rows, cols, layers, channels, tsteps)))), - std::vector{image_shape[image_shape.size()-2], image_shape[image_shape.size()-1]}); - case (128): - return std::make_tuple( - std::make_shared(std::move(*(GetImageDataTemplated(source, rows, cols, layers, channels, tsteps)))), - std::vector{image_shape[image_shape.size()-2], image_shape[image_shape.size()-1]}); - case (256): - return std::make_tuple( - std::make_shared(std::move(*(GetImageDataTemplated(source, rows, cols, layers, channels, tsteps)))), - std::vector{image_shape[image_shape.size()-2], image_shape[image_shape.size()-1]}); - case (512): - return std::make_tuple( - std::make_shared(std::move(*(GetImageDataTemplated(source, rows, cols, layers, channels, tsteps)))), - std::vector{image_shape[image_shape.size()-2], image_shape[image_shape.size()-1]}); - default: - break; - } -} +// switch (_data_type_code) +// { +// case (1): +// return std::make_tuple( +// std::make_shared(std::move(*(GetImageDataTemplated(source, rows, cols, layers, channels, tsteps)))), +// std::vector{image_shape[image_shape.size()-2], image_shape[image_shape.size()-1]}); +// case (2): +// return std::make_tuple( +// std::make_shared(std::move(*(GetImageDataTemplated(source, rows, cols, layers, channels, tsteps)))), +// std::vector{image_shape[image_shape.size()-2], image_shape[image_shape.size()-1]}); +// case (4): +// return std::make_tuple( +// std::make_shared(std::move(*(GetImageDataTemplated(source, rows, cols, layers, channels, tsteps)))), +// std::vector{image_shape[image_shape.size()-2], image_shape[image_shape.size()-1]}); +// case (8): +// return std::make_tuple( +// std::make_shared(std::move(*(GetImageDataTemplated(source, rows, cols, layers, channels, tsteps)))), +// std::vector{image_shape[image_shape.size()-2], image_shape[image_shape.size()-1]}); +// case (16): +// return std::make_tuple( +// std::make_shared(std::move(*(GetImageDataTemplated(source, rows, cols, layers, channels, tsteps)))), +// std::vector{image_shape[image_shape.size()-2], image_shape[image_shape.size()-1]}); +// case (32): +// return std::make_tuple( +// std::make_shared(std::move(*(GetImageDataTemplated(source, rows, cols, layers, channels, tsteps)))), +// std::vector{image_shape[image_shape.size()-2], image_shape[image_shape.size()-1]}); +// case (64): +// return std::make_tuple( +// std::make_shared(std::move(*(GetImageDataTemplated(source, rows, cols, layers, channels, tsteps)))), +// std::vector{image_shape[image_shape.size()-2], image_shape[image_shape.size()-1]}); +// case (128): +// return std::make_tuple( +// std::make_shared(std::move(*(GetImageDataTemplated(source, rows, cols, layers, channels, tsteps)))), +// std::vector{image_shape[image_shape.size()-2], image_shape[image_shape.size()-1]}); +// case (256): +// return std::make_tuple( +// std::make_shared(std::move(*(GetImageDataTemplated(source, rows, cols, layers, channels, tsteps)))), +// std::vector{image_shape[image_shape.size()-2], image_shape[image_shape.size()-1]}); +// case (512): +// return std::make_tuple( +// std::make_shared(std::move(*(GetImageDataTemplated(source, rows, cols, layers, channels, tsteps)))), +// std::vector{image_shape[image_shape.size()-2], image_shape[image_shape.size()-1]}); +// default: +// break; +// } +// } } // ns argolid \ No newline at end of file diff --git a/src/cpp/core/pyramid_compositor.h b/src/cpp/core/pyramid_compositor.h index 860e748..8d5029b 100644 --- a/src/cpp/core/pyramid_compositor.h +++ b/src/cpp/core/pyramid_compositor.h @@ -23,17 +23,60 @@ struct TupleHash { } }; -using image_data = std::variant, - std::vector, - std::vector, - std::vector, - std::vector, - std::vector, - std::vector, - std::vector, - std::vector, - std::vector>; - +// using image_data = std::variant, +// std::vector, +// std::vector, +// std::vector, +// std::vector, +// std::vector, +// std::vector, +// std::vector, +// std::vector, +// std::vector>; + +using image_data = std::variant, + tensorstore::Array, + tensorstore::Array, + tensorstore::Array, + tensorstore::Array, + tensorstore::Array, + tensorstore::Array, + tensorstore::Array, + tensorstore::Array, + tensorstore::Array, + + tensorstore::Array, + tensorstore::Array, + tensorstore::Array, + tensorstore::Array, + tensorstore::Array, + tensorstore::Array, + tensorstore::Array, + tensorstore::Array, + tensorstore::Array, + tensorstore::Array, + + tensorstore::Array, + tensorstore::Array, + tensorstore::Array, + tensorstore::Array, + tensorstore::Array, + tensorstore::Array, + tensorstore::Array, + tensorstore::Array, + tensorstore::Array, + tensorstore::Array, + + tensorstore::Array, + tensorstore::Array, + tensorstore::Array, + tensorstore::Array, + tensorstore::Array, + tensorstore::Array, + tensorstore::Array, + tensorstore::Array, + tensorstore::Array, + tensorstore::Array>; class Seq { private: @@ -56,7 +99,12 @@ class PyramidCompositor { void create_zattr_file(); void create_zgroup_file(); void create_auxiliary_files(); + + template + void _write_zarr_chunk(int level, int channel, int y_index, int x_index); + void write_zarr_chunk(int level, int channel, int y_index, int x_index); + private: template @@ -69,28 +117,27 @@ class PyramidCompositor { const std::optional& channels, const std::optional& tsteps); - template - std::shared_ptr> GetImageDataTemplated( - tensorstore::TensorStore& source, + auto ReadImageData( + const std::string& fname, const Seq& rows, const Seq& cols, - const Seq& layers, - const Seq& channels, - const Seq& tsteps); - - std::tuple, std::vector> ReadImageData( - const std::string& fname, - const Seq& rows, const Seq& cols, const Seq& layers = Seq(0,0), const Seq& channels = Seq(0,0), - const Seq& tsteps = Seq(0,0) - ); + const Seq& tsteps = Seq(0,0)); + + // std::tuple, std::vector> ReadImageData( + // const std::string& fname, + // const Seq& rows, const Seq& cols, + // const Seq& layers = Seq(0,0), + // const Seq& channels = Seq(0,0), + // const Seq& tsteps = Seq(0,0) + // ); image_data GetAssembledImageVector(int size); void setComposition(const std::unordered_map, std::string, TupleHash>& comp_map); - // members for pyramid composition + // members for pyramid composition std::string _input_pyramids_loc; std::string _output_pyramid_name; std::string _ome_metadata_file; @@ -98,6 +145,8 @@ class PyramidCompositor { // std::unordered_map> _zarr_arrays; std::unordered_map, std::vector>> _zarr_arrays; // tuple at 0 is the image and at 1 is the size + //std::unordered_map + std::unordered_map> _unit_image_shapes; std::unordered_map, std::string, TupleHash> _composition_map; @@ -105,6 +154,7 @@ class PyramidCompositor { int _pyramid_levels; int _num_channels; + tensorstore::DataType _image_ts_dtype; std::string _image_dtype; std::uint16_t _image_dtype_code; diff --git a/src/cpp/interface/pyramid_python_interface.cpp b/src/cpp/interface/pyramid_python_interface.cpp index 8a4031b..9f6d731 100644 --- a/src/cpp/interface/pyramid_python_interface.cpp +++ b/src/cpp/interface/pyramid_python_interface.cpp @@ -20,14 +20,22 @@ PYBIND11_MODULE(libargolid, m) { py::class_(m, "PyramidCompositorCPP") \ - .def(py::init) \ - .def("set_well_map", &argolid::PyramidCompositor::set_well_map) \ + .def(py::init()) \ .def("reset_composition", &argolid::PyramidCompositor::reset_composition) \ - .def("get_tile_data", &argolid::PyramidCompositor::get_tile_data) \ .def("create_xml", &argolid::PyramidCompositor::create_xml) \ .def("create_zattr_file", &argolid::PyramidCompositor::create_zattr_file) \ .def("create_auxiliary_files", &argolid::PyramidCompositor::create_auxiliary_files) \ - .def("write_zarr_chunk", &argolid::PyramidCompositor::write_zarr_chunk) ; + .def("write_zarr_chunk", &argolid::PyramidCompositor::write_zarr_chunk); + // .def("write_zarr_chunk", &argolid::PyramidCompositor::write_zarr_chunk) \ + // .def("write_zarr_chunk", &argolid::PyramidCompositor::write_zarr_chunk) \ + // .def("write_zarr_chunk", &argolid::PyramidCompositor::write_zarr_chunk) \ + // .def("write_zarr_chunk", &argolid::PyramidCompositor::write_zarr_chunk) \ + // .def("write_zarr_chunk", &argolid::PyramidCompositor::write_zarr_chunk) \ + // .def("write_zarr_chunk", &argolid::PyramidCompositor::write_zarr_chunk) \ + // .def("write_zarr_chunk", &argolid::PyramidCompositor::write_zarr_chunk) \ + // .def("write_zarr_chunk", &argolid::PyramidCompositor::write_zarr_chunk) \ + // .def("write_zarr_chunk", &argolid::PyramidCompositor::write_zarr_chunk) \ + // .def("write_zarr_chunk", &argolid::PyramidCompositor::write_zarr_chunk) ; py::enum_(m, "VisType") .value("NG_Zarr", argolid::VisType::NG_Zarr) diff --git a/src/python/argolid/__init__.py b/src/python/argolid/__init__.py index bb5933b..69b8d46 100644 --- a/src/python/argolid/__init__.py +++ b/src/python/argolid/__init__.py @@ -1,5 +1,5 @@ from .pyramid_generator import PyramidGenerartor, PyramidView, PlateVisualizationMetadata, Downsample -from .pyramid_compositor import PyramidCompositor +from .pyramid_compositor import PyramidCompositorPY, PyramidCompositor2 from .volume_generator import VolumeGenerator, PyramidGenerator3D from . import _version diff --git a/src/python/argolid/pyramid_compositor.py b/src/python/argolid/pyramid_compositor.py index 37b31e7..0189277 100644 --- a/src/python/argolid/pyramid_compositor.py +++ b/src/python/argolid/pyramid_compositor.py @@ -321,7 +321,7 @@ def set_composition(self, composition_map: dict) -> None: self._chunk_cache = set() for l in self._unit_image_shapes: level = int(l) - self._plate_image_shapes[level] = ( + self._plate_image_shapes[level] = (1, 1, num_channels, 1, From b6f40747dcc07e4a628592a8e4e14d431302c4b8 Mon Sep 17 00:00:00 2001 From: JesseMckinzie Date: Tue, 21 Jan 2025 13:17:17 -0700 Subject: [PATCH 04/24] Fix issues in cpp pyramid compositor --- src/cpp/core/pyramid_compositor.cpp | 296 ++++++++++-------- src/cpp/core/pyramid_compositor.h | 13 +- .../interface/pyramid_python_interface.cpp | 3 +- src/cpp/utilities/utilities.cpp | 21 +- src/cpp/utilities/utilities.h | 4 +- src/python/argolid/pyramid_compositor.py | 27 +- 6 files changed, 221 insertions(+), 143 deletions(-) diff --git a/src/cpp/core/pyramid_compositor.cpp b/src/cpp/core/pyramid_compositor.cpp index 3969e7a..5fab267 100644 --- a/src/cpp/core/pyramid_compositor.cpp +++ b/src/cpp/core/pyramid_compositor.cpp @@ -23,41 +23,41 @@ namespace fs = std::filesystem; namespace argolid { PyramidCompositor::PyramidCompositor(const std::string& input_pyramids_loc, const std::string& out_dir, const std::string& output_pyramid_name): _input_pyramids_loc(input_pyramids_loc), - _output_pyramid_name(out_dir + "/" + output_pyramid_name), + _output_pyramid_name(std::filesystem::path(out_dir) / output_pyramid_name), _ome_metadata_file(out_dir + "/" + output_pyramid_name + "/METADATA.ome.xml"), _pyramid_levels(-1), _num_channels(-1) { + std::cerr << "Constructor" << std::endl; + // tensorstore::TensorStore source; - tensorstore::TensorStore source; - - TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( - GetZarrSpecToRead(input_pyramids_loc), - tensorstore::OpenMode::open, - tensorstore::ReadWriteMode::read).result()); + // TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( + // GetZarrSpecToRead(input_pyramids_loc), + // tensorstore::OpenMode::open, + // tensorstore::ReadWriteMode::read).result()); - tensorstore::DataType _image_ts_dtype = source.dtype(); - _image_dtype = source.dtype().name(); - _image_dtype_code = GetDataTypeCode(_image_dtype); + // tensorstore::DataType _image_ts_dtype = source.dtype(); + // _image_dtype = source.dtype().name(); + // _image_dtype_code = GetDataTypeCode(_image_dtype); - auto read_spec = GetZarrSpecToRead(input_pyramids_loc); + // auto read_spec = GetZarrSpecToRead(input_pyramids_loc); - TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( - read_spec, - tensorstore::OpenMode::open, - tensorstore::ReadWriteMode::read).result()); + // TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( + // read_spec, + // tensorstore::OpenMode::open, + // tensorstore::ReadWriteMode::read).result()); - auto image_shape = source.domain().shape(); - const auto read_chunk_shape = source.chunk_layout().value().read_chunk_shape(); - - if (image_shape.size() == 5){ - _t_index.emplace(0); - _c_index.emplace(1); - _z_index.emplace(2); - } else { - assert(image_shape.size() >= 2); - std::tie(_t_index, _c_index, _z_index) = ParseMultiscaleMetadata("XYCTZ", image_shape.size()); - } + // auto image_shape = source.domain().shape(); + // const auto read_chunk_shape = source.chunk_layout().value().read_chunk_shape(); + + // if (image_shape.size() == 5){ + // _t_index.emplace(0); + // _c_index.emplace(1); + // _z_index.emplace(2); + // } else { + // assert(image_shape.size() >= 2); + // std::tie(_t_index, _c_index, _z_index) = ParseMultiscaleMetadata("XYCTZ", image_shape.size()); + // } // _data_type = source.dtype().name(); // _data_type_code = GetDataTypeCode(_data_type); @@ -68,83 +68,129 @@ PyramidCompositor::PyramidCompositor(const std::string& input_pyramids_loc, cons // tensorstore::OpenMode::delete_existing, // tensorstore::ReadWriteMode::write).result() // ); + + std::cerr << "End Constructor" << std::endl; } // Private Methods void PyramidCompositor::create_xml() { + std::cerr << "XML" << std::endl; CreateXML( this->_output_pyramid_name, this->_ome_metadata_file, this->_plate_image_shapes[0], this->_image_dtype ); + std::cerr << "End XML" << std::endl; } void PyramidCompositor::create_zattr_file() { + std::cerr << "zattr" << std::endl; WriteTSZattrFilePlateImage( this->_output_pyramid_name, - this->_output_pyramid_name + "/data.zarr/0", + this->_output_pyramid_name.u8string() + "/data.zarr/0", this->_plate_image_shapes ); + + std::cerr << "End zattr" << std::endl; } void PyramidCompositor::create_zgroup_file() { + std::cerr << "zgroup" << std::endl; WriteVivZgroupFiles(this->_output_pyramid_name); + std::cerr << "End zgroup" << std::endl; } void PyramidCompositor::create_auxiliary_files() { + std::cerr << "aux" << std::endl; create_xml(); create_zattr_file(); create_zgroup_file(); + std::cerr << "End aux" << std::endl; } void PyramidCompositor::write_zarr_chunk(int level, int channel, int y_int, int x_int) { - this->write_zarr_chunk(level, channel, y_int, x_int); -} + std::cerr << "zarr wrapper" << std::endl; -template -void PyramidCompositor::_write_zarr_chunk(int level, int channel, int y_int, int x_int) { - // Implementation of tile data writing logic - // Placeholder logic + // get datatype by opening + std::string input_file_name = _composition_map[{x_int, y_int, channel}]; + std::filesystem::path zarrArrayLoc = _output_pyramid_name / input_file_name / std::filesystem::path("data.zarr/0/") / std::to_string(level); - // std::tuple y_range = { - // y_int * CHUNK_SIZE, - // std::min((y_int+1) * CHUNK_SIZE, this->zarr_arrays_[level].size()[3]) - // }; + auto read_spec = GetZarrSpecToRead(zarrArrayLoc.u8string()); - // std::tuple x_range = { - // x_int * CHUNK_SIZE, - // min((x_int+1) * CHUNK_SIZE, this->zarr_arrays_[level].size()[4]) - // }; + std::cerr << "Read path 121: " << zarrArrayLoc.u8string() << std::endl; - // int assembled_width = std::get<1>(x_range) - std::get<0>(x_range); - // int assembled_height = std::get<1>(y_range) - std::get<0>(y_range); + tensorstore::TensorStore source; - // std::vector<> + TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( + read_spec, + tensorstore::OpenMode::open, + tensorstore::ReadWriteMode::read).result()); + + auto data_type = GetDataTypeCode(source.dtype().name()); + + switch(data_type){ + case 1: + _write_zarr_chunk(level, channel, y_int, x_int); + break; + case 2: + _write_zarr_chunk(level, channel, y_int, x_int); + break; + case 4: + _write_zarr_chunk(level, channel, y_int, x_int); + break; + case 8: + _write_zarr_chunk(level, channel, y_int, x_int); + break; + case 16: + _write_zarr_chunk(level, channel, y_int, x_int); + break; + case 32: + _write_zarr_chunk(level, channel, y_int, x_int); + break; + case 64: + _write_zarr_chunk(level, channel, y_int, x_int); + break; + case 128: + _write_zarr_chunk(level, channel, y_int, x_int); + break; + case 256: + _write_zarr_chunk(level, channel, y_int, x_int); + break; + case 512: + _write_zarr_chunk(level, channel, y_int, x_int); + break; + default: + break; + } + + std::cerr << "End zarr wrapper" << std::endl; +} + +template +void PyramidCompositor::_write_zarr_chunk(int level, int channel, int y_int, int x_int) { // Compute ranges in global coordinates + std::cerr << "write zarr" << std::endl; + auto image_shape = _zarr_arrays[level].domain().shape(); int y_start = y_int * CHUNK_SIZE; - int y_end = std::min((y_int + 1) * CHUNK_SIZE, (int)std::get<1>(_zarr_arrays[level])[0]); + int y_end = std::min((y_int + 1) * CHUNK_SIZE, (int)image_shape[3]); int x_start = x_int * CHUNK_SIZE; - int x_end = std::min((x_int + 1) * CHUNK_SIZE, (int)std::get<1>(_zarr_arrays[level])[1]); + int x_end = std::min((x_int + 1) * CHUNK_SIZE, (int)image_shape[4]); int assembled_width = x_end - x_start; int assembled_height = y_end - y_start; + std::vector assembled_image_vec(assembled_width*assembled_height); - // auto assembled_image = tensorstore::AllocateArray({ - // assembled_width, - // assembled_height - // }, tensorstore::c_order, - // tensorstore::value_init, - // _image_ts_dtype); - // Determine required input images int unit_image_height = _unit_image_shapes[level].first; int unit_image_width = _unit_image_shapes[level].second; int row_start_pos = y_start; + + std::cerr << "1" << std::endl; while (row_start_pos < y_end) { int row = row_start_pos / unit_image_height; int local_y_start = row_start_pos - y_start; @@ -159,16 +205,16 @@ void PyramidCompositor::_write_zarr_chunk(int level, int channel, int y_int, int int tile_x_start = col_start_pos - col * unit_image_width; int tile_x_dim = std::min((col + 1) * unit_image_width - col_start_pos, x_end - col_start_pos); int tile_x_end = tile_x_start + tile_x_dim; - + std::cerr << "2" << std::endl; // Fetch input file path from composition map std::string input_file_name = _composition_map[{col, row, channel}]; std::filesystem::path zarrArrayLoc = std::filesystem::path(input_file_name) / "data.zarr/0/" / std::to_string(level); - + std::cerr << "3" << std::endl; // Read inage //auto read_data = ReadImageData(zarrArrayLoc, Seq(tile_y_start, tile_y_end), Seq(tile_x_start, tile_x_end)); auto read_spec = GetZarrSpecToRead(zarrArrayLoc.u8string()); - + std::cerr << "4" << std::endl; tensorstore::TensorStore source; TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( @@ -181,11 +227,11 @@ void PyramidCompositor::_write_zarr_chunk(int level, int channel, int y_int, int auto array = tensorstore::Array(read_buffer.data(), {tile_y_end-tile_y_start, tile_x_end-tile_x_start}, tensorstore::c_order); tensorstore::IndexTransform<> read_transform = tensorstore::IdentityTransform(source.domain()); - + std::cerr << "5" << std::endl; Seq rows = Seq(y_start, y_end); Seq cols = Seq(x_start, x_end); Seq tsteps = Seq(0, 0); - Seq channels = Seq(channel, channel); + Seq channels = Seq(channel-1, channel); Seq layers = Seq(0, 0); int x_int=1, y_int=0; @@ -204,18 +250,20 @@ void PyramidCompositor::_write_zarr_chunk(int level, int channel, int y_int, int x_int++; y_int++; } - read_transform = (std::move(read_transform) | tensorstore::Dims(y_int).ClosedInterval(rows.Start(), rows.Stop()) | - tensorstore::Dims(x_int).ClosedInterval(cols.Start(), cols.Stop())).value(); - + read_transform = (std::move(read_transform) | tensorstore::Dims(3).ClosedInterval(rows.Start(), rows.Stop()-1) | + tensorstore::Dims(4).ClosedInterval(cols.Start(), cols.Stop()-1)).value(); + std::cerr << "rows start: " << rows.Start() << ", rows stop: " << rows.Stop() << std::endl; + std::cerr << "rows start: " << cols.Start() << ", rows stop: " << cols.Stop() << std::endl; + std::cerr << "6" << std::endl; tensorstore::Read(source | read_transform, tensorstore::UnownedToShared(array)).value(); - + std::cerr << "7" << std::endl; for (int i = local_y_start; i < local_y_start + tile_y_end - tile_y_start; ++i) { for (int j = local_x_start; j < local_x_start + tile_x_end - tile_x_start; ++j) { assembled_image_vec[i * assembled_width + j] = read_buffer[(i - local_y_start) * (tile_x_end - tile_x_start) + (j - local_x_start)]; } } - + std::cerr << "8" << std::endl; // // auto read_spec = GetZarrSpecToRead(zarrArrayLoc.u8string()); // // tensorstore::TensorStore source; @@ -242,56 +290,16 @@ void PyramidCompositor::_write_zarr_chunk(int level, int channel, int y_int, int // // tensorstore::Dims(3).ClosedInterval(0, image_height - 1) | // // tensorstore::Dims(4).ClosedInterval(0, image_width - 1), // // array).value(); - - //auto read_data = *std::get<0>(ReadImageData(zarrArrayLoc, Seq(tile_y_start, tile_y_end), Seq(tile_x_start, tile_x_end))).get(); - - // Use std::visit to handle the different types in the variant - // std::visit([&](auto& assembled_vec) { - // using T = typename std::decay_t::value_type; - - // for (int i = local_y_start; i < local_y_start + tile_y_end - tile_y_start; ++i) { - // for (int j = local_x_start; j < local_x_start + tile_x_end - tile_x_start; ++j) { - // // Use std::visit to extract the value from read_data - // std::visit([&](auto& value) { - // assembled_vec[i * assembled_width + j] = static_cast(value[(i - local_y_start) * (tile_x_end - tile_x_start) + (j - local_x_start)]); - // }, read_data); - // } - // } - // }, assembled_image); - - - // for (int i = local_y_start; i < local_y_start + tile_y_end - tile_y_start; ++i) { - // for (int j = local_x_start; i < local_x_start + tile_x_end - tile_x_start; ++j) { - // assembled_image[i*(tile_x_end - tile_x_start) + j] = read_data[(i-local_y_start)*(tile_x_end - tile_x_start) + (j-local_x_start)]; - // } - // } - - // // Open the Zarr array for reading with TensorStore - // auto zarrSpec = tensorstore::Context::Default() - // | tensorstore::OpenMode::read - // | tensorstore::dtype_v - // | tensorstore::Shape({1, 1, 1, assembledHeight, assembledWidth}); - - // auto inputZarrArray = tensorstore::Open(zarrSpec).result(); - - // TENSORSTORE_CHECK_OK_AND_ASSIGN() - - // if (!inputZarrArray.ok()) { - // std::cerr << "Failed to open Zarr array: " << inputZarrArray.status() << "\n"; - // return; - // } - - // // Read data from input Zarr file into assembledImage - // tensorstore::Read(*inputZarrArray, - // tensorstore::Span(assembledImage.data(), assembledImage.size())).result(); col_start_pos += tile_x_end - tile_x_start; } row_start_pos += tile_y_end - tile_y_start; } - + std::cerr << "9" << std::endl; + auto& write_source = _zarr_arrays[level]; + std::cerr << "9.1" << std::endl; WriteImageData( - _input_pyramids_loc, + write_source, assembled_image_vec, // Access the data from the std::shared_ptr Seq(y_start, y_end), Seq(x_start, x_end), @@ -300,6 +308,7 @@ void PyramidCompositor::_write_zarr_chunk(int level, int channel, int y_int, int Seq(0, 0) ); + std::cerr << "10" << std::endl; //auto data = *(std::get<0>(_zarr_arrays[level]).get()); // WriteImageData( // _input_pyramids_loc, @@ -325,6 +334,8 @@ void PyramidCompositor::_write_zarr_chunk(int level, int channel, int y_int, int // Seq(0, 0) // ); //}, data); // Dereference the shared_ptr to get the variant + + std::cerr << "end write zarr" << std::endl; } void PyramidCompositor::setComposition(const std::unordered_map, std::string, TupleHash>& comp_map) { @@ -355,9 +366,16 @@ void PyramidCompositor::setComposition(const std::unordered_map(coord.first)); num_cols = std::max(num_cols, std::get<0>(coord.first)); - _num_channels = std::max(num_channels, std::get<2>(coord.first)); + _num_channels = std::max(_num_channels, std::get<2>(coord.first)); } + num_cols += 1; num_rows += 1; - num_channels += 1; - _num_channels = num_channels; + _num_channels += 1; + + std::cerr << "NUM CHANNELS: " << _num_channels << std::endl; + + _plate_image_shapes.clear(); + _zarr_arrays.clear(); + + for (auto& [level, shape]: _unit_image_shapes) { + + _plate_image_shapes[level] = std::vector { + 1, + _num_channels, + 1, + num_rows * shape.first, + num_cols * shape.second + }; + } + for (const auto& [level, shape] : _unit_image_shapes) { + std::string path = _output_pyramid_name.u8string() + "/data.zarr/0/" + std::to_string(level); + + tensorstore::TensorStore source; + + TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( + GetZarrSpecToWrite( + path, + _plate_image_shapes[level], + {1,1,1, CHUNK_SIZE, CHUNK_SIZE}, + ChooseBaseDType(_image_ts_dtype).value().encoded_dtype + ), + tensorstore::OpenMode::create | + tensorstore::OpenMode::delete_existing, + tensorstore::ReadWriteMode::write).result()); + + _zarr_arrays[level] = source; + } + create_auxiliary_files(); + + // TODO: add in code to add size to _zarr_arrays and then use to fix seg fault + // when determining size at beinging of write_zarr // for (const auto& [level, shape] : _unit_image_shapes) { // std::string path = _output_pyramid_name + "/data.zarr/0/" + std::to_string(level); @@ -406,7 +462,7 @@ void PyramidCompositor::reset_composition() { template void PyramidCompositor::WriteImageData( - const std::string& path, + tensorstore::TensorStore& source, std::vector& image, const Seq& rows, const Seq& cols, @@ -416,14 +472,7 @@ void PyramidCompositor::WriteImageData( std::vector shape; - tensorstore::TensorStore source; - - auto result_array = tensorstore::Array(image->data(), {rows.Stop()-rows.Start(), cols.Stop()-cols.Start()}, tensorstore::c_order); - - TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( - GetZarrSpecToRead(path), - tensorstore::OpenMode::open, - tensorstore::ReadWriteMode::write).result()); + auto result_array = tensorstore::Array(image.data(), {rows.Stop()-rows.Start(), cols.Stop()-cols.Start()}, tensorstore::c_order); auto output_transform = tensorstore::IdentityTransform(source.domain()); @@ -442,13 +491,14 @@ void PyramidCompositor::WriteImageData( shape.emplace_back(layers.value().Stop() - layers.value().Start()+1); } - output_transform = (std::move(output_transform) | tensorstore::Dims(_y_index).ClosedInterval(rows.Start(), rows.Stop()) | - tensorstore::Dims(_x_index).ClosedInterval(cols.Start(), cols.Stop())).value(); + output_transform = (std::move(output_transform) | tensorstore::Dims(1).ClosedInterval(channels.value().Start(), channels.value().Stop())).value(); + + output_transform = (std::move(output_transform) | tensorstore::Dims(3).ClosedInterval(rows.Start(), rows.Stop()-1) | + tensorstore::Dims(4).ClosedInterval(cols.Start(), cols.Stop()-1)).value(); shape.emplace_back(rows.Stop() - rows.Start()+1); shape.emplace_back(cols.Stop() - cols.Start()+1); - // Write data array to TensorStore auto write_result = tensorstore::Write(tensorstore::UnownedToShared(result_array), source | output_transform).result(); diff --git a/src/cpp/core/pyramid_compositor.h b/src/cpp/core/pyramid_compositor.h index 8d5029b..bd0d0ee 100644 --- a/src/cpp/core/pyramid_compositor.h +++ b/src/cpp/core/pyramid_compositor.h @@ -105,11 +105,12 @@ class PyramidCompositor { void write_zarr_chunk(int level, int channel, int y_index, int x_index); + void setComposition(const std::unordered_map, std::string, TupleHash>& comp_map); private: template void WriteImageData( - const std::string& path, + tensorstore::TensorStore& source, std::vector& image, const Seq& rows, const Seq& cols, @@ -135,17 +136,15 @@ class PyramidCompositor { image_data GetAssembledImageVector(int size); - void setComposition(const std::unordered_map, std::string, TupleHash>& comp_map); - // members for pyramid composition std::string _input_pyramids_loc; - std::string _output_pyramid_name; + std::filesystem::path _output_pyramid_name; std::string _ome_metadata_file; - std::unordered_map> _plate_image_shapes; + std::unordered_map> _plate_image_shapes; // std::unordered_map> _zarr_arrays; - std::unordered_map, std::vector>> _zarr_arrays; // tuple at 0 is the image and at 1 is the size + //std::unordered_map>> _zarr_arrays; // tuple at 0 is the image and at 1 is the size - //std::unordered_map + std::unordered_map> _zarr_arrays; std::unordered_map> _unit_image_shapes; std::unordered_map, std::string, TupleHash> _composition_map; diff --git a/src/cpp/interface/pyramid_python_interface.cpp b/src/cpp/interface/pyramid_python_interface.cpp index 9f6d731..15d9220 100644 --- a/src/cpp/interface/pyramid_python_interface.cpp +++ b/src/cpp/interface/pyramid_python_interface.cpp @@ -25,7 +25,8 @@ PYBIND11_MODULE(libargolid, m) { .def("create_xml", &argolid::PyramidCompositor::create_xml) \ .def("create_zattr_file", &argolid::PyramidCompositor::create_zattr_file) \ .def("create_auxiliary_files", &argolid::PyramidCompositor::create_auxiliary_files) \ - .def("write_zarr_chunk", &argolid::PyramidCompositor::write_zarr_chunk); + .def("write_zarr_chunk", &argolid::PyramidCompositor::write_zarr_chunk) \ + .def("set_composition", &argolid::PyramidCompositor::setComposition); // .def("write_zarr_chunk", &argolid::PyramidCompositor::write_zarr_chunk) \ // .def("write_zarr_chunk", &argolid::PyramidCompositor::write_zarr_chunk) \ // .def("write_zarr_chunk", &argolid::PyramidCompositor::write_zarr_chunk) \ diff --git a/src/cpp/utilities/utilities.cpp b/src/cpp/utilities/utilities.cpp index e51187a..69db45a 100644 --- a/src/cpp/utilities/utilities.cpp +++ b/src/cpp/utilities/utilities.cpp @@ -205,7 +205,7 @@ void WriteTSZattrFile(const std::string& tiff_file_name, const std::string& zarr void WriteTSZattrFilePlateImage( const std::string& tiff_file_name, const std::string& zarr_root_dir, - const std::unordered_map>& plate_image_shape){ + const std::unordered_map>& plate_image_shape){ json zarr_multiscale_axes; zarr_multiscale_axes = json::parse(R"([ @@ -366,9 +366,14 @@ void GenerateOmeXML(const std::string& image_name, const std::string& output_fil void CreateXML( const std::string& image_name, const std::string& output_file, - const std::tuple& image_shape, + const std::vector& image_shape, const std::string& image_dtype) { + + if (image_shape.size() != 5) { + throw std::invalid_argument("image_shape must have length 5."); + } + pugi::xml_document doc; // Create the root element @@ -393,15 +398,15 @@ void CreateXML( pixelsNode.append_attribute("DimensionOrder") = "XYZCT"; pixelsNode.append_attribute("ID") = "Pixels:0"; pixelsNode.append_attribute("Interleaved") = "false"; - pixelsNode.append_attribute("SizeC") = std::to_string(std::get<1>(image_shape)).c_str(); - pixelsNode.append_attribute("SizeZ") = std::to_string(std::get<2>(image_shape)).c_str(); - pixelsNode.append_attribute("SizeT") = std::to_string(std::get<0>(image_shape)).c_str(); - pixelsNode.append_attribute("SizeX") = std::to_string(std::get<4>(image_shape)).c_str(); - pixelsNode.append_attribute("SizeY") = std::to_string(std::get<3>(image_shape)).c_str(); + pixelsNode.append_attribute("SizeC") = std::to_string(image_shape[1]).c_str(); + pixelsNode.append_attribute("SizeZ") = std::to_string(image_shape[2]).c_str(); + pixelsNode.append_attribute("SizeT") = std::to_string(image_shape[0]).c_str(); + pixelsNode.append_attribute("SizeX") = std::to_string(image_shape[4]).c_str(); + pixelsNode.append_attribute("SizeY") = std::to_string(image_shape[3]).c_str(); pixelsNode.append_attribute("Type") = image_dtype.c_str(); // Create the elements - for(std::int64_t i=0; i(image_shape); ++i){ + for(std::int64_t i=0; i> GetTiffDims (const std:: void CreateXML( const std::string& image_name, const std::string& output_file, - const std::tuple& image_shape, + const std::vector& image_shape, const std::string& image_dtype); void WriteTSZattrFilePlateImage( const std::string& tiff_file_name, const std::string& zarr_root_dir, - const std::unordered_map>& plate_image_shape); + const std::unordered_map>& plate_image_shape); std::tuple, std::optional, std::optional>ParseMultiscaleMetadata(const std::string& axes_list, int len); diff --git a/src/python/argolid/pyramid_compositor.py b/src/python/argolid/pyramid_compositor.py index 0189277..6be477d 100644 --- a/src/python/argolid/pyramid_compositor.py +++ b/src/python/argolid/pyramid_compositor.py @@ -419,7 +419,7 @@ def __init__( output_pyramid_name (str): The name of the zarr pyramid file. """ self._pyramid_compositor_cpp = PyramidCompositorCPP(input_pyramids_loc, out_dir, output_pyramid_name) - + self._chunk_cache: set = set() def _create_xml(self) -> None: """ @@ -468,12 +468,14 @@ def set_composition(self, composition_map: dict) -> None: Args: composition_map (dict): A dictionary mapping composition images to file paths. """ + self._chunk_cache = set() self._pyramid_compositor_cpp.set_composition(composition_map) def reset_composition(self) -> None: """ Resets the pyramid composition by removing the pyramid file and clearing internal data structures. """ + self._chunk_cache = None self._pyramid_compositor_cpp.reset_composition() def get_zarr_chunk( @@ -496,4 +498,25 @@ def get_zarr_chunk( the requested channel does not exist, or the requested y_index or x_index is out of bounds. """ - self._pyramid_compositor_cpp.get_zarr_chunk(level, channel, y_index, x_index) + # TODO: Add this to C++ side + # if self._composition_map is None: + # raise ValueError("No composition map is set. Unable to generate pyramid") + + # if level not in self._unit_image_shapes: + # raise ValueError(f"Requested level ({level}) does not exist") + + # if channel >= self._num_channels: + # raise ValueError(f"Requested channel ({channel}) does not exist") + + # if y_index > (self._plate_image_shapes[level][3] // CHUNK_SIZE): + # raise ValueError(f"Requested y index ({y_index}) does not exist") + + # if x_index > (self._plate_image_shapes[level][4] // CHUNK_SIZE): + # raise ValueError(f"Requested y index ({x_index}) does not exist") + + if (level, channel, y_index, x_index) in self._chunk_cache: + return + else: + self._pyramid_compositor_cpp.write_zarr_chunk(level, channel, y_index, x_index) + self._chunk_cache.add((level, channel, y_index, x_index)) + return From 2cfee3dd62b8c03270baf469238dfe5126fa0f84 Mon Sep 17 00:00:00 2001 From: JesseMckinzie Date: Wed, 22 Jan 2025 08:53:44 -0700 Subject: [PATCH 05/24] Clean cpp pyramid composition --- src/cpp/core/pyramid_compositor.cpp | 380 +----------------- src/cpp/core/pyramid_compositor.h | 29 +- .../interface/pyramid_python_interface.cpp | 12 +- 3 files changed, 28 insertions(+), 393 deletions(-) diff --git a/src/cpp/core/pyramid_compositor.cpp b/src/cpp/core/pyramid_compositor.cpp index 5fab267..9085621 100644 --- a/src/cpp/core/pyramid_compositor.cpp +++ b/src/cpp/core/pyramid_compositor.cpp @@ -26,91 +26,37 @@ PyramidCompositor::PyramidCompositor(const std::string& input_pyramids_loc, cons _output_pyramid_name(std::filesystem::path(out_dir) / output_pyramid_name), _ome_metadata_file(out_dir + "/" + output_pyramid_name + "/METADATA.ome.xml"), _pyramid_levels(-1), - _num_channels(-1) { + _num_channels(-1) {} - std::cerr << "Constructor" << std::endl; - // tensorstore::TensorStore source; - // TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( - // GetZarrSpecToRead(input_pyramids_loc), - // tensorstore::OpenMode::open, - // tensorstore::ReadWriteMode::read).result()); - - // tensorstore::DataType _image_ts_dtype = source.dtype(); - // _image_dtype = source.dtype().name(); - // _image_dtype_code = GetDataTypeCode(_image_dtype); - - // auto read_spec = GetZarrSpecToRead(input_pyramids_loc); - - // TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( - // read_spec, - // tensorstore::OpenMode::open, - // tensorstore::ReadWriteMode::read).result()); - - // auto image_shape = source.domain().shape(); - // const auto read_chunk_shape = source.chunk_layout().value().read_chunk_shape(); - - // if (image_shape.size() == 5){ - // _t_index.emplace(0); - // _c_index.emplace(1); - // _z_index.emplace(2); - // } else { - // assert(image_shape.size() >= 2); - // std::tie(_t_index, _c_index, _z_index) = ParseMultiscaleMetadata("XYCTZ", image_shape.size()); - // } - - // _data_type = source.dtype().name(); - // _data_type_code = GetDataTypeCode(_data_type); - - // TENSORSTORE_CHECK_OK_AND_ASSIGN(___mb_cur_max, tensorstore::Open( - // GetZarrSpecToWrite(_filename, _image_shape, _chunk_shape, GetEncodedType(_dtype_code)), - // tensorstore::OpenMode::create | - // tensorstore::OpenMode::delete_existing, - // tensorstore::ReadWriteMode::write).result() - // ); - - std::cerr << "End Constructor" << std::endl; -} - -// Private Methods void PyramidCompositor::create_xml() { - std::cerr << "XML" << std::endl; CreateXML( this->_output_pyramid_name, this->_ome_metadata_file, this->_plate_image_shapes[0], this->_image_dtype ); - std::cerr << "End XML" << std::endl; } void PyramidCompositor::create_zattr_file() { - std::cerr << "zattr" << std::endl; WriteTSZattrFilePlateImage( this->_output_pyramid_name, this->_output_pyramid_name.u8string() + "/data.zarr/0", this->_plate_image_shapes ); - - std::cerr << "End zattr" << std::endl; } void PyramidCompositor::create_zgroup_file() { - std::cerr << "zgroup" << std::endl; WriteVivZgroupFiles(this->_output_pyramid_name); - std::cerr << "End zgroup" << std::endl; } void PyramidCompositor::create_auxiliary_files() { - std::cerr << "aux" << std::endl; create_xml(); create_zattr_file(); create_zgroup_file(); - std::cerr << "End aux" << std::endl; } void PyramidCompositor::write_zarr_chunk(int level, int channel, int y_int, int x_int) { - std::cerr << "zarr wrapper" << std::endl; // get datatype by opening std::string input_file_name = _composition_map[{x_int, y_int, channel}]; @@ -118,8 +64,6 @@ void PyramidCompositor::write_zarr_chunk(int level, int channel, int y_int, int auto read_spec = GetZarrSpecToRead(zarrArrayLoc.u8string()); - std::cerr << "Read path 121: " << zarrArrayLoc.u8string() << std::endl; - tensorstore::TensorStore source; TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( @@ -163,15 +107,12 @@ void PyramidCompositor::write_zarr_chunk(int level, int channel, int y_int, int default: break; } - - std::cerr << "End zarr wrapper" << std::endl; } template void PyramidCompositor::_write_zarr_chunk(int level, int channel, int y_int, int x_int) { // Compute ranges in global coordinates - std::cerr << "write zarr" << std::endl; auto image_shape = _zarr_arrays[level].domain().shape(); int y_start = y_int * CHUNK_SIZE; int y_end = std::min((y_int + 1) * CHUNK_SIZE, (int)image_shape[3]); @@ -180,7 +121,6 @@ void PyramidCompositor::_write_zarr_chunk(int level, int channel, int y_int, int int assembled_width = x_end - x_start; int assembled_height = y_end - y_start; - std::vector assembled_image_vec(assembled_width*assembled_height); @@ -190,8 +130,8 @@ void PyramidCompositor::_write_zarr_chunk(int level, int channel, int y_int, int int row_start_pos = y_start; - std::cerr << "1" << std::endl; while (row_start_pos < y_end) { + int row = row_start_pos / unit_image_height; int local_y_start = row_start_pos - y_start; int tile_y_start = row_start_pos - row * unit_image_height; @@ -199,22 +139,21 @@ void PyramidCompositor::_write_zarr_chunk(int level, int channel, int y_int, int int tile_y_end = tile_y_start + tile_y_dim; int col_start_pos = x_start; + while (col_start_pos < x_end) { + int col = col_start_pos / unit_image_width; int local_x_start = col_start_pos - x_start; int tile_x_start = col_start_pos - col * unit_image_width; int tile_x_dim = std::min((col + 1) * unit_image_width - col_start_pos, x_end - col_start_pos); int tile_x_end = tile_x_start + tile_x_dim; - std::cerr << "2" << std::endl; + // Fetch input file path from composition map std::string input_file_name = _composition_map[{col, row, channel}]; std::filesystem::path zarrArrayLoc = std::filesystem::path(input_file_name) / "data.zarr/0/" / std::to_string(level); - std::cerr << "3" << std::endl; - // Read inage - //auto read_data = ReadImageData(zarrArrayLoc, Seq(tile_y_start, tile_y_end), Seq(tile_x_start, tile_x_end)); auto read_spec = GetZarrSpecToRead(zarrArrayLoc.u8string()); - std::cerr << "4" << std::endl; + tensorstore::TensorStore source; TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( @@ -222,84 +161,37 @@ void PyramidCompositor::_write_zarr_chunk(int level, int channel, int y_int, int tensorstore::OpenMode::open, tensorstore::ReadWriteMode::read).result()); - std::vector read_buffer((tile_x_end-tile_x_start)*(tile_y_end-tile_y_start)); auto array = tensorstore::Array(read_buffer.data(), {tile_y_end-tile_y_start, tile_x_end-tile_x_start}, tensorstore::c_order); tensorstore::IndexTransform<> read_transform = tensorstore::IdentityTransform(source.domain()); - std::cerr << "5" << std::endl; + Seq rows = Seq(y_start, y_end); Seq cols = Seq(x_start, x_end); Seq tsteps = Seq(0, 0); Seq channels = Seq(channel-1, channel); Seq layers = Seq(0, 0); + + read_transform = (std::move(read_transform) | tensorstore::Dims(_y_index).ClosedInterval(rows.Start(), rows.Stop()-1) | + tensorstore::Dims(_x_index).ClosedInterval(cols.Start(), cols.Stop()-1)).value(); - int x_int=1, y_int=0; - if (_t_index.has_value()){ - read_transform = (std::move(read_transform) | tensorstore::Dims(_t_index.value()).ClosedInterval(tsteps.Start(), tsteps.Stop())).value(); - x_int++; - y_int++; - } - if (_c_index.has_value()){ - read_transform = (std::move(read_transform) | tensorstore::Dims(_c_index.value()).ClosedInterval(channels.Start(), channels.Stop())).value(); - x_int++; - y_int++; - } - if (_z_index.has_value()){ - read_transform = (std::move(read_transform) | tensorstore::Dims(_z_index.value()).ClosedInterval(layers.Start(), layers.Stop())).value(); - x_int++; - y_int++; - } - read_transform = (std::move(read_transform) | tensorstore::Dims(3).ClosedInterval(rows.Start(), rows.Stop()-1) | - tensorstore::Dims(4).ClosedInterval(cols.Start(), cols.Stop()-1)).value(); - - std::cerr << "rows start: " << rows.Start() << ", rows stop: " << rows.Stop() << std::endl; - std::cerr << "rows start: " << cols.Start() << ", rows stop: " << cols.Stop() << std::endl; - std::cerr << "6" << std::endl; tensorstore::Read(source | read_transform, tensorstore::UnownedToShared(array)).value(); - std::cerr << "7" << std::endl; + for (int i = local_y_start; i < local_y_start + tile_y_end - tile_y_start; ++i) { for (int j = local_x_start; j < local_x_start + tile_x_end - tile_x_start; ++j) { assembled_image_vec[i * assembled_width + j] = read_buffer[(i - local_y_start) * (tile_x_end - tile_x_start) + (j - local_x_start)]; } } - std::cerr << "8" << std::endl; - // // auto read_spec = GetZarrSpecToRead(zarrArrayLoc.u8string()); - - // // tensorstore::TensorStore source; - - // // TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( - // // read_spec, - // // tensorstore::OpenMode::open, - // // tensorstore::ReadWriteMode::read).result()); - - // // auto image_shape = source.domain().shape(); - // // const auto read_chunk_shape = source.chunk_layout().value().read_chunk_shape(); - - // // auto image_width = image_shape[4]; - // // auto image_height = image_shape[3]; - - // // auto array = tensorstore::AllocateArray({ - // // image_height, - // // image_width - // // }, tensorstore::c_order, - // // tensorstore::value_init, source.dtype()); - - // // // initiate a read - // // tensorstore::Read(source | - // // tensorstore::Dims(3).ClosedInterval(0, image_height - 1) | - // // tensorstore::Dims(4).ClosedInterval(0, image_width - 1), - // // array).value(); col_start_pos += tile_x_end - tile_x_start; } row_start_pos += tile_y_end - tile_y_start; } - std::cerr << "9" << std::endl; - auto& write_source = _zarr_arrays[level]; - std::cerr << "9.1" << std::endl; + + //auto& write_source = _zarr_arrays[level]; + WriteImageData( - write_source, + _zarr_arrays[level], assembled_image_vec, // Access the data from the std::shared_ptr Seq(y_start, y_end), Seq(x_start, x_end), @@ -307,38 +199,9 @@ void PyramidCompositor::_write_zarr_chunk(int level, int channel, int y_int, int Seq(channel, channel), Seq(0, 0) ); - - std::cerr << "10" << std::endl; - //auto data = *(std::get<0>(_zarr_arrays[level]).get()); - // WriteImageData( - // _input_pyramids_loc, - // data, - // Seq(y_start,y_end), Seq(x_start, x_end), - // Seq(0,0), Seq(channel, channel), Seq(0,0) - // ); - - // Assuming _zarr_arrays[level] is of type std::tuple>, std::vector> - //auto& [variant_ptr, some_vector] = _zarr_arrays[level]; // Structured binding (C++17) - - // std::visit([&](const auto& array) { - // using T = std::decay_t; - - // // Use external variables and call WriteImageData - // WriteImageData( - // _input_pyramids_loc, - // array, // Access the data from the std::shared_ptr - // Seq(y_start, y_end), - // Seq(x_start, x_end), - // Seq(0, 0), - // Seq(channel, channel), - // Seq(0, 0) - // ); - //}, data); // Dereference the shared_ptr to get the variant - - std::cerr << "end write zarr" << std::endl; } -void PyramidCompositor::setComposition(const std::unordered_map, std::string, TupleHash>& comp_map) { +void PyramidCompositor::set_composition(const std::unordered_map, std::string, TupleHash>& comp_map) { _composition_map = comp_map; for (const auto& coord : comp_map) { @@ -370,12 +233,6 @@ void PyramidCompositor::setComposition(const std::unordered_map(zarr_array_data); - // } + create_auxiliary_files(); } void PyramidCompositor::reset_composition() { - // Remove pyramid file and clear internal data structures fs::remove_all(_output_pyramid_name); _composition_map.clear(); _plate_image_shapes.clear(); @@ -476,25 +314,10 @@ void PyramidCompositor::WriteImageData( auto output_transform = tensorstore::IdentityTransform(source.domain()); - if (_t_index.has_value() && tsteps.has_value()) { - output_transform = (std::move(output_transform) | tensorstore::Dims(_t_index.value()).ClosedInterval(tsteps.value().Start(), tsteps.value().Stop())).value(); - shape.emplace_back(tsteps.value().Stop() - tsteps.value().Start()+1); - } + output_transform = (std::move(output_transform) | tensorstore::Dims(_c_index).ClosedInterval(channels.value().Start(), channels.value().Stop())).value(); - if (_c_index.has_value() && channels.has_value()) { - output_transform = (std::move(output_transform) | tensorstore::Dims(_c_index.value()).ClosedInterval(channels.value().Start(), channels.value().Stop())).value(); - shape.emplace_back(channels.value().Stop() - channels.value().Start()+1); - } - - if (_z_index.has_value() && layers.has_value()) { - output_transform = (std::move(output_transform) | tensorstore::Dims(_z_index.value()).ClosedInterval(layers.value().Start(), layers.value().Stop())).value(); - shape.emplace_back(layers.value().Stop() - layers.value().Start()+1); - } - - output_transform = (std::move(output_transform) | tensorstore::Dims(1).ClosedInterval(channels.value().Start(), channels.value().Stop())).value(); - - output_transform = (std::move(output_transform) | tensorstore::Dims(3).ClosedInterval(rows.Start(), rows.Stop()-1) | - tensorstore::Dims(4).ClosedInterval(cols.Start(), cols.Stop()-1)).value(); + output_transform = (std::move(output_transform) | tensorstore::Dims(_y_index).ClosedInterval(rows.Start(), rows.Stop()-1) | + tensorstore::Dims(_x_index).ClosedInterval(cols.Start(), cols.Stop()-1)).value(); shape.emplace_back(rows.Stop() - rows.Start()+1); shape.emplace_back(cols.Stop() - cols.Start()+1); @@ -506,167 +329,4 @@ void PyramidCompositor::WriteImageData( std::cerr << "Error writing image: " << write_result.status() << std::endl; } } - - -auto PyramidCompositor::ReadImageData( - const std::string& fname, - const Seq& rows, - const Seq& cols, - const Seq& layers, - const Seq& channels, - const Seq& tsteps) { - - auto read_spec = GetZarrSpecToRead(fname); - - tensorstore::TensorStore source; - - TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( - read_spec, - tensorstore::OpenMode::open, - tensorstore::ReadWriteMode::read).result()); - - auto image_shape = source.domain().shape(); - const auto read_chunk_shape = source.chunk_layout().value().read_chunk_shape(); - - const auto data_height = rows.Stop() - rows.Start() + 1; - const auto data_width = cols.Stop() - cols.Start() + 1; - const auto data_depth = layers.Stop() - layers.Start() + 1; - const auto data_num_channels = channels.Stop() - channels.Start() + 1; - const auto data_tsteps = tsteps.Stop() - tsteps.Start() + 1; - - auto array = tensorstore::AllocateArray( - {data_tsteps, data_num_channels, data_depth, data_height, data_width}, - tensorstore::c_order, - tensorstore::value_init, - source.dtype() - ); - - tensorstore::IndexTransform<> read_transform = tensorstore::IdentityTransform(source.domain()); - - int x_int=1, y_int=0; - if (_t_index.has_value()){ - read_transform = (std::move(read_transform) | tensorstore::Dims(_t_index.value()).ClosedInterval(tsteps.Start(), tsteps.Stop())).value(); - x_int++; - y_int++; - } - if (_c_index.has_value()){ - read_transform = (std::move(read_transform) | tensorstore::Dims(_c_index.value()).ClosedInterval(channels.Start(), channels.Stop())).value(); - x_int++; - y_int++; - } - if (_z_index.has_value()){ - read_transform = (std::move(read_transform) | tensorstore::Dims(_z_index.value()).ClosedInterval(layers.Start(), layers.Stop())).value(); - x_int++; - y_int++; - } - read_transform = (std::move(read_transform) | tensorstore::Dims(y_int).ClosedInterval(rows.Start(), rows.Stop()) | - tensorstore::Dims(x_int).ClosedInterval(cols.Start(), cols.Stop())).value(); - - - tensorstore::Read(source | read_transform, tensorstore::UnownedToShared(array)).value(); - - return array; -} - -// std::tuple, std::vector> PyramidCompositor::ReadImageData( -// const std::string& fname, -// const Seq& rows, const Seq& cols, -// const Seq& layers, -// const Seq& channels, -// const Seq& tsteps) { - -// auto read_spec = GetZarrSpecToRead(fname); - -// tensorstore::TensorStore source; - -// TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( -// read_spec, -// tensorstore::OpenMode::open, -// tensorstore::ReadWriteMode::read).result()); - -// auto image_shape = source.domain().shape(); -// const auto read_chunk_shape = source.chunk_layout().value().read_chunk_shape(); - -// // if (image_shape.size() == 5){ -// // _image_height = image_shape[3]; -// // _image_width = image_shape[4]; -// // _image_depth = image_shape[2]; -// // _num_channels = image_shape[1]; -// // _num_tsteps = image_shape[0]; -// // _t_index.emplace(0); -// // _c_index.emplace(1); -// // _z_index.emplace(2); -// // } else { -// // assert(image_shape.size() >= 2); -// // std::tie(_t_index, _c_index, _z_index) = ParseMultiscaleMetadata(axes_list, image_shape.size()); -// // _image_height = image_shape[image_shape.size()-2]; -// // _image_width = image_shape[image_shape.size()-1]; -// // if (_t_index.has_value()) { -// // _num_tsteps = image_shape[_t_index.value()]; -// // } else { -// // _num_tsteps = 1; -// // } - -// // if (_c_index.has_value()) { -// // _num_channels = image_shape[_c_index.value()]; -// // } else { -// // _num_channels = 1; -// // } - -// // if (_z_index.has_value()) { -// // _image_depth = image_shape[_z_index.value()]; -// // } else { -// // _image_depth = 1; -// // } -// // } - -// _image_dtype = source.dtype().name(); -// _image_dtype_code = GetDataTypeCode(_image_dtype); - -// switch (_data_type_code) -// { -// case (1): -// return std::make_tuple( -// std::make_shared(std::move(*(GetImageDataTemplated(source, rows, cols, layers, channels, tsteps)))), -// std::vector{image_shape[image_shape.size()-2], image_shape[image_shape.size()-1]}); -// case (2): -// return std::make_tuple( -// std::make_shared(std::move(*(GetImageDataTemplated(source, rows, cols, layers, channels, tsteps)))), -// std::vector{image_shape[image_shape.size()-2], image_shape[image_shape.size()-1]}); -// case (4): -// return std::make_tuple( -// std::make_shared(std::move(*(GetImageDataTemplated(source, rows, cols, layers, channels, tsteps)))), -// std::vector{image_shape[image_shape.size()-2], image_shape[image_shape.size()-1]}); -// case (8): -// return std::make_tuple( -// std::make_shared(std::move(*(GetImageDataTemplated(source, rows, cols, layers, channels, tsteps)))), -// std::vector{image_shape[image_shape.size()-2], image_shape[image_shape.size()-1]}); -// case (16): -// return std::make_tuple( -// std::make_shared(std::move(*(GetImageDataTemplated(source, rows, cols, layers, channels, tsteps)))), -// std::vector{image_shape[image_shape.size()-2], image_shape[image_shape.size()-1]}); -// case (32): -// return std::make_tuple( -// std::make_shared(std::move(*(GetImageDataTemplated(source, rows, cols, layers, channels, tsteps)))), -// std::vector{image_shape[image_shape.size()-2], image_shape[image_shape.size()-1]}); -// case (64): -// return std::make_tuple( -// std::make_shared(std::move(*(GetImageDataTemplated(source, rows, cols, layers, channels, tsteps)))), -// std::vector{image_shape[image_shape.size()-2], image_shape[image_shape.size()-1]}); -// case (128): -// return std::make_tuple( -// std::make_shared(std::move(*(GetImageDataTemplated(source, rows, cols, layers, channels, tsteps)))), -// std::vector{image_shape[image_shape.size()-2], image_shape[image_shape.size()-1]}); -// case (256): -// return std::make_tuple( -// std::make_shared(std::move(*(GetImageDataTemplated(source, rows, cols, layers, channels, tsteps)))), -// std::vector{image_shape[image_shape.size()-2], image_shape[image_shape.size()-1]}); -// case (512): -// return std::make_tuple( -// std::make_shared(std::move(*(GetImageDataTemplated(source, rows, cols, layers, channels, tsteps)))), -// std::vector{image_shape[image_shape.size()-2], image_shape[image_shape.size()-1]}); -// default: -// break; -// } -// } } // ns argolid \ No newline at end of file diff --git a/src/cpp/core/pyramid_compositor.h b/src/cpp/core/pyramid_compositor.h index bd0d0ee..e10ef15 100644 --- a/src/cpp/core/pyramid_compositor.h +++ b/src/cpp/core/pyramid_compositor.h @@ -105,7 +105,8 @@ class PyramidCompositor { void write_zarr_chunk(int level, int channel, int y_index, int x_index); - void setComposition(const std::unordered_map, std::string, TupleHash>& comp_map); + void set_composition(const std::unordered_map, std::string, TupleHash>& comp_map); + private: template @@ -118,35 +119,18 @@ class PyramidCompositor { const std::optional& channels, const std::optional& tsteps); - auto ReadImageData( - const std::string& fname, - const Seq& rows, - const Seq& cols, - const Seq& layers = Seq(0,0), - const Seq& channels = Seq(0,0), - const Seq& tsteps = Seq(0,0)); - - // std::tuple, std::vector> ReadImageData( - // const std::string& fname, - // const Seq& rows, const Seq& cols, - // const Seq& layers = Seq(0,0), - // const Seq& channels = Seq(0,0), - // const Seq& tsteps = Seq(0,0) - // ); - image_data GetAssembledImageVector(int size); - // members for pyramid composition std::string _input_pyramids_loc; std::filesystem::path _output_pyramid_name; std::string _ome_metadata_file; + std::unordered_map> _plate_image_shapes; - // std::unordered_map> _zarr_arrays; - //std::unordered_map>> _zarr_arrays; // tuple at 0 is the image and at 1 is the size std::unordered_map> _zarr_arrays; std::unordered_map> _unit_image_shapes; + std::unordered_map, std::string, TupleHash> _composition_map; std::set> _tile_cache; @@ -157,7 +141,8 @@ class PyramidCompositor { std::string _image_dtype; std::uint16_t _image_dtype_code; - std::optional_z_index, _c_index, _t_index; - int _x_index, _y_index; + int _x_index = 4; + int _y_index = 3; + int _c_index = 1; }; } // ns argolid diff --git a/src/cpp/interface/pyramid_python_interface.cpp b/src/cpp/interface/pyramid_python_interface.cpp index 15d9220..29c81d3 100644 --- a/src/cpp/interface/pyramid_python_interface.cpp +++ b/src/cpp/interface/pyramid_python_interface.cpp @@ -26,17 +26,7 @@ PYBIND11_MODULE(libargolid, m) { .def("create_zattr_file", &argolid::PyramidCompositor::create_zattr_file) \ .def("create_auxiliary_files", &argolid::PyramidCompositor::create_auxiliary_files) \ .def("write_zarr_chunk", &argolid::PyramidCompositor::write_zarr_chunk) \ - .def("set_composition", &argolid::PyramidCompositor::setComposition); - // .def("write_zarr_chunk", &argolid::PyramidCompositor::write_zarr_chunk) \ - // .def("write_zarr_chunk", &argolid::PyramidCompositor::write_zarr_chunk) \ - // .def("write_zarr_chunk", &argolid::PyramidCompositor::write_zarr_chunk) \ - // .def("write_zarr_chunk", &argolid::PyramidCompositor::write_zarr_chunk) \ - // .def("write_zarr_chunk", &argolid::PyramidCompositor::write_zarr_chunk) \ - // .def("write_zarr_chunk", &argolid::PyramidCompositor::write_zarr_chunk) \ - // .def("write_zarr_chunk", &argolid::PyramidCompositor::write_zarr_chunk) \ - // .def("write_zarr_chunk", &argolid::PyramidCompositor::write_zarr_chunk) \ - // .def("write_zarr_chunk", &argolid::PyramidCompositor::write_zarr_chunk) \ - // .def("write_zarr_chunk", &argolid::PyramidCompositor::write_zarr_chunk) ; + .def("set_composition", &argolid::PyramidCompositor::set_composition); py::enum_(m, "VisType") .value("NG_Zarr", argolid::VisType::NG_Zarr) From d1d5d557ed085b9b13f4e516413eb89247aff4d6 Mon Sep 17 00:00:00 2001 From: JesseMckinzie Date: Wed, 22 Jan 2025 09:56:38 -0700 Subject: [PATCH 06/24] Move error handling and cache to cpp compositor --- src/cpp/core/pyramid_compositor.cpp | 73 ++++++++++++++++++------ src/cpp/core/pyramid_compositor.h | 4 +- src/python/argolid/pyramid_compositor.py | 23 +------- 3 files changed, 56 insertions(+), 44 deletions(-) diff --git a/src/cpp/core/pyramid_compositor.cpp b/src/cpp/core/pyramid_compositor.cpp index 9085621..9cb9fb7 100644 --- a/src/cpp/core/pyramid_compositor.cpp +++ b/src/cpp/core/pyramid_compositor.cpp @@ -56,10 +56,42 @@ void PyramidCompositor::create_auxiliary_files() { create_zgroup_file(); } -void PyramidCompositor::write_zarr_chunk(int level, int channel, int y_int, int x_int) { +void PyramidCompositor::write_zarr_chunk(int level, int channel, int y_index, int x_index) { + + std::tuple chunk = std::make_tuple(level, channel, y_index, x_index); + + if (_chunk_cache.find(chunk) != _chunk_cache.end()) { + return; + } + + if (_composition_map.empty()) { + throw std::runtime_error("No composition map is set. Unable to generate pyramid"); + } + + if (_unit_image_shapes.find(level) == _unit_image_shapes.end()) { + throw std::invalid_argument("Requested level (" + std::to_string(level) + ") does not exist"); + } + + if (channel >= _num_channels) { + throw std::invalid_argument("Requested channel (" + std::to_string(channel) + ") does not exist"); + } + + auto plate_shape_it = _plate_image_shapes.find(level); + if (plate_shape_it == _plate_image_shapes.end()) { + throw std::invalid_argument("No plate image shapes found for level (" + std::to_string(level) + ")"); + } + + const auto& plate_shape = plate_shape_it->second; + if (plate_shape.size() < 4 || y_index > (plate_shape[3] / CHUNK_SIZE)) { + throw std::invalid_argument("Requested y index (" + std::to_string(y_index) + ") does not exist"); + } + + if (plate_shape.size() < 5 || x_index > (plate_shape[4] / CHUNK_SIZE)) { + throw std::invalid_argument("Requested x index (" + std::to_string(x_index) + ") does not exist"); + } // get datatype by opening - std::string input_file_name = _composition_map[{x_int, y_int, channel}]; + std::string input_file_name = _composition_map[{x_index, y_index, channel}]; std::filesystem::path zarrArrayLoc = _output_pyramid_name / input_file_name / std::filesystem::path("data.zarr/0/") / std::to_string(level); auto read_spec = GetZarrSpecToRead(zarrArrayLoc.u8string()); @@ -75,49 +107,51 @@ void PyramidCompositor::write_zarr_chunk(int level, int channel, int y_int, int switch(data_type){ case 1: - _write_zarr_chunk(level, channel, y_int, x_int); + _write_zarr_chunk(level, channel, y_index, x_index); break; case 2: - _write_zarr_chunk(level, channel, y_int, x_int); + _write_zarr_chunk(level, channel, y_index, x_index); break; case 4: - _write_zarr_chunk(level, channel, y_int, x_int); + _write_zarr_chunk(level, channel, y_index, x_index); break; case 8: - _write_zarr_chunk(level, channel, y_int, x_int); + _write_zarr_chunk(level, channel, y_index, x_index); break; case 16: - _write_zarr_chunk(level, channel, y_int, x_int); + _write_zarr_chunk(level, channel, y_index, x_index); break; case 32: - _write_zarr_chunk(level, channel, y_int, x_int); + _write_zarr_chunk(level, channel, y_index, x_index); break; case 64: - _write_zarr_chunk(level, channel, y_int, x_int); + _write_zarr_chunk(level, channel, y_index, x_index); break; case 128: - _write_zarr_chunk(level, channel, y_int, x_int); + _write_zarr_chunk(level, channel, y_index, x_index); break; case 256: - _write_zarr_chunk(level, channel, y_int, x_int); + _write_zarr_chunk(level, channel, y_index, x_index); break; case 512: - _write_zarr_chunk(level, channel, y_int, x_int); + _write_zarr_chunk(level, channel, y_index, x_index); break; default: break; } + + _chunk_cache.insert(chunk); } template -void PyramidCompositor::_write_zarr_chunk(int level, int channel, int y_int, int x_int) { +void PyramidCompositor::_write_zarr_chunk(int level, int channel, int y_index, int x_index) { // Compute ranges in global coordinates auto image_shape = _zarr_arrays[level].domain().shape(); - int y_start = y_int * CHUNK_SIZE; - int y_end = std::min((y_int + 1) * CHUNK_SIZE, (int)image_shape[3]); - int x_start = x_int * CHUNK_SIZE; - int x_end = std::min((x_int + 1) * CHUNK_SIZE, (int)image_shape[4]); + int y_start = y_index * CHUNK_SIZE; + int y_end = std::min((y_index + 1) * CHUNK_SIZE, (int)image_shape[3]); + int x_start = x_index * CHUNK_SIZE; + int x_end = std::min((x_index + 1) * CHUNK_SIZE, (int)image_shape[4]); int assembled_width = x_end - x_start; int assembled_height = y_end - y_start; @@ -242,7 +276,8 @@ void PyramidCompositor::set_composition(const std::unordered_map(coord.first)); num_cols = std::max(num_cols, std::get<0>(coord.first)); @@ -293,7 +328,7 @@ void PyramidCompositor::reset_composition() { fs::remove_all(_output_pyramid_name); _composition_map.clear(); _plate_image_shapes.clear(); - _tile_cache.clear(); + _chunk_cache.clear(); _plate_image_shapes.clear(); _zarr_arrays.clear(); } diff --git a/src/cpp/core/pyramid_compositor.h b/src/cpp/core/pyramid_compositor.h index e10ef15..6ef411c 100644 --- a/src/cpp/core/pyramid_compositor.h +++ b/src/cpp/core/pyramid_compositor.h @@ -119,8 +119,6 @@ class PyramidCompositor { const std::optional& channels, const std::optional& tsteps); - image_data GetAssembledImageVector(int size); - std::string _input_pyramids_loc; std::filesystem::path _output_pyramid_name; std::string _ome_metadata_file; @@ -133,7 +131,7 @@ class PyramidCompositor { std::unordered_map, std::string, TupleHash> _composition_map; - std::set> _tile_cache; + std::set> _chunk_cache; int _pyramid_levels; int _num_channels; diff --git a/src/python/argolid/pyramid_compositor.py b/src/python/argolid/pyramid_compositor.py index 6be477d..ae5e5f7 100644 --- a/src/python/argolid/pyramid_compositor.py +++ b/src/python/argolid/pyramid_compositor.py @@ -468,14 +468,12 @@ def set_composition(self, composition_map: dict) -> None: Args: composition_map (dict): A dictionary mapping composition images to file paths. """ - self._chunk_cache = set() self._pyramid_compositor_cpp.set_composition(composition_map) def reset_composition(self) -> None: """ Resets the pyramid composition by removing the pyramid file and clearing internal data structures. """ - self._chunk_cache = None self._pyramid_compositor_cpp.reset_composition() def get_zarr_chunk( @@ -498,25 +496,6 @@ def get_zarr_chunk( the requested channel does not exist, or the requested y_index or x_index is out of bounds. """ - # TODO: Add this to C++ side - # if self._composition_map is None: - # raise ValueError("No composition map is set. Unable to generate pyramid") - - # if level not in self._unit_image_shapes: - # raise ValueError(f"Requested level ({level}) does not exist") - # if channel >= self._num_channels: - # raise ValueError(f"Requested channel ({channel}) does not exist") - - # if y_index > (self._plate_image_shapes[level][3] // CHUNK_SIZE): - # raise ValueError(f"Requested y index ({y_index}) does not exist") - - # if x_index > (self._plate_image_shapes[level][4] // CHUNK_SIZE): - # raise ValueError(f"Requested y index ({x_index}) does not exist") + self._pyramid_compositor_cpp.write_zarr_chunk(level, channel, y_index, x_index) - if (level, channel, y_index, x_index) in self._chunk_cache: - return - else: - self._pyramid_compositor_cpp.write_zarr_chunk(level, channel, y_index, x_index) - self._chunk_cache.add((level, channel, y_index, x_index)) - return From abbd55c40b4702b0c97ca70853c8f9eb032a61f2 Mon Sep 17 00:00:00 2001 From: JesseMckinzie Date: Wed, 22 Jan 2025 10:20:22 -0700 Subject: [PATCH 07/24] Clear cache in set_composition --- src/cpp/core/pyramid_compositor.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/cpp/core/pyramid_compositor.cpp b/src/cpp/core/pyramid_compositor.cpp index 9cb9fb7..7a1cde3 100644 --- a/src/cpp/core/pyramid_compositor.cpp +++ b/src/cpp/core/pyramid_compositor.cpp @@ -290,6 +290,7 @@ void PyramidCompositor::set_composition(const std::unordered_map Date: Wed, 22 Jan 2025 10:20:48 -0700 Subject: [PATCH 08/24] Remove python compositor --- src/python/argolid/pyramid_compositor.py | 394 +---------------------- 1 file changed, 1 insertion(+), 393 deletions(-) diff --git a/src/python/argolid/pyramid_compositor.py b/src/python/argolid/pyramid_compositor.py index ae5e5f7..df8eaa3 100644 --- a/src/python/argolid/pyramid_compositor.py +++ b/src/python/argolid/pyramid_compositor.py @@ -10,399 +10,8 @@ from .libargolid import PyramidCompositorCPP -CHUNK_SIZE: int = 1024 -OME_DTYPE = { - "uint8": ome_types.model.PixelType.UINT8, - "int8": ome_types.model.PixelType.INT8, - "uint16": ome_types.model.PixelType.UINT16, - "int16": ome_types.model.PixelType.INT16, - "uint32": ome_types.model.PixelType.UINT32, - "int32": ome_types.model.PixelType.INT32, - "float": ome_types.model.PixelType.FLOAT, - "double": ome_types.model.PixelType.DOUBLE, -} - - -def get_zarr_read_spec(file_path: str) -> dict: - """ - Returns a dictionary containing the specification for reading a Zarr file. - - Args: - file_path (str): The path to the Zarr file. - - Returns: - dict: A dictionary containing the specification for reading the Zarr file. - """ - return { - "driver": "zarr", - "kvstore": { - "driver": "file", - "path": file_path, - }, - "open": True, - "context": { - "cache_pool": {}, - "data_copy_concurrency": {"limit": os.cpu_count() // 2}, - "file_io_concurrency": {"limit": os.cpu_count() // 2}, - "file_io_sync": False, - }, - } - - -def get_zarr_write_spec( - file_path: str, chunk_size: int, base_shape: tuple, dtype: str -) -> dict: - """ - Returns a dictionary containing the specification for writing a Zarr file. - - Args: - file_path (str): The path to the Zarr file. - chunk_size (int): The size of the chunks in the Zarr file. - base_shape (tuple): The base shape of the Zarr file. - dtype (str): The data type of the Zarr file. - - Returns: - dict: A dictionary containing the specification for writing the Zarr file. - """ - return { - "driver": "zarr", - "kvstore": { - "driver": "file", - "path": file_path, - }, - "create": True, - "delete_existing": False, - "open": True, - "metadata": { - "shape": base_shape, - "chunks": [1, 1, 1, chunk_size, chunk_size], - "dtype": np.dtype(dtype).str, - "compressor": { - "id": "blosc", - "cname": "zstd", - "clevel": 1, - "shuffle": 1, - "blocksize": 0, - }, - }, - "context": { - "cache_pool": {}, - "data_copy_concurrency": {"limit": os.cpu_count() // 2}, - "file_io_concurrency": {"limit": os.cpu_count() // 2}, - "file_io_sync": False, - }, - } - - -class PyramidCompositorPY: - """ - A class for composing a group of pyramid images into an assembled pyramid structure. - """ - - def __init__( - self, input_pyramids_loc: str, out_dir: str, output_pyramid_name: str - ) -> None: - """ - Initializes the PyramidCompositor object. - - Args: - input_pyramids_loc (str): The location of the input pyramid images. - out_dir (str): The output directory for the composed zarr pyramid file. - output_pyramid_name (str): The name of the zarr pyramid file. - """ - self._input_pyramids_loc: str = input_pyramids_loc - self._chunk_cache: set = set() - self._output_pyramid_name: str = f"{out_dir}/{output_pyramid_name}" - self._ome_metadata_file: str = f"{out_dir}/{output_pyramid_name}/METADATA.ome.xml" - self._composition_map: dict = None - self._plate_image_shapes: dict = {} - self._zarr_arrays: dict = {} - self._unit_image_shapes: dict = {} - self._pyramid_levels: int = None - self._image_dtype: np.dtype = None - self._num_channels: int = None - - def _create_xml(self) -> None: - """ - Writes an OME-XML metadata file for the pyramid. - """ - ome_metadata = ome_types.model.OME() - ome_metadata.images.append( - ome_types.model.Image( - id="Image:0", - pixels=ome_types.model.Pixels( - id="Pixels:0", - dimension_order="XYZCT", - big_endian=False, - size_c=self._plate_image_shapes[0][1], - size_z=self._plate_image_shapes[0][2], - size_t=self._plate_image_shapes[0][0], - size_x=self._plate_image_shapes[0][4], - size_y=self._plate_image_shapes[0][3], - channels=[ - ome_types.model.Channel( - id=f"Channel:{i}", - samples_per_pixel=1, - ) - for i in range(self._plate_image_shapes[0][1]) - ], - type=OME_DTYPE[str(self._image_dtype)], - ), - ) - ) - - with open(self._ome_metadata_file, "w") as fw: - fw.write(str(ome_metadata.to_xml())) - - def _create_zattr_file(self) -> None: - """ - Creates a .zattrs file for the zarr pyramid. - """ - attr_dict: dict = {} - multiscale_metadata: list = [] - for key in self._plate_image_shapes: - multiscale_metadata.append({"path": str(key)}) - attr_dict["datasets"] = multiscale_metadata - attr_dict["version"] = "0.1" - attr_dict["name"] = self._output_pyramid_name - attr_dict["metadata"] = {"method": "mean"} - - final_attr_dict = {"multiscales": [attr_dict]} - - with open(f"{self._output_pyramid_name}/data.zarr/0/.zattrs", "w") as json_file: - json.dump(final_attr_dict, json_file) - - def _create_zgroup_file(self) -> None: - """ - Creates .zgroup files for the zarr pyramid. - """ - zgroup_dict = {"zarr_format": 2} - - with open(f"{self._output_pyramid_name}/data.zarr/0/.zgroup", "w") as json_file: - json.dump(zgroup_dict, json_file) - - with open(f"{self._output_pyramid_name}/data.zarr/.zgroup", "w") as json_file: - json.dump(zgroup_dict, json_file) - - def _create_auxilary_files(self) -> None: - """ - Creates auxiliary files (OME-XML, .zattrs, .zgroup) for the pyramid. - """ - self._create_xml() - self._create_zattr_file() - self._create_zgroup_file() - - def _write_zarr_chunk( - self, level: int, channel: int, y_index: int, x_index: int - ) -> None: - """ - Writes the chunk file at the specified level, channel, y_index, and x_index. - - Args: - level (int): The level of the pyramid. - channel (int): The channel of the pyramid. - y_index (int): The y-index of the tile. - x_index (int): The x-index of the tile. - """ - - # x_range and y_range are in global coordinates at the corresponding level - y_range = [ - y_index * CHUNK_SIZE, - min((y_index + 1) * CHUNK_SIZE, self._zarr_arrays[level].shape[3]), - ] - x_range = [ - x_index * CHUNK_SIZE, - min((x_index + 1) * CHUNK_SIZE, self._zarr_arrays[level].shape[4]), - ] - - assembled_width = x_range[1] - x_range[0] - assembled_height = y_range[1] - y_range[0] - assembled_image = np.zeros( - (assembled_height, assembled_width), dtype=self._image_dtype - ) - - # find what input images are needed - - # are we at the begining of the input image? - unit_image_height = self._unit_image_shapes[level][0] - unit_image_width = self._unit_image_shapes[level][1] - - row_start_pos = y_range[0] - while row_start_pos < y_range[1]: - # row and col are unit map coordinates - row = row_start_pos // unit_image_height - # local_y* and local_x* are coordintes in the chunk - local_y_start = row_start_pos - y_range[0] - tile_y_start = row_start_pos - row * unit_image_height - tile_y_dim = min((row + 1) * unit_image_height - row_start_pos, y_range[1] - row_start_pos) - tile_y_end = tile_y_start + tile_y_dim - col_start_pos = x_range[0] - while col_start_pos < x_range[1]: - col = col_start_pos // unit_image_width - local_x_start = col_start_pos - x_range[0] - tile_x_start = col_start_pos - col * unit_image_width - tile_x_dim = min((col + 1) * unit_image_width - col_start_pos, x_range[1] - col_start_pos) - tile_x_end = tile_x_start + tile_x_dim - # read input zarr file - input_file_name = self._composition_map.get((col, row, channel)) - zarr_file_loc = Path(input_file_name) / "data.zarr/0/" - zarr_array_loc = zarr_file_loc / str(level) - zarr_file = ts.open(get_zarr_read_spec(str(zarr_array_loc))).result() - - # copy data - assembled_image[ - local_y_start : local_y_start + tile_y_end - tile_y_start, - local_x_start : local_x_start + tile_x_end - tile_x_start, - ] = ( - zarr_file[0, 0, 0, tile_y_start:tile_y_end, tile_x_start:tile_x_end] - .read() - .result() - ) - col_start_pos += tile_x_end - tile_x_start # update col index - - row_start_pos += tile_y_end - tile_y_start - - zarr_array = self._zarr_arrays[level] - zarr_array[ - 0, channel, 0, y_range[0] : y_range[1], x_range[0] : x_range[1] - ].write(assembled_image).result() - - def set_composition(self, composition_map: dict) -> None: - """ - Sets the composition for the pyramid. - - Args: - composition_map (dict): A dictionary mapping composition images to file paths. - """ - self._composition_map = composition_map - self._unit_image_shapes = {} - for coord in composition_map: - file = composition_map[coord] - zarr_file_loc = Path(file) / "data.zarr/0/" - attr_file_loc = Path(file) / "data.zarr/0/.zattrs" - if attr_file_loc.exists(): - with open(str(attr_file_loc), "r") as f: - attrs = json.load(f) - mutliscale_metadata = attrs["multiscales"][0]["datasets"] - self._pyramid_levels = len(mutliscale_metadata) - for dic in mutliscale_metadata: - res_key = dic["path"] - zarr_array_loc = zarr_file_loc / res_key - zarr_file = ts.open( - get_zarr_read_spec(str(zarr_array_loc)) - ).result() - self._unit_image_shapes[int(res_key)] = ( - zarr_file.shape[-2], - zarr_file.shape[-1], - ) - if res_key == "0": - self._image_dtype = zarr_file.dtype.numpy_dtype - - break - - num_rows = 0 - num_cols = 0 - num_channels = 0 - - for coord in composition_map: - num_rows = max(num_rows, coord[1]) - num_cols = max(num_cols, coord[0]) - num_channels = max(num_channels, coord[2]) - - num_cols += 1 - num_rows += 1 - num_channels += 1 - - self._num_channels = num_channels - - self._plate_image_shapes = {} - self._zarr_arrays = {} - self._chunk_cache = set() - for l in self._unit_image_shapes: - level = int(l) - self._plate_image_shapes[level] = (1, - 1, - num_channels, - 1, - num_rows * self._unit_image_shapes[level][0], - num_cols * self._unit_image_shapes[level][1], - ) - num_row_tiles = math.ceil( - 1.0 * num_rows * self._unit_image_shapes[level][0] / CHUNK_SIZE - ) - num_col_tiles = math.ceil( - 1.0 * num_cols * self._unit_image_shapes[level][1] / CHUNK_SIZE - ) - if num_row_tiles == 0: - num_row_tiles = 1 - if num_col_tiles == 0: - num_col_tiles == 1 - self._zarr_arrays[level] = ts.open( - get_zarr_write_spec( - f"{self._output_pyramid_name}/data.zarr/0/{level}", - CHUNK_SIZE, - self._plate_image_shapes[level], - np.dtype(self._image_dtype).str, - ) - ).result() - - self._create_auxilary_files() - - def reset_composition(self) -> None: - """ - Resets the pyramid composition by removing the pyramid file and clearing internal data structures. - """ - shutil.rmtree(self._output_pyramid_name) - self._composition_map = None - self._plate_image_shapes = None - self._chunk_cache = None - self._plate_image_shapes = {} - self._zarr_arrays = {} - - def get_zarr_chunk( - self, level: int, channel: int, y_index: int, x_index: int - ) -> None: - """ - Retrieves zarr chunk data from the pyramid at the specified level, channel, y_index, and x_index. - - This method checks if the requested zarr chunk is already in the cache. If not, it calls - the `_write_zarr_chunk` method to generate the chunk and add it to the cache. - - Args: - level (int): The level of the pyramid. - channel (int): The channel of the pyramid. - y_index (int): The y-index of the tile. - x_index (int): The x-index of the tile. - - Raises: - ValueError: If the composition map is not set, the requested level does not exist, - the requested channel does not exist, or the requested y_index or x_index - is out of bounds. - """ - if self._composition_map is None: - raise ValueError("No composition map is set. Unable to generate pyramid") - - if level not in self._unit_image_shapes: - raise ValueError(f"Requested level ({level}) does not exist") - - if channel >= self._num_channels: - raise ValueError(f"Requested channel ({channel}) does not exist") - - if y_index > (self._plate_image_shapes[level][3] // CHUNK_SIZE): - raise ValueError(f"Requested y index ({y_index}) does not exist") - - if x_index > (self._plate_image_shapes[level][4] // CHUNK_SIZE): - raise ValueError(f"Requested y index ({x_index}) does not exist") - - if (level, channel, y_index, x_index) in self._chunk_cache: - return - else: - self._write_zarr_chunk(level, channel, y_index, x_index) - self._chunk_cache.add((level, channel, y_index, x_index)) - return - -class PyramidCompositor2: +class PyramidCompositor: """ A class for composing a group of pyramid images into an assembled pyramid structure. """ @@ -419,7 +28,6 @@ def __init__( output_pyramid_name (str): The name of the zarr pyramid file. """ self._pyramid_compositor_cpp = PyramidCompositorCPP(input_pyramids_loc, out_dir, output_pyramid_name) - self._chunk_cache: set = set() def _create_xml(self) -> None: """ From cb47a72ca6efe84654301887fd6042fcf06b15bd Mon Sep 17 00:00:00 2001 From: JesseMckinzie Date: Wed, 22 Jan 2025 10:21:39 -0700 Subject: [PATCH 09/24] Remove outdated files --- .../core/pyramid_compositor_from_trash.cpp | 514 ------------------ src/cpp/core/pyramid_compositor_old.cpp | 514 ------------------ 2 files changed, 1028 deletions(-) delete mode 100644 src/cpp/core/pyramid_compositor_from_trash.cpp delete mode 100644 src/cpp/core/pyramid_compositor_old.cpp diff --git a/src/cpp/core/pyramid_compositor_from_trash.cpp b/src/cpp/core/pyramid_compositor_from_trash.cpp deleted file mode 100644 index 444db5a..0000000 --- a/src/cpp/core/pyramid_compositor_from_trash.cpp +++ /dev/null @@ -1,514 +0,0 @@ -#include "pyramid_compositor.h" -#include "../utilities/utilities.h" - -#include -#include -#include -#include - -#include "tensorstore/context.h" -#include "tensorstore/array.h" -#include "tensorstore/driver/zarr/dtype.h" -#include "tensorstore/index_space/dim_expression.h" -#include "tensorstore/kvstore/kvstore.h" -#include "tensorstore/open.h" - -#include - -using ::tensorstore::internal_zarr::ChooseBaseDType; -using json = nlohmann::json; - -namespace fs = std::filesystem; - -namespace argolid { -PyramidCompositor::PyramidCompositor(const std::string& input_pyramids_loc, const std::string& out_dir, const std::string& output_pyramid_name): - _input_pyramids_loc(input_pyramids_loc), - _output_pyramid_name(out_dir + "/" + output_pyramid_name), - _ome_metadata_file(out_dir + "/" + output_pyramid_name + "/METADATA.ome.xml"), - _pyramid_levels(-1), - _num_channels(-1) { - - - tensorstore::TensorStore source; - - TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( - GetZarrSpecToRead(input_pyramids_loc), - tensorstore::OpenMode::open, - tensorstore::ReadWriteMode::read).result()); - - _image_dtype = source.dtype().name(); - _data_type_code = GetDataTypeCode(_image_dtype); - - // auto read_spec = GetZarrSpecToRead(input_pyramids_loc) - - // TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( - // read_spec, - // tensorstore::OpenMode::open, - // tensorstore::ReadWriteMode::read).result()); - - // auto image_shape = source.domain().shape(); - // const auto read_chunk_shape = source.chunk_layout().value().read_chunk_shape(); - - // if (image_shape.size() == 5){ - // _image_height = image_shape[3]; - // _image_width = image_shape[4]; - // _image_depth = image_shape[2]; - // _num_channels = image_shape[1]; - // _num_tsteps = image_shape[0]; - // _t_int.emplace(0); - // _c_int.emplace(1); - // _z_int.emplace(2); - // } else { - // assert(image_shape.size() >= 2); - // std::tie(_t_int, _c_int, _z_int) = ParseMultiscaleMetadata(axes_list, image_shape.size()); - // _image_height = image_shape[image_shape.size()-2]; - // _image_width = image_shape[image_shape.size()-1]; - // if (_t_int.has_value()) { - // _num_tsteps = image_shape[_t_int.value()]; - // } else { - // _num_tsteps = 1; - // } - - // if (_c_int.has_value()) { - // _num_channels = image_shape[_c_int.value()]; - // } else { - // _num_channels = 1; - // } - - // if (_z_int.has_value()) { - // _image_depth = image_shape[_z_int.value()]; - // } else { - // _image_depth = 1; - // } - // } - - // _data_type = source.dtype().name(); - // _data_type_code = GetDataTypeCode(_data_type); - - // TENSORSTORE_CHECK_OK_AND_ASSIGN(___mb_cur_max, tensorstore::Open( - // GetZarrSpecToWrite(_filename, _image_shape, _chunk_shape, GetEncodedType(_dtype_code)), - // tensorstore::OpenMode::create | - // tensorstore::OpenMode::delete_existing, - // tensorstore::ReadWriteMode::write).result() - // ); -} - -// Private Methods -void PyramidCompositor::create_xml() { - CreateXML( - this->_output_pyramid_name, - this->_ome_metadata_file, - this->_plate_image_shapes[0], - this->_image_dtype - ); -} - -void PyramidCompositor::create_zattr_file() { - WriteTSZattrFilePlateImage( - this->_output_pyramid_name, - this->_output_pyramid_name + "/data.zarr/0", - this->_plate_image_shapes - ); -} - -void PyramidCompositor::create_zgroup_file() { - WriteVivZgroupFiles(this->_output_pyramid_name); -} - -void PyramidCompositor::create_auxiliary_files() { - create_xml(); - create_zattr_file(); - create_zgroup_file(); -} - -void PyramidCompositor::write_zarr_chunk(int level, int channel, int y_int, int x_int) { - // Implementation of tile data writing logic - // Placeholder logic - - // std::tuple y_range = { - // y_int * CHUNK_SIZE, - // std::min((y_int+1) * CHUNK_SIZE, this->zarr_arrays_[level].size()[3]) - // }; - - // std::tuple x_range = { - // x_int * CHUNK_SIZE, - // min((x_int+1) * CHUNK_SIZE, this->zarr_arrays_[level].size()[4]) - // }; - - // int assembled_width = std::get<1>(x_range) - std::get<0>(x_range); - // int assembled_height = std::get<1>(y_range) - std::get<0>(y_range); - - // std::vector<> - - // Compute ranges in global coordinates - int y_start = y_int * CHUNK_SIZE; - int y_end = std::min((y_int + 1) * CHUNK_SIZE, std::get<1>(_zarr_arrays[level])[0]); - int x_start = x_int * CHUNK_SIZE; - int x_end = std::min((x_int + 1) * CHUNK_SIZE, std::get<1>(_zarr_arrays[level])[1]); - - int assembled_width = x_end - x_start; - int assembled_height = y_end - y_start; - - image_data assembled_image = GetAssembledImageVector(assembled_height*assembled_width); - - // Determine required input images - int unit_image_height = _unit_image_shapes[level].first; - int unit_image_width = _unit_image_shapes[level].second; - - int row_start_pos = y_start; - while (row_start_pos < y_end) { - int row = row_start_pos / unit_image_height; - int local_y_start = row_start_pos - y_start; - int tile_y_start = row_start_pos - row * unit_image_height; - int tile_y_dim = std::min((row + 1) * unit_image_height - row_start_pos, y_end - row_start_pos); - int tile_y_end = tile_y_start + tile_y_dim; - - int col_start_pos = x_start; - while (col_start_pos < x_end) { - int col = col_start_pos / unit_image_width; - int local_x_start = col_start_pos - x_start; - int tile_x_start = col_start_pos - col * unit_image_width; - int tile_x_dim = std::min((col + 1) * unit_image_width - col_start_pos, x_end - col_start_pos); - int tile_x_end = tile_x_start + tile_x_dim; - - // Fetch input file path from composition map - std::string input_file_name = _composition_map[{col, row, channel}]; - std::filesystem::path zarrArrayLoc = std::filesystem::path(input_file_name) / "data.zarr/0/" / std::to_string(level); - - auto read_data = std::get<0>(ReadImageData(zarrArrayLoc, Seq(tile_y_start, tile_y_end), Seq(tile_x_start, tile_x_end))).get(); - - for (int i = local_y_start; i < local_y_start + tile_y_end - tile_y_start; ++i) { - for (int j = local_x_start; i < local_x_start + tile_x_end - tile_x_start; ++j) { - assembled_image[i*(tile_x_end - tile_x_start) + j] = read_data[(i-local_y_start)*(tile_x_end - tile_x_start) + (j-local_x_start)]; - } - } - - // // Open the Zarr array for reading with TensorStore - // auto zarrSpec = tensorstore::Context::Default() - // | tensorstore::OpenMode::read - // | tensorstore::dtype_v - // | tensorstore::Shape({1, 1, 1, assembledHeight, assembledWidth}); - - // auto inputZarrArray = tensorstore::Open(zarrSpec).result(); - - // TENSORSTORE_CHECK_OK_AND_ASSIGN() - - // if (!inputZarrArray.ok()) { - // std::cerr << "Failed to open Zarr array: " << inputZarrArray.status() << "\n"; - // return; - // } - - // // Read data from input Zarr file into assembledImage - // tensorstore::Read(*inputZarrArray, - // tensorstore::Span(assembledImage.data(), assembledImage.size())).result(); - - col_start_pos += tile_x_end - tile_x_start; - } - row_start_pos += tile_y_end - tile_y_start; - } - - - WriteImageData( - _input_pyramids_loc, - std::get<0>(_zarr_arrays[level]).get(), - Seq(y_start,y_end), Seq(x_start, x_end), - Seq(0,0), Seq(channel, channel), Seq(0,0) - ); -} - -void PyramidCompositor::setComposition(const std::unordered_map, std::string>& comp_map) { - _composition_map = comp_map; - - for (const auto& coord : comp_map) { - std::string file_path = coord.second; - std::filesystem::path attr_file_loc = std::filesystem::path(file_path) / "data.zarr/0/.zattrs"; - - if (std::filesystem::exists(attr_file_loc)) { - std::ifstream attr_file(attr_file_loc); - json attrs; - attr_file >> attrs; - - const auto& multiscale_metadata = attrs["multiscales"][0]["datasets"]; - _pyramid_levels = multiscale_metadata.size(); - for (const auto& dic : multiscale_metadata) { - std::string res_key = dic["path"]; - std::filesystem::path zarr_array_loc = std::filesystem::path(file_path) / "data.zarr/0/" / res_key; - - auto read_spec = GetZarrSpecToRead(zarr_array_loc); - - tensorstore::TensorStore source; - - TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( - read_spec, - tensorstore::OpenMode::open, - tensorstore::ReadWriteMode::read).result()); - - auto image_shape = source.domain().shape(); - _image_dtype = source.dtype().name(); - _data_type_code = GetDataTypeCode(_image_dtype); - - _unit_image_shapes[std::stoi(res_key)] = { - image_shape[image_shape.size()-2], - image_shape[image_shape.size()-1] - }; - } - break; - } - } - - int num_rows = 0, num_cols = 0, num_channels = 0; - for (const auto& coord : comp_map) { - num_rows = std::max(num_rows, std::get<1>(coord.first)); - num_cols = std::max(num_cols, std::get<0>(coord.first)); - _num_channels = std::max(num_channels, std::get<2>(coord.first)); - } - num_cols += 1; - num_rows += 1; - num_channels += 1; - _num_channels = num_channels; - - for (const auto& [level, shape] : _unit_image_shapes) { - std::string path = _output_pyramid_name + "/data.zarr/0/" + std::to_string(level); - - auto zarr_array_data = ReadImageData( - path, - Seq(0, num_rows * shape.first), - Seq(0, num_cols * shape.second), - Seq(0,0), - Seq(0, num_channels), - Seq(0,0) - ); - - _zarr_arrays[level] = zarr_array_data; - } -} - -void PyramidCompositor::reset_composition() { - // Remove pyramid file and clear internal data structures - fs::remove_all(_output_pyramid_name); - _composition_map.clear(); - _plate_image_shapes.clear(); - _tile_cache.clear(); - _plate_image_shapes.clear(); - _zarr_arrays.clear(); -} - -image_data PyramidCompositor::GetAssembledImageVector(int size){ - switch (_image_dtype_code) - { - case (1): - return std::vector(size); - case (2): - return std::vector(size); - case (4): - return std::vector(size); - case (8): - return std::vector(size); - case (16): - return std::vector(size); - case (32): - return std::vector(size); - case (64): - return std::vector(size); - case (128): - return std::vector(size); - case (256): - return std::vector(size); - case (512): - return std::vector(size); - default: - throw std::runtime_error("Invalid data type."); - } -} - -void PyramidCompositor::WriteImageData( - const std::string& path, - image_data& image, - const Seq& rows, - const Seq& cols, - const std::optional& layers, - const std::optional& channels, - const std::optional& tsteps) { - - std::vector shape; - - tensorstore::TensorStore source; - - TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( - GetZarrSpecToRead(path), - tensorstore::OpenMode::open, - tensorstore::ReadWriteMode::write).result()); - - auto output_transform = tensorstore::IdentityTransform(source.domain()); - - if (_t_index.has_value() && tsteps.has_value()) { - output_transform = (std::move(output_transform) | tensorstore::Dims(_t_int.value()).ClosedInterval(tsteps.value().Start(), tsteps.value().Stop())).value(); - shape.emplace_back(tsteps.value().Stop() - tsteps.value().Start()+1); - } - - if (_c_int.has_value() && channels.has_value()) { - output_transform = (std::move(output_transform) | tensorstore::Dims(_c_int.value()).ClosedInterval(channels.value().Start(), channels.value().Stop())).value(); - shape.emplace_back(channels.value().Stop() - channels.value().Start()+1); - } - - if (_z_int.has_value() && layers.has_value()) { - output_transform = (std::move(output_transform) | tensorstore::Dims(_z_int.value()).ClosedInterval(layers.value().Start(), layers.value().Stop())).value(); - shape.emplace_back(layers.value().Stop() - layers.value().Start()+1); - } - - output_transform = (std::move(output_transform) | tensorstore::Dims(_y_int).ClosedInterval(rows.Start(), rows.Stop()) | - tensorstore::Dims(_x_int).ClosedInterval(cols.Start(), cols.Stop())).value(); - - shape.emplace_back(rows.Stop() - rows.Start()+1); - shape.emplace_back(cols.Stop() - cols.Start()+1); - - auto data_array = tensorstore::Array(image.data(), shape, tensorstore::c_order); - - // Write data array to TensorStore - auto write_result = tensorstore::Write(tensorstore::UnownedToShared(data_array), source | output_transform).result(); - - if (!write_result.ok()) { - std::cerr << "Error writing image: " << write_result.status() << std::endl; - } -} - -template -std::shared_ptr> PyramidCompositor::GetImageDataTemplated(const Seq& rows, const Seq& cols, const Seq& layers, const Seq& channels, const Seq& tsteps){ - - const auto data_height = rows.Stop() - rows.Start() + 1; - const auto data_width = cols.Stop() - cols.Start() + 1; - const auto data_depth = layers.Stop() - layers.Start() + 1; - const auto data_num_channels = channels.Stop() - channels.Start() + 1; - const auto data_tsteps = tsteps.Stop() - tsteps.Start() + 1; - - auto read_buffer = std::make_shared>(data_height*data_width*data_depth*data_num_channels*data_tsteps); - auto array = tensorstore::Array(read_buffer->data(), {data_tsteps, data_num_channels, data_depth, data_height, data_width}, tensorstore::c_order); - tensorstore::intTransform<> read_transform = tensorstore::IdentityTransform(source.domain()); - - int x_int=1, y_int=0; - if (_t_int.has_value()){ - read_transform = (std::move(read_transform) | tensorstore::Dims(_t_int.value()).ClosedInterval(tsteps.Start(), tsteps.Stop())).value(); - x_int++; - y_int++; - } - if (_c_int.has_value()){ - read_transform = (std::move(read_transform) | tensorstore::Dims(_c_int.value()).ClosedInterval(channels.Start(), channels.Stop())).value(); - x_int++; - y_int++; - } - if (_z_int.has_value()){ - read_transform = (std::move(read_transform) | tensorstore::Dims(_z_int.value()).ClosedInterval(layers.Start(), layers.Stop())).value(); - x_int++; - y_int++; - } - read_transform = (std::move(read_transform) | tensorstore::Dims(y_int).ClosedInterval(rows.Start(), rows.Stop()) | - tensorstore::Dims(x_int).ClosedInterval(cols.Start(), cols.Stop())).value(); - - - tensorstore::Read(source | read_transform, tensorstore::UnownedToShared(array)).value(); - - return read_buffer; -} - -std::tuple, std::vector> PyramidCompositor::ReadImageData( - const std::string& fname, - const Seq& rows, const Seq& cols, - const Seq& layers = Seq(0,0), - const Seq& channels = Seq(0,0), - const Seq& tsteps = Seq(0,0)) { - - auto read_spec = GetZarrSpecToRead(fname); - - tensorstore::TensorStore source; - - TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( - read_spec, - tensorstore::OpenMode::open, - tensorstore::ReadWriteMode::read).result()); - - auto image_shape = source.domain().shape(); - const auto read_chunk_shape = source.chunk_layout().value().read_chunk_shape(); - - // if (image_shape.size() == 5){ - // _image_height = image_shape[3]; - // _image_width = image_shape[4]; - // _image_depth = image_shape[2]; - // _num_channels = image_shape[1]; - // _num_tsteps = image_shape[0]; - // _t_index.emplace(0); - // _c_index.emplace(1); - // _z_index.emplace(2); - // } else { - // assert(image_shape.size() >= 2); - // std::tie(_t_index, _c_index, _z_index) = ParseMultiscaleMetadata(axes_list, image_shape.size()); - // _image_height = image_shape[image_shape.size()-2]; - // _image_width = image_shape[image_shape.size()-1]; - // if (_t_index.has_value()) { - // _num_tsteps = image_shape[_t_index.value()]; - // } else { - // _num_tsteps = 1; - // } - - // if (_c_index.has_value()) { - // _num_channels = image_shape[_c_index.value()]; - // } else { - // _num_channels = 1; - // } - - // if (_z_index.has_value()) { - // _image_depth = image_shape[_z_index.value()]; - // } else { - // _image_depth = 1; - // } - // } - - _image_dtype = source.dtype().name(); - _data_type_code = GetDataTypeCode(_image_dtype); - - switch (_data_type_code) - { - case (1): - return std::make_tuple( - std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), - std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); - case (2): - return std::make_tuple( - std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), - std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); - case (4): - return std::make_tuple( - std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), - std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); - case (8): - return std::make_tuple( - std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), - std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); - case (16): - return std::make_tuple( - std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), - std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); - case (32): - return std::make_tuple( - std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), - std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); - case (64): - return std::make_tuple( - std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), - std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); - case (128): - return std::make_tuple( - std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), - std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); - case (256): - return std::make_tuple( - std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), - std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); - case (512): - return std::make_tuple( - std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), - std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); - default: - break; - } -} -} // ns argolid \ No newline at end of file diff --git a/src/cpp/core/pyramid_compositor_old.cpp b/src/cpp/core/pyramid_compositor_old.cpp deleted file mode 100644 index 40fa5ff..0000000 --- a/src/cpp/core/pyramid_compositor_old.cpp +++ /dev/null @@ -1,514 +0,0 @@ -#include "pyramid_compositor.h" -#include "../utilities/utilities.h" - -#include -#include -#include -#include - -#include "tensorstore/context.h" -#include "tensorstore/array.h" -#include "tensorstore/driver/zarr/dtype.h" -#include "tensorstore/index_space/dim_expression.h" -#include "tensorstore/kvstore/kvstore.h" -#include "tensorstore/open.h" - -#include - -using ::tensorstore::internal_zarr::ChooseBaseDType; -using json = nlohmann::json; - -namespace fs = std::filesystem; - -namespace argolid { -PyramidCompositor::PyramidCompositor(const std::string& input_pyramids_loc, const std::string& out_dir, const std::string& output_pyramid_name): - _input_pyramids_loc(input_pyramids_loc), - _output_pyramid_name(out_dir + "/" + output_pyramid_name), - _ome_metadata_file(out_dir + "/" + output_pyramid_name + "/METADATA.ome.xml"), - _pyramid_levels(-1), - _num_channels(-1) { - - - tensorstore::TensorStore source; - - TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( - GetZarrSpecToRead(input_pyramids_loc), - tensorstore::OpenMode::open, - tensorstore::ReadWriteMode::read).result()); - - _image_dtype = source.dtype().name(); - _data_type_code = GetDataTypeCode(_image_dtype); - - // auto read_spec = GetZarrSpecToRead(input_pyramids_loc) - - // TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( - // read_spec, - // tensorstore::OpenMode::open, - // tensorstore::ReadWriteMode::read).result()); - - // auto image_shape = source.domain().shape(); - // const auto read_chunk_shape = source.chunk_layout().value().read_chunk_shape(); - - // if (image_shape.size() == 5){ - // _image_height = image_shape[3]; - // _image_width = image_shape[4]; - // _image_depth = image_shape[2]; - // _num_channels = image_shape[1]; - // _num_tsteps = image_shape[0]; - // _t_index.emplace(0); - // _c_index.emplace(1); - // _z_index.emplace(2); - // } else { - // assert(image_shape.size() >= 2); - // std::tie(_t_index, _c_index, _z_index) = ParseMultiscaleMetadata(axes_list, image_shape.size()); - // _image_height = image_shape[image_shape.size()-2]; - // _image_width = image_shape[image_shape.size()-1]; - // if (_t_index.has_value()) { - // _num_tsteps = image_shape[_t_index.value()]; - // } else { - // _num_tsteps = 1; - // } - - // if (_c_index.has_value()) { - // _num_channels = image_shape[_c_index.value()]; - // } else { - // _num_channels = 1; - // } - - // if (_z_index.has_value()) { - // _image_depth = image_shape[_z_index.value()]; - // } else { - // _image_depth = 1; - // } - // } - - // _data_type = source.dtype().name(); - // _data_type_code = GetDataTypeCode(_data_type); - - // TENSORSTORE_CHECK_OK_AND_ASSIGN(___mb_cur_max, tensorstore::Open( - // GetZarrSpecToWrite(_filename, _image_shape, _chunk_shape, GetEncodedType(_dtype_code)), - // tensorstore::OpenMode::create | - // tensorstore::OpenMode::delete_existing, - // tensorstore::ReadWriteMode::write).result() - // ); -} - -// Private Methods -void PyramidCompositor::create_xml() { - CreateXML( - this->_output_pyramid_name, - this->_ome_metadata_file, - this->_plate_image_shapes[0], - this->_image_dtype - ); -} - -void PyramidCompositor::create_zattr_file() { - WriteTSZattrFilePlateImage( - this->_output_pyramid_name, - this->_output_pyramid_name + "/data.zarr/0", - this->_plate_image_shapes - ); -} - -void PyramidCompositor::create_zgroup_file() { - WriteVivZgroupFiles(this->_output_pyramid_name); -} - -void PyramidCompositor::create_auxiliary_files() { - create_xml(); - create_zattr_file(); - create_zgroup_file(); -} - -void PyramidCompositor::write_zarr_chunk(int level, int channel, int y_int, int x_int) { - // Implementation of tile data writing logic - // Placeholder logic - - // std::tuple y_range = { - // y_int * CHUNK_SIZE, - // std::min((y_int+1) * CHUNK_SIZE, this->zarr_arrays_[level].size()[3]) - // }; - - // std::tuple x_range = { - // x_int * CHUNK_SIZE, - // min((x_int+1) * CHUNK_SIZE, this->zarr_arrays_[level].size()[4]) - // }; - - // int assembled_width = std::get<1>(x_range) - std::get<0>(x_range); - // int assembled_height = std::get<1>(y_range) - std::get<0>(y_range); - - // std::vector<> - - // Compute ranges in global coordinates - int y_start = y_int * CHUNK_SIZE; - int y_end = std::min((y_int + 1) * CHUNK_SIZE, std::get<1>(_zarr_arrays[level])[0]); - int x_start = x_int * CHUNK_SIZE; - int x_end = std::min((x_int + 1) * CHUNK_SIZE, std::get<1>(_zarr_arrays[level])[1]); - - int assembled_width = x_end - x_start; - int assembled_height = y_end - y_start; - - image_data assembled_image = GetAssembledImageVector(assembled_height*assembled_width); - - // Determine required input images - int unit_image_height = _unit_image_shapes[level].first; - int unit_image_width = _unit_image_shapes[level].second; - - int row_start_pos = y_start; - while (row_start_pos < y_end) { - int row = row_start_pos / unit_image_height; - int local_y_start = row_start_pos - y_start; - int tile_y_start = row_start_pos - row * unit_image_height; - int tile_y_dim = std::min((row + 1) * unit_image_height - row_start_pos, y_end - row_start_pos); - int tile_y_end = tile_y_start + tile_y_dim; - - int col_start_pos = x_start; - while (col_start_pos < x_end) { - int col = col_start_pos / unit_image_width; - int local_x_start = col_start_pos - x_start; - int tile_x_start = col_start_pos - col * unit_image_width; - int tile_x_dim = std::min((col + 1) * unit_image_width - col_start_pos, x_end - col_start_pos); - int tile_x_end = tile_x_start + tile_x_dim; - - // Fetch input file path from composition map - std::string input_file_name = _composition_map[{col, row, channel}]; - std::filesystem::path zarrArrayLoc = std::filesystem::path(input_file_name) / "data.zarr/0/" / std::to_string(level); - - auto read_data = std::get<0>(ReadImageData(zarrArrayLoc, Seq(tile_y_start, tile_y_end), Seq(tile_x_start, tile_x_end))).get(); - - for (int i = local_y_start; i < local_y_start + tile_y_end - tile_y_start; ++i) { - for (int j = local_x_start; i < local_x_start + tile_x_end - tile_x_start; ++j) { - assembled_image[i*(tile_x_end - tile_x_start) + j] = read_data[(i-local_y_start)*(tile_x_end - tile_x_start) + (j-local_x_start)]; - } - } - - // // Open the Zarr array for reading with TensorStore - // auto zarrSpec = tensorstore::Context::Default() - // | tensorstore::OpenMode::read - // | tensorstore::dtype_v - // | tensorstore::Shape({1, 1, 1, assembledHeight, assembledWidth}); - - // auto inputZarrArray = tensorstore::Open(zarrSpec).result(); - - // TENSORSTORE_CHECK_OK_AND_ASSIGN() - - // if (!inputZarrArray.ok()) { - // std::cerr << "Failed to open Zarr array: " << inputZarrArray.status() << "\n"; - // return; - // } - - // // Read data from input Zarr file into assembledImage - // tensorstore::Read(*inputZarrArray, - // tensorstore::Span(assembledImage.data(), assembledImage.size())).result(); - - col_start_pos += tile_x_end - tile_x_start; - } - row_start_pos += tile_y_end - tile_y_start; - } - - - WriteImageData( - _input_pyramids_loc, - std::get<0>(_zarr_arrays[level]).get(), - Seq(y_start,y_end), Seq(x_start, x_end), - Seq(0,0), Seq(channel, channel), Seq(0,0) - ); -} - -void PyramidCompositor::setComposition(const std::unordered_map, std::string>& comp_map) { - _composition_map = comp_map; - - for (const auto& coord : comp_map) { - std::string file_path = coord.second; - std::filesystem::path attr_file_loc = std::filesystem::path(file_path) / "data.zarr/0/.zattrs"; - - if (std::filesystem::exists(attr_file_loc)) { - std::ifstream attr_file(attr_file_loc); - json attrs; - attr_file >> attrs; - - const auto& multiscale_metadata = attrs["multiscales"][0]["datasets"]; - _pyramid_levels = multiscale_metadata.size(); - for (const auto& dic : multiscale_metadata) { - std::string res_key = dic["path"]; - std::filesystem::path zarr_array_loc = std::filesystem::path(file_path) / "data.zarr/0/" / res_key; - - auto read_spec = GetZarrSpecToRead(zarr_array_loc); - - tensorstore::TensorStore source; - - TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( - read_spec, - tensorstore::OpenMode::open, - tensorstore::ReadWriteMode::read).result()); - - auto image_shape = source.domain().shape(); - _image_dtype = source.dtype().name(); - _data_type_code = GetDataTypeCode(_image_dtype); - - _unit_image_shapes[std::stoi(res_key)] = { - image_shape[image_shape.size()-2], - image_shape[image_shape.size()-1] - }; - } - break; - } - } - - int num_rows = 0, num_cols = 0, num_channels = 0; - for (const auto& coord : comp_map) { - num_rows = std::max(num_rows, std::get<1>(coord.first)); - num_cols = std::max(num_cols, std::get<0>(coord.first)); - _num_channels = std::max(num_channels, std::get<2>(coord.first)); - } - num_cols += 1; - num_rows += 1; - num_channels += 1; - _num_channels = num_channels; - - for (const auto& [level, shape] : _unit_image_shapes) { - std::string path = _output_pyramid_name + "/data.zarr/0/" + std::to_string(level); - - auto zarr_array_data = ReadImageData( - path, - Seq(0, num_rows * shape.first), - Seq(0, num_cols * shape.second), - Seq(0,0), - Seq(0, num_channels), - Seq(0,0) - ); - - _zarr_arrays[level] = zarr_array_data; - } -} - -void PyramidCompositor::reset_composition() { - // Remove pyramid file and clear internal data structures - fs::remove_all(_output_pyramid_name); - _composition_map.clear(); - _plate_image_shapes.clear(); - _tile_cache.clear(); - _plate_image_shapes.clear(); - _zarr_arrays.clear(); -} - -image_data PyramidCompositor::GetAssembledImageVector(int size){ - switch (_image_dtype_code) - { - case (1): - return std::vector(size); - case (2): - return std::vector(size); - case (4): - return std::vector(size); - case (8): - return std::vector(size); - case (16): - return std::vector(size); - case (32): - return std::vector(size); - case (64): - return std::vector(size); - case (128): - return std::vector(size); - case (256): - return std::vector(size); - case (512): - return std::vector(size); - default: - throw std::runtime_error("Invalid data type."); - } -} - -void PyramidCompositor::WriteImageData( - const std::string& path, - image_data& image, - const Seq& rows, - const Seq& cols, - const std::optional& layers, - const std::optional& channels, - const std::optional& tsteps) { - - std::vector shape; - - tensorstore::TensorStore source; - - TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( - GetZarrSpecToRead(path), - tensorstore::OpenMode::open, - tensorstore::ReadWriteMode::write).result()); - - auto output_transform = tensorstore::IdentityTransform(source.domain()); - - if (_t_index.has_value() && tsteps.has_value()) { - output_transform = (std::move(output_transform) | tensorstore::Dims(_t_index.value()).ClosedInterval(tsteps.value().Start(), tsteps.value().Stop())).value(); - shape.emplace_back(tsteps.value().Stop() - tsteps.value().Start()+1); - } - - if (_c_index.has_value() && channels.has_value()) { - output_transform = (std::move(output_transform) | tensorstore::Dims(_c_index.value()).ClosedInterval(channels.value().Start(), channels.value().Stop())).value(); - shape.emplace_back(channels.value().Stop() - channels.value().Start()+1); - } - - if (_z_index.has_value() && layers.has_value()) { - output_transform = (std::move(output_transform) | tensorstore::Dims(_z_index.value()).ClosedInterval(layers.value().Start(), layers.value().Stop())).value(); - shape.emplace_back(layers.value().Stop() - layers.value().Start()+1); - } - - output_transform = (std::move(output_transform) | tensorstore::Dims(_y_index).ClosedInterval(rows.Start(), rows.Stop()) | - tensorstore::Dims(_x_index).ClosedInterval(cols.Start(), cols.Stop())).value(); - - shape.emplace_back(rows.Stop() - rows.Start()+1); - shape.emplace_back(cols.Stop() - cols.Start()+1); - - auto data_array = tensorstore::Array(image.data(), shape, tensorstore::c_order); - - // Write data array to TensorStore - auto write_result = tensorstore::Write(tensorstore::UnownedToShared(data_array), source | output_transform).result(); - - if (!write_result.ok()) { - std::cerr << "Error writing image: " << write_result.status() << std::endl; - } -} - -template -std::shared_ptr> PyramidCompositor::GetImageDataTemplated(const Seq& rows, const Seq& cols, const Seq& layers, const Seq& channels, const Seq& tsteps){ - - const auto data_height = rows.Stop() - rows.Start() + 1; - const auto data_width = cols.Stop() - cols.Start() + 1; - const auto data_depth = layers.Stop() - layers.Start() + 1; - const auto data_num_channels = channels.Stop() - channels.Start() + 1; - const auto data_tsteps = tsteps.Stop() - tsteps.Start() + 1; - - auto read_buffer = std::make_shared>(data_height*data_width*data_depth*data_num_channels*data_tsteps); - auto array = tensorstore::Array(read_buffer->data(), {data_tsteps, data_num_channels, data_depth, data_height, data_width}, tensorstore::c_order); - tensorstore::intTransform<> read_transform = tensorstore::IdentityTransform(source.domain()); - - int x_int=1, y_int=0; - if (_t_index.has_value()){ - read_transform = (std::move(read_transform) | tensorstore::Dims(_t_index.value()).ClosedInterval(tsteps.Start(), tsteps.Stop())).value(); - x_int++; - y_int++; - } - if (_c_index.has_value()){ - read_transform = (std::move(read_transform) | tensorstore::Dims(_c_index.value()).ClosedInterval(channels.Start(), channels.Stop())).value(); - x_int++; - y_int++; - } - if (_z_index.has_value()){ - read_transform = (std::move(read_transform) | tensorstore::Dims(_z_index.value()).ClosedInterval(layers.Start(), layers.Stop())).value(); - x_int++; - y_int++; - } - read_transform = (std::move(read_transform) | tensorstore::Dims(y_int).ClosedInterval(rows.Start(), rows.Stop()) | - tensorstore::Dims(x_int).ClosedInterval(cols.Start(), cols.Stop())).value(); - - - tensorstore::Read(source | read_transform, tensorstore::UnownedToShared(array)).value(); - - return read_buffer; -} - -std::tuple, std::vector> PyramidCompositor::ReadImageData( - const std::string& fname, - const Seq& rows, const Seq& cols, - const Seq& layers = Seq(0,0), - const Seq& channels = Seq(0,0), - const Seq& tsteps = Seq(0,0)) { - - auto read_spec = GetZarrSpecToRead(fname); - - tensorstore::TensorStore source; - - TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( - read_spec, - tensorstore::OpenMode::open, - tensorstore::ReadWriteMode::read).result()); - - auto image_shape = source.domain().shape(); - const auto read_chunk_shape = source.chunk_layout().value().read_chunk_shape(); - - // if (image_shape.size() == 5){ - // _image_height = image_shape[3]; - // _image_width = image_shape[4]; - // _image_depth = image_shape[2]; - // _num_channels = image_shape[1]; - // _num_tsteps = image_shape[0]; - // _t_index.emplace(0); - // _c_index.emplace(1); - // _z_index.emplace(2); - // } else { - // assert(image_shape.size() >= 2); - // std::tie(_t_index, _c_index, _z_index) = ParseMultiscaleMetadata(axes_list, image_shape.size()); - // _image_height = image_shape[image_shape.size()-2]; - // _image_width = image_shape[image_shape.size()-1]; - // if (_t_index.has_value()) { - // _num_tsteps = image_shape[_t_index.value()]; - // } else { - // _num_tsteps = 1; - // } - - // if (_c_index.has_value()) { - // _num_channels = image_shape[_c_index.value()]; - // } else { - // _num_channels = 1; - // } - - // if (_z_index.has_value()) { - // _image_depth = image_shape[_z_index.value()]; - // } else { - // _image_depth = 1; - // } - // } - - _image_dtype = source.dtype().name(); - _data_type_code = GetDataTypeCode(_image_dtype); - - switch (_data_type_code) - { - case (1): - return std::make_tuple( - std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), - std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); - case (2): - return std::make_tuple( - std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), - std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); - case (4): - return std::make_tuple( - std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), - std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); - case (8): - return std::make_tuple( - std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), - std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); - case (16): - return std::make_tuple( - std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), - std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); - case (32): - return std::make_tuple( - std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), - std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); - case (64): - return std::make_tuple( - std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), - std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); - case (128): - return std::make_tuple( - std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), - std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); - case (256): - return std::make_tuple( - std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), - std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); - case (512): - return std::make_tuple( - std::make_shared(std::move(*(GetImageDataTemplated(rows, cols, layers, channels, tsteps)))), - std::make_tuple(image_shape[image_shape.size()-2], image_shape[image_shape.size()-1])); - default: - break; - } -} -} // ns argolid \ No newline at end of file From 3a8f28935c775c0d78fdbe4a8d97bebad7c12b05 Mon Sep 17 00:00:00 2001 From: JesseMckinzie Date: Wed, 22 Jan 2025 10:24:09 -0700 Subject: [PATCH 10/24] Remove duplicate --- CMakeLists.txt | 1 - 1 file changed, 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 3376210..ccc1201 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -33,7 +33,6 @@ set(SOURCE src/cpp/core/ome_tiff_to_chunked_pyramid.cpp src/cpp/core/pyramid_compositor.cpp src/cpp/core/pyramid_view.cpp - src/cpp/core/pyramid_view.cpp src/cpp/utilities/utilities.cpp ) From 50589e34e996bafbffb5d542518d61d0898f7eb0 Mon Sep 17 00:00:00 2001 From: JesseMckinzie Date: Wed, 22 Jan 2025 10:24:51 -0700 Subject: [PATCH 11/24] Remove unused variant --- src/cpp/core/pyramid_compositor.h | 11 ----------- 1 file changed, 11 deletions(-) diff --git a/src/cpp/core/pyramid_compositor.h b/src/cpp/core/pyramid_compositor.h index 6ef411c..4f3cd6f 100644 --- a/src/cpp/core/pyramid_compositor.h +++ b/src/cpp/core/pyramid_compositor.h @@ -23,17 +23,6 @@ struct TupleHash { } }; -// using image_data = std::variant, -// std::vector, -// std::vector, -// std::vector, -// std::vector, -// std::vector, -// std::vector, -// std::vector, -// std::vector, -// std::vector>; - using image_data = std::variant, tensorstore::Array, tensorstore::Array, From 7b2a7e08ed357fec74dd87847c21be78cc59f760 Mon Sep 17 00:00:00 2001 From: JesseMckinzie Date: Wed, 22 Jan 2025 10:30:37 -0700 Subject: [PATCH 12/24] Update __init__.py --- src/python/argolid/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/python/argolid/__init__.py b/src/python/argolid/__init__.py index 69b8d46..bb5933b 100644 --- a/src/python/argolid/__init__.py +++ b/src/python/argolid/__init__.py @@ -1,5 +1,5 @@ from .pyramid_generator import PyramidGenerartor, PyramidView, PlateVisualizationMetadata, Downsample -from .pyramid_compositor import PyramidCompositorPY, PyramidCompositor2 +from .pyramid_compositor import PyramidCompositor from .volume_generator import VolumeGenerator, PyramidGenerator3D from . import _version From 19a1e2a2efd330a1741eb8f9fefab2e8da454476 Mon Sep 17 00:00:00 2001 From: JesseMckinzie Date: Wed, 22 Jan 2025 10:34:26 -0700 Subject: [PATCH 13/24] Remove unused variant --- src/cpp/core/pyramid_compositor.h | 43 ------------------------------- 1 file changed, 43 deletions(-) diff --git a/src/cpp/core/pyramid_compositor.h b/src/cpp/core/pyramid_compositor.h index 4f3cd6f..b55b62a 100644 --- a/src/cpp/core/pyramid_compositor.h +++ b/src/cpp/core/pyramid_compositor.h @@ -23,49 +23,6 @@ struct TupleHash { } }; -using image_data = std::variant, - tensorstore::Array, - tensorstore::Array, - tensorstore::Array, - tensorstore::Array, - tensorstore::Array, - tensorstore::Array, - tensorstore::Array, - tensorstore::Array, - tensorstore::Array, - - tensorstore::Array, - tensorstore::Array, - tensorstore::Array, - tensorstore::Array, - tensorstore::Array, - tensorstore::Array, - tensorstore::Array, - tensorstore::Array, - tensorstore::Array, - tensorstore::Array, - - tensorstore::Array, - tensorstore::Array, - tensorstore::Array, - tensorstore::Array, - tensorstore::Array, - tensorstore::Array, - tensorstore::Array, - tensorstore::Array, - tensorstore::Array, - tensorstore::Array, - - tensorstore::Array, - tensorstore::Array, - tensorstore::Array, - tensorstore::Array, - tensorstore::Array, - tensorstore::Array, - tensorstore::Array, - tensorstore::Array, - tensorstore::Array, - tensorstore::Array>; class Seq { private: From a2904a23ece74436dd0bb39300f8c1d47b84754a Mon Sep 17 00:00:00 2001 From: JesseMckinzie Date: Wed, 22 Jan 2025 11:16:18 -0700 Subject: [PATCH 14/24] Move TupleHash to utilities --- src/cpp/core/pyramid_compositor.h | 12 ++---------- src/cpp/utilities/utilities.h | 9 +++++++++ 2 files changed, 11 insertions(+), 10 deletions(-) diff --git a/src/cpp/core/pyramid_compositor.h b/src/cpp/core/pyramid_compositor.h index b55b62a..31ad2a6 100644 --- a/src/cpp/core/pyramid_compositor.h +++ b/src/cpp/core/pyramid_compositor.h @@ -9,20 +9,12 @@ #include "tensorstore/tensorstore.h" #include "tensorstore/spec.h" +#include "../utilities/utilities.h" + namespace fs = std::filesystem; namespace argolid { static constexpr int CHUNK_SIZE = 1024; -// custom hashing for using std::tuple as key -struct TupleHash { - template - std::size_t operator()(const std::tuple& t) const { - return std::apply([](const auto&... args) { - return (std::hash>{}(args) ^ ...); - }, t); - } -}; - class Seq { private: diff --git a/src/cpp/utilities/utilities.h b/src/cpp/utilities/utilities.h index 3554956..0c74ca8 100644 --- a/src/cpp/utilities/utilities.h +++ b/src/cpp/utilities/utilities.h @@ -74,4 +74,13 @@ void WriteTSZattrFilePlateImage( std::tuple, std::optional, std::optional>ParseMultiscaleMetadata(const std::string& axes_list, int len); +// custom hashing for using std::tuple as key +struct TupleHash { + template + std::size_t operator()(const std::tuple& t) const { + return std::apply([](const auto&... args) { + return (std::hash>{}(args) ^ ...); + }, t); + } +}; } // ns argolid \ No newline at end of file From 4e3960d9aaeadfb92522cd55e3e60f5fc6fa8796 Mon Sep 17 00:00:00 2001 From: JesseMckinzie Date: Wed, 22 Jan 2025 11:16:26 -0700 Subject: [PATCH 15/24] Fix windows build error --- src/cpp/core/pyramid_compositor.cpp | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/src/cpp/core/pyramid_compositor.cpp b/src/cpp/core/pyramid_compositor.cpp index 7a1cde3..cd86d10 100644 --- a/src/cpp/core/pyramid_compositor.cpp +++ b/src/cpp/core/pyramid_compositor.cpp @@ -1,5 +1,4 @@ #include "pyramid_compositor.h" -#include "../utilities/utilities.h" #include #include @@ -31,7 +30,7 @@ PyramidCompositor::PyramidCompositor(const std::string& input_pyramids_loc, cons void PyramidCompositor::create_xml() { CreateXML( - this->_output_pyramid_name, + this->_output_pyramid_name.u8string(), this->_ome_metadata_file, this->_plate_image_shapes[0], this->_image_dtype @@ -40,14 +39,14 @@ void PyramidCompositor::create_xml() { void PyramidCompositor::create_zattr_file() { WriteTSZattrFilePlateImage( - this->_output_pyramid_name, + this->_output_pyramid_name.u8string(), this->_output_pyramid_name.u8string() + "/data.zarr/0", this->_plate_image_shapes ); } void PyramidCompositor::create_zgroup_file() { - WriteVivZgroupFiles(this->_output_pyramid_name); + WriteVivZgroupFiles(this->_output_pyramid_name.u8string()); } void PyramidCompositor::create_auxiliary_files() { @@ -253,7 +252,7 @@ void PyramidCompositor::set_composition(const std::unordered_map source; From 50113a16c05beaa89a89954f0e38abecbf5d0b5a Mon Sep 17 00:00:00 2001 From: JesseMckinzie Date: Thu, 23 Jan 2025 12:40:47 -0700 Subject: [PATCH 16/24] Add multithreading to write_zarr_chunk --- src/cpp/core/pyramid_compositor.cpp | 79 +++++++++++++++++++---------- src/cpp/core/pyramid_compositor.h | 9 ++-- 2 files changed, 58 insertions(+), 30 deletions(-) diff --git a/src/cpp/core/pyramid_compositor.cpp b/src/cpp/core/pyramid_compositor.cpp index cd86d10..f0b2c6b 100644 --- a/src/cpp/core/pyramid_compositor.cpp +++ b/src/cpp/core/pyramid_compositor.cpp @@ -163,6 +163,8 @@ void PyramidCompositor::_write_zarr_chunk(int level, int channel, int y_index, i int row_start_pos = y_start; + //int row, local_y_start, tile_y_start, tile_y_dim, tile_y_end, col_start_pos; + //int col, local_x_start, tile_x_start, tile_x_dim, tile_x_end; while (row_start_pos < y_end) { int row = row_start_pos / unit_image_height; @@ -181,47 +183,70 @@ void PyramidCompositor::_write_zarr_chunk(int level, int channel, int y_index, i int tile_x_dim = std::min((col + 1) * unit_image_width - col_start_pos, x_end - col_start_pos); int tile_x_end = tile_x_start + tile_x_dim; - // Fetch input file path from composition map - std::string input_file_name = _composition_map[{col, row, channel}]; - std::filesystem::path zarrArrayLoc = std::filesystem::path(input_file_name) / "data.zarr/0/" / std::to_string(level); + _th_pool.detach_task([ + row, + col, + channel, + level, + x_start, + x_end, + local_x_start, + tile_x_start, + tile_x_dim, + tile_x_end, + local_y_start, + tile_y_start, + tile_y_dim, + tile_y_end, + y_start, + y_end, + &assembled_image_vec, + assembled_width, + this + ]() { + // Fetch input file path from composition map + std::string input_file_name = _composition_map[{col, row, channel}]; + std::filesystem::path zarrArrayLoc = std::filesystem::path(input_file_name) / "data.zarr/0/" / std::to_string(level); + + auto read_spec = GetZarrSpecToRead(zarrArrayLoc.u8string()); - auto read_spec = GetZarrSpecToRead(zarrArrayLoc.u8string()); - - tensorstore::TensorStore source; + tensorstore::TensorStore source; - TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( - read_spec, - tensorstore::OpenMode::open, - tensorstore::ReadWriteMode::read).result()); + TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( + read_spec, + tensorstore::OpenMode::open, + tensorstore::ReadWriteMode::read).result()); - std::vector read_buffer((tile_x_end-tile_x_start)*(tile_y_end-tile_y_start)); - auto array = tensorstore::Array(read_buffer.data(), {tile_y_end-tile_y_start, tile_x_end-tile_x_start}, tensorstore::c_order); + std::vector read_buffer((tile_x_end-tile_x_start)*(tile_y_end-tile_y_start)); + auto array = tensorstore::Array(read_buffer.data(), {tile_y_end-tile_y_start, tile_x_end-tile_x_start}, tensorstore::c_order); - tensorstore::IndexTransform<> read_transform = tensorstore::IdentityTransform(source.domain()); + tensorstore::IndexTransform<> read_transform = tensorstore::IdentityTransform(source.domain()); - Seq rows = Seq(y_start, y_end); - Seq cols = Seq(x_start, x_end); - Seq tsteps = Seq(0, 0); - Seq channels = Seq(channel-1, channel); - Seq layers = Seq(0, 0); - - read_transform = (std::move(read_transform) | tensorstore::Dims(_y_index).ClosedInterval(rows.Start(), rows.Stop()-1) | - tensorstore::Dims(_x_index).ClosedInterval(cols.Start(), cols.Stop()-1)).value(); + Seq rows = Seq(y_start, y_end); + Seq cols = Seq(x_start, x_end); + Seq tsteps = Seq(0, 0); + Seq channels = Seq(channel-1, channel); + Seq layers = Seq(0, 0); + + read_transform = (std::move(read_transform) | tensorstore::Dims(_y_index).ClosedInterval(rows.Start(), rows.Stop()-1) | + tensorstore::Dims(_x_index).ClosedInterval(cols.Start(), cols.Stop()-1)).value(); - tensorstore::Read(source | read_transform, tensorstore::UnownedToShared(array)).value(); + tensorstore::Read(source | read_transform, tensorstore::UnownedToShared(array)).value(); - for (int i = local_y_start; i < local_y_start + tile_y_end - tile_y_start; ++i) { - for (int j = local_x_start; j < local_x_start + tile_x_end - tile_x_start; ++j) { - assembled_image_vec[i * assembled_width + j] = read_buffer[(i - local_y_start) * (tile_x_end - tile_x_start) + (j - local_x_start)]; + for (int i = local_y_start; i < local_y_start + tile_y_end - tile_y_start; ++i) { + for (int j = local_x_start; j < local_x_start + tile_x_end - tile_x_start; ++j) { + assembled_image_vec[i * assembled_width + j] = read_buffer[(i - local_y_start) * (tile_x_end - tile_x_start) + (j - local_x_start)]; + } } - } + }); col_start_pos += tile_x_end - tile_x_start; } row_start_pos += tile_y_end - tile_y_start; } - //auto& write_source = _zarr_arrays[level]; + // wait for threads to finish assembling vector before writing + _th_pool.wait(); WriteImageData( _zarr_arrays[level], diff --git a/src/cpp/core/pyramid_compositor.h b/src/cpp/core/pyramid_compositor.h index 31ad2a6..916edcd 100644 --- a/src/cpp/core/pyramid_compositor.h +++ b/src/cpp/core/pyramid_compositor.h @@ -5,6 +5,7 @@ #include #include #include +#include "BS_thread_pool.hpp" #include "tensorstore/tensorstore.h" #include "tensorstore/spec.h" @@ -38,14 +39,14 @@ class PyramidCompositor { void create_zgroup_file(); void create_auxiliary_files(); - template - void _write_zarr_chunk(int level, int channel, int y_index, int x_index); - void write_zarr_chunk(int level, int channel, int y_index, int x_index); void set_composition(const std::unordered_map, std::string, TupleHash>& comp_map); private: + + template + void _write_zarr_chunk(int level, int channel, int y_index, int x_index); template void WriteImageData( @@ -80,5 +81,7 @@ class PyramidCompositor { int _x_index = 4; int _y_index = 3; int _c_index = 1; + + BS::thread_pool _th_pool; }; } // ns argolid From 8124eb79ca75f13e51ae4874baf593d2bfe3c0ac Mon Sep 17 00:00:00 2001 From: JesseMckinzie Date: Thu, 23 Jan 2025 13:47:02 -0700 Subject: [PATCH 17/24] Add multithreading for set_composition and store reader in unordered_map --- src/cpp/core/pyramid_compositor.cpp | 103 ++++++++++++++++++---------- src/cpp/core/pyramid_compositor.h | 2 + 2 files changed, 67 insertions(+), 38 deletions(-) diff --git a/src/cpp/core/pyramid_compositor.cpp b/src/cpp/core/pyramid_compositor.cpp index f0b2c6b..ea67f2e 100644 --- a/src/cpp/core/pyramid_compositor.cpp +++ b/src/cpp/core/pyramid_compositor.cpp @@ -212,10 +212,16 @@ void PyramidCompositor::_write_zarr_chunk(int level, int channel, int y_index, i tensorstore::TensorStore source; - TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( - read_spec, - tensorstore::OpenMode::open, - tensorstore::ReadWriteMode::read).result()); + if (_zarr_readers.find(zarrArrayLoc.u8string()) != _zarr_readers.end()) { + source = _zarr_readers[zarrArrayLoc.u8string()]; + } else { + TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( + read_spec, + tensorstore::OpenMode::open, + tensorstore::ReadWriteMode::read).result()); + + _zarr_readers[zarrArrayLoc.u8string()] = source; + } std::vector read_buffer((tile_x_end-tile_x_start)*(tile_y_end-tile_y_start)); auto array = tensorstore::Array(read_buffer.data(), {tile_y_end-tile_y_start, tile_x_end-tile_x_start}, tensorstore::c_order); @@ -274,32 +280,45 @@ void PyramidCompositor::set_composition(const std::unordered_map source; - - TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( - read_spec, - tensorstore::OpenMode::open, - tensorstore::ReadWriteMode::read).result()); - - auto image_shape = source.domain().shape(); - _image_ts_dtype = source.dtype(); - _image_dtype = source.dtype().name(); - _image_dtype_code = GetDataTypeCode(_image_dtype); - - _unit_image_shapes[std::stoi(res_key)] = { - image_shape[image_shape.size()-2], - image_shape[image_shape.size()-1] - }; + _th_pool.detach_task([ + dic, + file_path, + this + ](){ + std::string res_key = dic["path"]; + + std::filesystem::path zarr_array_loc = std::filesystem::path(file_path) / "data.zarr/0/" / res_key; + + auto read_spec = GetZarrSpecToRead(zarr_array_loc.u8string()); + + tensorstore::TensorStore source; + + TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( + read_spec, + tensorstore::OpenMode::open, + tensorstore::ReadWriteMode::read).result()); + + // store reader to use later + _zarr_readers[zarr_array_loc.u8string()] = source; + + auto image_shape = source.domain().shape(); + _image_ts_dtype = source.dtype(); + _image_dtype = source.dtype().name(); + _image_dtype_code = GetDataTypeCode(_image_dtype); + + _unit_image_shapes[std::stoi(res_key)] = { + image_shape[image_shape.size()-2], + image_shape[image_shape.size()-1] + }; + }); } break; } } + _th_pool.wait(); + int num_rows = 0, num_cols = 0; _num_channels = 0; for (const auto& coord : comp_map) { @@ -330,22 +349,30 @@ void PyramidCompositor::set_composition(const std::unordered_map source; - - TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( - GetZarrSpecToWrite( - path, - _plate_image_shapes[level], - {1,1,1, CHUNK_SIZE, CHUNK_SIZE}, - ChooseBaseDType(_image_ts_dtype).value().encoded_dtype - ), - tensorstore::OpenMode::create | - tensorstore::OpenMode::delete_existing, - tensorstore::ReadWriteMode::write).result()); - - _zarr_arrays[level] = source; + _th_pool.detach_task([ + path, + level=level, + this + ](){ + tensorstore::TensorStore source; + + TENSORSTORE_CHECK_OK_AND_ASSIGN(source, tensorstore::Open( + GetZarrSpecToWrite( + path, + _plate_image_shapes[level], + {1,1,1, CHUNK_SIZE, CHUNK_SIZE}, + ChooseBaseDType(_image_ts_dtype).value().encoded_dtype + ), + tensorstore::OpenMode::create | + tensorstore::OpenMode::delete_existing, + tensorstore::ReadWriteMode::write).result()); + + _zarr_arrays[level] = source; + }); } + _th_pool.wait(); + create_auxiliary_files(); } diff --git a/src/cpp/core/pyramid_compositor.h b/src/cpp/core/pyramid_compositor.h index 916edcd..df37f96 100644 --- a/src/cpp/core/pyramid_compositor.h +++ b/src/cpp/core/pyramid_compositor.h @@ -66,6 +66,8 @@ class PyramidCompositor { std::unordered_map> _zarr_arrays; + std::unordered_map> _zarr_readers; + std::unordered_map> _unit_image_shapes; std::unordered_map, std::string, TupleHash> _composition_map; From a3c5805d69ab2ddde6d8dffa61b9dfc1fd2f34b0 Mon Sep 17 00:00:00 2001 From: JesseMckinzie Date: Thu, 23 Jan 2025 14:39:58 -0700 Subject: [PATCH 18/24] Update to upload-artifact@v4 --- .github/workflows/wheel_build.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/wheel_build.yml b/.github/workflows/wheel_build.yml index 26bd3b1..e73dce6 100644 --- a/.github/workflows/wheel_build.yml +++ b/.github/workflows/wheel_build.yml @@ -78,7 +78,7 @@ jobs: CIBW_TEST_COMMAND: python -W default -m unittest discover -s {project}/tests -v - name: Upload Artifact - uses: actions/upload-artifact@v3 + uses: actions/upload-artifact@v4 with: name: argolid-wheels path: dist/*.whl @@ -129,7 +129,7 @@ jobs: CIBW_TEST_COMMAND: python -W default -m unittest discover -s {project}/tests -v - name: Upload Artifact - uses: actions/upload-artifact@v3 + uses: actions/upload-artifact@v4 with: name: argolid-wheels-apple-arm64 path: dist/*.whl From 555e708829b66adf2ba4da9303d0bfe2239f4bd7 Mon Sep 17 00:00:00 2001 From: JesseMckinzie Date: Fri, 24 Jan 2025 08:18:15 -0700 Subject: [PATCH 19/24] Update artifact names --- .github/workflows/wheel_build.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/wheel_build.yml b/.github/workflows/wheel_build.yml index e73dce6..8807da7 100644 --- a/.github/workflows/wheel_build.yml +++ b/.github/workflows/wheel_build.yml @@ -80,7 +80,7 @@ jobs: - name: Upload Artifact uses: actions/upload-artifact@v4 with: - name: argolid-wheels + name: argolid-wheels-${{ matrix.os }}-${{ matrix.cibw_archs }}-${{ matrix.cibw_build }} path: dist/*.whl retention-days: 1 @@ -131,6 +131,6 @@ jobs: - name: Upload Artifact uses: actions/upload-artifact@v4 with: - name: argolid-wheels-apple-arm64 + name: argolid-wheels-${{ matrix.os }}-${{ matrix.cibw_archs }}-${{ matrix.cibw_build }} path: dist/*.whl retention-days: 1 \ No newline at end of file From 0075b109a5853ae736e2f0ed79eddb2a293449fa Mon Sep 17 00:00:00 2001 From: JesseMckinzie Date: Fri, 24 Jan 2025 08:49:09 -0700 Subject: [PATCH 20/24] Update wheel_build.yml --- .github/workflows/wheel_build.yml | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/.github/workflows/wheel_build.yml b/.github/workflows/wheel_build.yml index 8807da7..897043d 100644 --- a/.github/workflows/wheel_build.yml +++ b/.github/workflows/wheel_build.yml @@ -14,7 +14,7 @@ jobs: matrix: os: [ubuntu-20.04, macos-13, windows-latest] cibw_archs: ["auto64"] - cibw_build: ["cp38-*", "cp39-*", "cp310-*", "cp311-*"] + cibw_build: ["cp38", "cp39", "cp310", "cp311"] steps: - uses: actions/checkout@v3 @@ -34,7 +34,7 @@ jobs: - uses: actions/setup-python@v4 name: Install Python with: - python-version: '3.9' + python-version: '3.10' - name: Install cibuildwheel run: | @@ -44,7 +44,7 @@ jobs: run: | python -m cibuildwheel --output-dir dist env: - CIBW_BUILD: ${{ matrix.cibw_build }} + CIBW_BUILD: ${{ matrix.cibw_build }}-* CIBW_BUILD_VERBOSITY: 3 CIBW_SKIP: "*musllinux*" CIBW_MANYLINUX_X86_64_IMAGE: manylinux2014 @@ -93,7 +93,7 @@ jobs: matrix: os: [macos-13-xlarge] cibw_archs: ["arm64"] - cibw_build: ["cp39-*", "cp310-*", "cp311-*"] + cibw_build: ["cp39", "cp310", "cp311"] steps: - uses: actions/checkout@v3 @@ -112,7 +112,7 @@ jobs: run: | python -m cibuildwheel --output-dir dist env: - CIBW_BUILD: ${{ matrix.cibw_build }} + CIBW_BUILD: ${{ matrix.cibw_build }}-* CIBW_BUILD_VERBOSITY: 3 CIBW_ARCHS_MACOS: arm64 CIBW_BEFORE_ALL_MACOS: brew install nasm libjpeg-turbo && From 9e7c42d0a8f091dc70f93e3fe2c724d5a9ac4b09 Mon Sep 17 00:00:00 2001 From: JesseMckinzie Date: Mon, 24 Feb 2025 14:06:10 -0700 Subject: [PATCH 21/24] Replace throw with logging --- src/cpp/core/pyramid_compositor.cpp | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/src/cpp/core/pyramid_compositor.cpp b/src/cpp/core/pyramid_compositor.cpp index ea67f2e..5e01be4 100644 --- a/src/cpp/core/pyramid_compositor.cpp +++ b/src/cpp/core/pyramid_compositor.cpp @@ -64,29 +64,35 @@ void PyramidCompositor::write_zarr_chunk(int level, int channel, int y_index, in } if (_composition_map.empty()) { - throw std::runtime_error("No composition map is set. Unable to generate pyramid"); + std::cerr << "No composition map is set. Unable to generate pyramid" << std::endl; + return; } if (_unit_image_shapes.find(level) == _unit_image_shapes.end()) { - throw std::invalid_argument("Requested level (" + std::to_string(level) + ") does not exist"); + std::cerr << "Requested level (" + std::to_string(level) + ") does not exist" << std::endl; + return; } if (channel >= _num_channels) { - throw std::invalid_argument("Requested channel (" + std::to_string(channel) + ") does not exist"); + std::cerr << "Requested channel (" + std::to_string(channel) + ") does not exist" << std::endl; + return; } auto plate_shape_it = _plate_image_shapes.find(level); if (plate_shape_it == _plate_image_shapes.end()) { - throw std::invalid_argument("No plate image shapes found for level (" + std::to_string(level) + ")"); + std::cerr << "No plate image shapes found for level (" + std::to_string(level) + ")" << std::endl; + return; } const auto& plate_shape = plate_shape_it->second; if (plate_shape.size() < 4 || y_index > (plate_shape[3] / CHUNK_SIZE)) { - throw std::invalid_argument("Requested y index (" + std::to_string(y_index) + ") does not exist"); + std::cerr << "Requested y index (" + std::to_string(y_index) + ") does not exist" << std::endl; + return; } if (plate_shape.size() < 5 || x_index > (plate_shape[4] / CHUNK_SIZE)) { - throw std::invalid_argument("Requested x index (" + std::to_string(x_index) + ") does not exist"); + tstd::cerr << "Requested x index (" + std::to_string(x_index) + ") does not exist" << std::endl; + return; } // get datatype by opening From 72978ce95f34dd01291392f9b57794cae5084486 Mon Sep 17 00:00:00 2001 From: JesseMckinzie Date: Mon, 24 Feb 2025 14:10:23 -0700 Subject: [PATCH 22/24] Typo fix --- src/cpp/core/pyramid_compositor.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/cpp/core/pyramid_compositor.cpp b/src/cpp/core/pyramid_compositor.cpp index 5e01be4..e66e4fa 100644 --- a/src/cpp/core/pyramid_compositor.cpp +++ b/src/cpp/core/pyramid_compositor.cpp @@ -91,7 +91,7 @@ void PyramidCompositor::write_zarr_chunk(int level, int channel, int y_index, in } if (plate_shape.size() < 5 || x_index > (plate_shape[4] / CHUNK_SIZE)) { - tstd::cerr << "Requested x index (" + std::to_string(x_index) + ") does not exist" << std::endl; + std::cerr << "Requested x index (" + std::to_string(x_index) + ") does not exist" << std::endl; return; } From e50721278d78a4e451ac4928c5d0ed2cb6b4e210 Mon Sep 17 00:00:00 2001 From: sameeul Date: Tue, 25 Feb 2025 09:51:58 -0500 Subject: [PATCH 23/24] Update wheel building workflow --- .github/workflows/wheel_build.yml | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/.github/workflows/wheel_build.yml b/.github/workflows/wheel_build.yml index d5f8325..93993f7 100644 --- a/.github/workflows/wheel_build.yml +++ b/.github/workflows/wheel_build.yml @@ -14,8 +14,7 @@ jobs: matrix: os: [ubuntu-20.04, macos-13, windows-latest] cibw_archs: ["auto64"] - cibw_build: ["cp38", "cp39", "cp310", "cp311"] - + cibw_build: ["cp39", "cp310", "cp311", "cp312", "cp313"] steps: - uses: actions/checkout@v3 @@ -94,8 +93,7 @@ jobs: matrix: os: [macos-13-xlarge] cibw_archs: ["arm64"] - cibw_build: ["cp39", "cp310", "cp311"] - + cibw_build: ["cp39", "cp310", "cp311", "cp312", "cp313"] steps: - uses: actions/checkout@v3 @@ -134,6 +132,6 @@ jobs: - name: Upload Artifact uses: actions/upload-artifact@v4 with: - name: argolid-wheels-${{ matrix.os }}-${{ matrix.cibw_archs }}-${{ matrix.cibw_build }} + name: argolid-wheels-apple-arm64-${{ matrix.os }}-${{ matrix.cibw_archs }}-${{ matrix.cibw_build }} path: dist/*.whl retention-days: 1 \ No newline at end of file From cdf49e50597419ab5fa8f8a196dc4f4dcda55278 Mon Sep 17 00:00:00 2001 From: JesseMckinzie Date: Tue, 25 Feb 2025 08:34:57 -0700 Subject: [PATCH 24/24] Use std::memcpy for copying buffer to assembled image --- src/cpp/core/pyramid_compositor.cpp | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/cpp/core/pyramid_compositor.cpp b/src/cpp/core/pyramid_compositor.cpp index e66e4fa..95b2f5b 100644 --- a/src/cpp/core/pyramid_compositor.cpp +++ b/src/cpp/core/pyramid_compositor.cpp @@ -245,10 +245,12 @@ void PyramidCompositor::_write_zarr_chunk(int level, int channel, int y_index, i tensorstore::Read(source | read_transform, tensorstore::UnownedToShared(array)).value(); - for (int i = local_y_start; i < local_y_start + tile_y_end - tile_y_start; ++i) { - for (int j = local_x_start; j < local_x_start + tile_x_end - tile_x_start; ++j) { - assembled_image_vec[i * assembled_width + j] = read_buffer[(i - local_y_start) * (tile_x_end - tile_x_start) + (j - local_x_start)]; - } + for (int i = 0; i < tile_y_end - tile_y_start; ++i) { + std::memcpy( + &assembled_image_vec[(local_y_start + i) * assembled_width + local_x_start], // Destination + &read_buffer[i * (tile_x_end - tile_x_start)], // Source + (tile_x_end - tile_x_start) * sizeof(assembled_image_vec[0]) // Number of bytes + ); } });