Add masking utility to AnnotationRegion

dlup/annotations.py

@@ -14,7 +14,7 @@
 - HaloXML
 """
 from __future__ import annotations
-
+import time
 import copy
 import errno
 import functools
@@ -1110,6 +1110,7 @@ def _affine_coords(coords: npt.NDArray[np.float_]) -> npt.NDArray[np.float_]:
  for annotation in cropped_annotations:
  annotation = transform(annotation, _affine_coords)
  output.append(annotation)
+
  return output
 
  def __contains__(self, item: Union[str, AnnotationClass]) -> bool:

dlup/annotations_experimental.py

@@ -213,9 +213,7 @@ def as_geojson(self) -> GeoJsonDict:
  return data
 
  def read_region(self, coordinates: tuple[int, int], scaling: float, size: tuple[int, int]):
- start_time = time.time()
  region = self._layers.read_region(coordinates, scaling, size)
- print(f"Time to read region (dlup v0.8.0.beta): {(time.time() - start_time):.5f}s")
  return region
 
  def scale(self, scaling: float) -> None:

gen_polygons.py

@@ -22,7 +22,8 @@
 # Bounding box:
 bbox = annotations.bounding_box
 print(f"Bounding box: {bbox}")
-region_start = (500, 0)
+region_start = (0, 0)
+
 
 # Let's get the region
 start_time = time.time()
@@ -93,12 +94,44 @@
 np.asarray((56630.2124, 69640.6535)) * 0.02
 region_size = (1393, 1133)
 
+
+
 _, mask, _ = convert_annotations(region, region_size=region_size, index_map=index_map)
 print(mask.shape)
 
 
 PIL.Image.fromarray(LUT[mask]).resize((1133 // 2, 1393 // 2)).save("dlup_original.png")
 
+mask_ = region2.to_mask(region_size, index_map, 0)
+PIL.Image.fromarray(LUT[mask_]).resize((1133 // 2, 1393 // 2)).save("dlup_new_opencv.png")
 
 mask3 = convert_annotations_new(region2.polygons, region_size=region_size, index_map=index_map)
 PIL.Image.fromarray(LUT[mask3]).resize((1133 // 2, 1393 // 2)).save("dlup_new.png")
+
+print()
+# Let's time everything separately.
+print("Benchmark\n=========")
+
+annotations = WsiAnnotations.from_geojson(fn, sorting="NONE")
+bbox = annotations.bounding_box
+
+start_time = time.time()
+region = annotations.read_region(region_start, 0.02, bbox[1])
+print(f"Time to read region (dlup v0.7.0): {(time.time() - start_time) * 1000:.2f}ms")
+start_time2 = time.time()
+_, mask, _ = convert_annotations(region, region_size=region_size, index_map=index_map)
+print(f"Time to convert annotations to mask (dlup v0.7.0): {(time.time() - start_time2) * 1000:.2f}ms")
+total_time = (time.time() - start_time)
+print(f"Total time to read region and convert to mask (dlup v0.7.0): {total_time * 1000:.2f}ms")
+print()
+annotations2 = WsiAnnotations2.from_geojson(fn)
+start_time = time.time()
+region2 = annotations2.read_region(region_start, 0.02, bbox[1])
+print(f"Time to read region (dlup v0.8.0.beta): {(time.time() - start_time) * 1000:.2f}ms")
+start_time2 = time.time()
+mask_ = region2.to_mask(region_size, index_map, 0)
+print(f"Time to convert annotations to mask (dlup v0.8.0.beta): {(time.time() - start_time2) * 1000:.2f}ms")
+total_time2 = (time.time() - start_time)
+print(f"Total time to read region and convert to mask (dlup v0.8.0.beta): {total_time2 * 1000:.2f}ms")
+
+print(f"\nSpeedup: {total_time/total_time2:.3f} times")

meson.build

@@ -51,6 +51,9 @@ incdir_pybind11 = include_directories(pybind11_include)
 boost_modules = ['system', 'serialization']
 boost_dep = dependency('boost', modules : boost_modules, required : true)
 
+# OpenCV
+opencv_dep = dependency('opencv4', required : true)
+
 ### End Includes ### 
 
 
@@ -84,4 +87,4 @@ _geometry = py.extension_module('_geometry',
  include_directories : [incdir_pybind11],
  install : true,
  cpp_args : base_cpp_args,
- dependencies : base_deps + boost_dep)
+ dependencies : base_deps + boost_dep + opencv_dep)

src/geometry.cpp

@@ -8,11 +8,16 @@
 #include "exceptions.h"
 #include "geometry.h"
 #include <memory>
+#include <opencv2/imgproc.hpp>
+#include <opencv2/opencv.hpp>
+#include <pybind11/numpy.h>
+#include <pybind11/pybind11.h>
 #include <stdexcept>
 #include <string>
+#include <unordered_map>
 #include <vector>
 
-#define DLUPDEBUG
+// #define DLUPDEBUG
 
 namespace bg = boost::geometry;
 namespace bgi = boost::geometry::index;
@@ -25,6 +30,8 @@ using BoostRing = bg::model::ring<BoostPoint>;
 using BoostLineString = bg::model::linestring<BoostPoint>;
 using BoostMultiPolygon = bg::model::multi_polygon<BoostPolygon>;
 
+namespace py = pybind11;
+
 class FactoryGuard {
 public:
  FactoryGuard(py::function &factory_ref, py::function new_factory)
@@ -282,6 +289,134 @@ class Point : public BaseGeometry {
  }
 };
 
+cv::Mat generateMaskFromAnnotations(const std::vector<std::shared_ptr<Polygon>> &annotations, cv::Size region_size,
+ const std::unordered_map<std::string, int> &index_map, int default_value) {
+ // Create the mask and initialize with the default value
+ cv::Mat mask(region_size, CV_32S, cv::Scalar(default_value));
+
+ std::vector<cv::Point> exterior_cv_points;
+ std::vector<std::vector<cv::Point>> interiors_cv_points;
+
+ for (const auto &annotation : annotations) {
+ int index_value = index_map.at(annotation->getField("label")->cast<std::string>());
+
+ // Convert exterior points
+ exterior_cv_points.clear();
+ const auto &exterior = annotation->getExterior();
+ exterior_cv_points.reserve(exterior.size());
+ for (const auto &[x, y] : exterior) {
+ exterior_cv_points.emplace_back(static_cast<int>(std::round(x)), static_cast<int>(std::round(y)));
+ }
+
+ // Convert interior points
+ interiors_cv_points.clear();
+ const auto &interiors = annotation->getInteriors();
+ interiors_cv_points.reserve(interiors.size());
+ for (const auto &interior : interiors) {
+ std::vector<cv::Point> interior_cv;
+ interior_cv.reserve(interior.size());
+ for (const auto &[x, y] : interior) {
+ interior_cv.emplace_back(static_cast<int>(std::round(x)), static_cast<int>(std::round(y)));
+ }
+ interiors_cv_points.push_back(std::move(interior_cv));
+ }
+
+ // Only clone mask if necessary
+ cv::Mat original_values;
+ if (!interiors_cv_points.empty()) {
+ original_values = mask.clone();
+ }
+
+ // Create a mask for holes if necessary
+ cv::Mat holes_mask;
+ if (!interiors_cv_points.empty()) {
+ holes_mask = cv::Mat::zeros(region_size, CV_8U);
+ cv::fillPoly(holes_mask, interiors_cv_points, cv::Scalar(1));
+ }
+
+ // Fill the exterior polygon in the mask
+ cv::fillPoly(mask, std::vector<std::vector<cv::Point>>{exterior_cv_points}, cv::Scalar(index_value));
+
+ // If interiors exist, reset the holes in the mask using the backup
+ if (!interiors_cv_points.empty()) {
+ original_values.copyTo(mask, holes_mask);
+ }
+ }
+
+ return mask;
+}
+
+// cv::Mat generateMaskFromAnnotations(const std::vector<std::shared_ptr<Polygon>> &annotations, cv::Size region_size,
+// const std::unordered_map<std::string, int> &index_map, int default_value) {
+// // Create the mask and initialize with the default value
+// cv::Mat mask(region_size, CV_32S, cv::Scalar(default_value));
+
+// for (const auto &annotation : annotations) {
+// // Extract the label and map it to an index value
+// int index_value = index_map.at(annotation->getField("label")->cast<std::string>());
+
+// // Convert the exterior and interiors to OpenCV points
+// std::vector<cv::Point> exterior_cv_points;
+// for (const auto &[x, y] : annotation->getExterior()) {
+// exterior_cv_points.emplace_back(static_cast<int>(std::round(x)), static_cast<int>(std::round(y)));
+// }
+
+// std::vector<std::vector<cv::Point>> interiors_cv_points;
+// for (const auto &interior : annotation->getInteriors()) {
+// std::vector<cv::Point> interior_cv;
+// for (const auto &[x, y] : interior) {
+// interior_cv.emplace_back(static_cast<int>(std::round(x)), static_cast<int>(std::round(y)));
+// }
+// interiors_cv_points.push_back(std::move(interior_cv));
+// }
+
+// // Backup original mask values where holes will be drawn
+// cv::Mat original_values = mask.clone();
+// cv::Mat holes_mask = cv::Mat::zeros(region_size, CV_8U);
+// if (!interiors_cv_points.empty()) {
+// cv::fillPoly(holes_mask, interiors_cv_points, cv::Scalar(1));
+// }
+
+// #ifdef DLUPDEBUG
+// // Debug: Check matrix types and sizes before the setTo operation
+// std::cout << "mask type: " << mask.type() << ", size: " << mask.size << std::endl;
+// std::cout << "original_values type: " << original_values.type() << ", size: " << original_values.size <<
+// std::endl; std::cout << "holes_mask type: " << holes_mask.type() << ", size: " << holes_mask.size <<
+// std::endl;
+// #endif
+
+// // Fill the exterior polygon in the mask
+// cv::fillPoly(mask, std::vector<std::vector<cv::Point>>{exterior_cv_points}, cv::Scalar(index_value));
+
+// // If interiors exist, reset the holes in the mask using the backup
+// if (!interiors_cv_points.empty()) {
+// original_values.copyTo(mask, holes_mask);
+
+// }
+// }
+
+// return mask;
+// }
+
+py::array_t<int> maskToPyArray(const cv::Mat &mask) {
+ // Ensure the mask is of type CV_32S (int type)
+ if (mask.type() != CV_32S) {
+ throw std::runtime_error("Mask must be of type CV_32S (int).");
+ }
+
+ // Create a buffer info that describes the numpy array
+ py::buffer_info buf_info(mask.data, // Pointer to buffer
+ sizeof(int), // Size of one scalar element
+ py::format_descriptor<int>::format(), // Python struct-style format descriptor
+ 2, // Number of dimensions
+ {mask.rows, mask.cols}, // Buffer dimensions
+ {sizeof(int) * mask.cols, sizeof(int)} // Strides (in bytes) for each dimension
+ );
+
+ // Create the numpy array from the buffer info
+ return py::array_t<int>(buf_info);
+}
+
 class AnnotationRegion {
 public:
  AnnotationRegion(std::vector<std::shared_ptr<Polygon>> polygons, std::vector<std::shared_ptr<Point>> points)
@@ -328,6 +463,22 @@ class AnnotationRegion {
  return py_points;
  }
 
+ py::array_t<int> toMask(std::tuple<int, int> mask_size, const std::unordered_map<std::string, int> &index_map,
+ int default_value = 0) const {
+#ifdef DLUPDEBUG
+ std::chrono::steady_clock::time_point begin = std::chrono::steady_clock::now();
+#endif
+ cv::Size region_size(std::get<1>(mask_size), std::get<0>(mask_size));
+ cv::Mat mask = generateMaskFromAnnotations(polygons_, region_size, index_map, default_value);
+#ifdef DLUPDEBUG
+ std::cout
+ << "AnnotationRegion::toMask: mask generated in "
+ << std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::steady_clock::now() - begin).count()
+ << " ms" << std::endl;
+#endif
+ return maskToPyArray(mask);
+ }
+
 private:
  std::vector<std::shared_ptr<Polygon>> polygons_;
  std::vector<std::shared_ptr<Point>> points_;
@@ -641,7 +792,9 @@ PYBIND11_MODULE(_geometry, m) {
 
  py::class_<AnnotationRegion, std::shared_ptr<AnnotationRegion>>(m, "RegionResult")
  .def_property_readonly("polygons", &AnnotationRegion::getPolygons)
- .def_property_readonly("points", &AnnotationRegion::getPoints);
+ .def_property_readonly("points", &AnnotationRegion::getPoints)
+ .def("to_mask", &AnnotationRegion::toMask, py::arg("mask_size"), py::arg("index_map"),
+ py::arg("default_value") = 0);
 
  py::register_exception<GeometryError>(m, "GeometryError");
  py::register_exception<GeometryIntersectionError>(m, "GeometryIntersectionError");