Implement findBlockBoundaryEdges WIP

include/nifty/distributed/edge_morphology.hxx

@@ -1,5 +1,6 @@
 #pragma once
 
+#include "nifty/tools/blocking.hxx"
 #include "nifty/distributed/mergeable_features.hxx"
 
 namespace fs = boost::filesystem;
@@ -10,19 +11,24 @@ namespace distributed {
 
 
     template<class COORD>
-    inline void loadBB(const COORD & begIn, const COORD & endIn,
-                       std::vector<std::size_t> & begOut, std::vector<std::size_t> & endOut) {
-
-        std::copy(begIn.begin(), begIn.end(), begOut.begin());
-        std::copy(endIn.begin(), endIn.end(), endOut.begin());
-
+    inline void increaseRoi(COORD & roiBegin) {
         // we increase the roi by decreasing roiBegin by 1.
         // to match what was done in the graph extraction when increaseRoi is true
         for(int axis = 0; axis < 3; ++axis) {
-            if(begOut[axis] > 0) {
-                --begOut[axis];
+            if(roiBegin[axis] > 0) {
+                --roiBegin[axis];
             }
         }
+
+    }
+
+
+    template<class COORD>
+    inline void loadBB(const COORD & begIn, const COORD & endIn,
+                       std::vector<std::size_t> & begOut, std::vector<std::size_t> & endOut) {
+        std::copy(begIn.begin(), begIn.end(), begOut.begin());
+        std::copy(endIn.begin(), endIn.end(), endOut.begin());
+        increaseRoi(begOut);
     }
 
 
@@ -134,5 +140,123 @@ namespace distributed {
     }
 
 
+    template<class GRAPH, class LABELS, class BLOCKING, class OUT>
+    inline void findBlockBoundaryEdgesBlock(const GRAPH & graph,
+                                            const LABELS & labels,
+                                            const std::vector<EdgeIndexType> & edgeIndices,
+                                            const BLOCKING & blocking,
+                                            const std::size_t blockId,
+                                            const bool ignoreLabel,
+                                            OUT & out) {
+
+
+        typedef nifty::array::StaticArray<int64_t, 3> CoordType;
+        CoordType shape;
+        std::copy(labels.shape().begin(), labels.shape().end(), shape.begin());
+
+        // iterate over the block faces
+        for(unsigned axis = 0; axis < 3; ++axis) {
+            // check if this is a border block == has no lower neighbor in this axis
+            if(blocking.getNeighborId(blockId, axis, true) == -1) {
+                continue;
+            }
+
+            // get the shape of this face
+            CoordType faceShape;
+            for(unsigned d = 0; d < 3; ++d) {
+                faceShape[d] = (d == axis) ? 1 : shape[d];
+            }
+
+            nifty::tools::forEachCoordinate(faceShape, [&](const CoordType & coord) {
+                const NodeType lU = xtensor::read(labels, coord.asStdArray());
+                // check for ignore label
+                if(lU == 0 && ignoreLabel) {
+                    return;
+                }
+
+                CoordType coord2 = coord;
+                coord2[axis] = 1;
+
+                const NodeType lV = xtensor::read(labels, coord2.asStdArray());
+                if(lV == 0 && ignoreLabel) {
+                    return;
+                }
+                if(lU == lV){
+                    return;
+                }
+
+                const auto edgeId = graph.findEdge(lU, lV);
+                const auto globalEdgeId = edgeIndices[edgeId];
+                out(globalEdgeId) = true;
+            });
+        }
+    }
+
+
+    template<class OUT>
+    inline void findBlockBoundaryEdges(const std::string & blockPrefix,
+                                       const std::string & labelPath,
+                                       const std::string & labelKey,
+                                       const std::vector<std::size_t> & blockShape,
+                                       const std::vector<std::size_t> & blockIds,
+                                       OUT & out) {
+        typedef xt::xtensor<NodeType, 3> LabelArray;
+        typedef nifty::array::StaticArray<int64_t, 3> VectorType;
+
+        fs::path labelsSetPath(labelPath);
+        labelsSetPath /= labelKey;
+
+        auto labelDs = z5::openDataset(labelsSetPath.string());
+        const std::vector<std::string> keys = {"ignoreLabel"};
+
+        VectorType shape, blockShapeVec;
+        std::copy(labelDs->shape().begin(), labelDs->shape().end(), shape.begin());
+        std::copy(blockShape.begin(), blockShape.end(), blockShapeVec.begin());
+
+        tools::Blocking<3> blocking(VectorType({0, 0, 0}),
+                                    shape, blockShapeVec);
+
+        for(const std::size_t blockId : blockIds) {
+
+            // load the sub-graph corresponding to this block
+            const std::string blockPath = blockPrefix + std::to_string(blockId);
+            Graph graph(blockPath);
+            // continue if we don't have edges in this graph
+            if(graph.numberOfEdges() == 0) {
+                continue;
+            }
+
+            const auto block = blocking.getBlock(blockId);
+            auto roiBegin = block.begin();
+            const auto & roiEnd = block.end();
+            increaseRoi(roiBegin);
+
+            // load the bounding box information
+            z5::handle::Group group(blockPath);
+            nlohmann::json j;
+            z5::readAttributes(group, keys, j);
+
+            // get the shape and create array
+            Shape3Type thisShape;
+            for(unsigned axis = 0; axis < 3; ++axis) {
+                thisShape[axis] = roiEnd[axis] - roiBegin[axis];
+            }
+            LabelArray labels(thisShape);
+
+            // load the labels from the bounding box
+            z5::multiarray::readSubarray<NodeType>(labelDs, labels, roiBegin.begin());
+            const bool ignoreLabel = j[keys[0]];
+
+            // load the global edge indices for this block
+            std::vector<EdgeIndexType> edgeIndices;
+            loadEdgeIndices(blockPath, edgeIndices, 0);
+
+            // find which edges are 1d
+            findBlockBoundaryEdgesBlock(graph, labels, edgeIndices, blocking, blockId,
+                                        ignoreLabel, out);
+        }
+    }
+
+
 }
 }

src/python/lib/distributed/edge_morphology.cxx

@@ -20,13 +20,29 @@ namespace distributed {
             {
                 py::gil_scoped_release allowThreads;
                 find1DEdges(blockPrefix, labelPath, labelKey, blockIds, out);
-                std::cout << "blub" << std::endl;
             }
-            std::cout << "blob" << std::endl;
             return out;
         }, py::arg("blockPrefix"),
            py::arg("labelPath"), py::arg("labelKey"),
            py::arg("numberOfEdges"), py::arg("blockIds"));
+
+
+        module.def("findBlockBoundaryEdges", [](const std::string & blockPrefix,
+                                                const std::string & labelPath,
+                                                const std::string & labelKey,
+                                                const std::size_t numberOfEdges,
+                                                const std::vector<std::size_t> & blockShape,
+                                                const std::vector<std::size_t> & blockIds) {
+            xt::pytensor<bool, 1> out = xt::zeros<bool>({numberOfEdges});
+            {
+                py::gil_scoped_release allowThreads;
+                findBlockBoundaryEdges(blockPrefix, labelPath, labelKey, blockShape, blockIds, out);
+            }
+            return out;
+        }, py::arg("blockPrefix"),
+           py::arg("labelPath"), py::arg("labelKey"),
+           py::arg("numberOfEdges"), py::arg("blockShape"),
+           py::arg("blockIds"));
     }
 
 }