Merge pull request #126 from constantinpape/master

Fix pybindings issues + stuff

include/nifty/distributed/distributed_utils.hxx

@@ -77,6 +77,7 @@ namespace distributed {
         }
 
         // write serialization
+        // FIXME not parallelization save ???
         auto ds = z5::openDataset(dsPath);
         ds->writeChunk(chunkId, &serialization[0], true, serSize);
     }

include/nifty/distributed/edge_morphology.hxx

@@ -0,0 +1,262 @@
+#pragma once
+
+#include "nifty/tools/blocking.hxx"
+#include "nifty/distributed/mergeable_features.hxx"
+
+namespace fs = boost::filesystem;
+
+
+namespace nifty {
+namespace distributed {
+
+
+    template<class COORD>
+    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(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);
+    }
+
+
+    template<class GRAPH, class LABELS, class OUT>
+    inline void find1DEdgesBlock(const GRAPH & graph, const LABELS & labels,
+                                 const std::vector<EdgeIndexType> & edgeIndices,
+                                 const bool ignoreLabel, OUT & out, std::vector<uint8_t> & edgeAxes) {
+        typedef nifty::array::StaticArray<int64_t, 3> CoordType;
+        CoordType shape;
+        std::copy(labels.shape().begin(), labels.shape().end(), shape.begin());
+
+        nifty::tools::forEachCoordinate(shape, [&](const CoordType & coord) {
+            const NodeType lU = xtensor::read(labels, coord.asStdArray());
+
+            // check for ignore label
+            if(lU == 0 && ignoreLabel) {
+                return;
+            }
+
+            CoordType coord2;
+            for(std::size_t axis = 0; axis < 3; ++axis){
+                makeCoord2(coord, coord2, axis);
+
+                // check if this is an actual edge
+                if(coord2[axis] >= shape[axis]){
+                    continue;
+                }
+                const NodeType lV = xtensor::read(labels, coord2.asStdArray());
+                if(lV == 0 && ignoreLabel) {
+                    continue;
+                }
+                if(lU == lV){
+                    continue;
+                }
+
+                // yes it is!
+                const auto edgeId = graph.findEdge(lU, lV);
+                const auto globalEdgeId = edgeIndices[edgeId];
+
+                // did we already visit the edge - and if so, what's the prev state?
+                auto & prevState = out(globalEdgeId);
+                // prev state is 0 -> the edge was not visited yet -> we just set our current axis
+                if(prevState == 0) {
+                    edgeAxes[globalEdgeId] = axis;
+                    prevState = 1;  // set prev state to visited and 1d axis
+                } else if(prevState == 1) {  // prev state is 1 -> edge was visited and is still 1d axes
+                    // if the previous axis does not correspond to the current axis, set
+                    // edge status to 2 = mixed axes
+                    if(edgeAxes[globalEdgeId] != axis) {
+                        prevState = 2;
+                    }
+                }
+                // otherwise prevState is 2 -> edge is already mixed axes
+            }
+        });
+    }
+
+
+    template<class OUT>
+    inline void find1DEdges(const std::string & blockPrefix,
+                            const std::string & labelPath,
+                            const std::string & labelKey,
+                            const std::vector<std::size_t> & blockIds,
+                            OUT & out) {
+        typedef xt::xtensor<NodeType, 3> LabelArray;
+        std::vector<uint8_t> edgeAxes(out.size());
+
+        fs::path labelsSetPath(labelPath);
+        labelsSetPath /= labelKey;
+
+        auto labelDs = z5::openDataset(labelsSetPath.string());
+        const std::vector<std::string> keys = {"roiBegin", "roiEnd", "ignoreLabel"};
+
+        std::vector<std::size_t> roiBegin(3), roiEnd(3);
+        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;
+            }
+
+            // load the bounding box information
+            z5::handle::Group group(blockPath);
+            nlohmann::json j;
+            z5::readAttributes(group, keys, j);
+            loadBB(j[keys[0]], j[keys[1]], roiBegin, roiEnd);
+
+            // get the shape and create array
+            Shape3Type shape;
+            for(unsigned axis = 0; axis < 3; ++axis) {
+                shape[axis] = roiEnd[axis] - roiBegin[axis];
+            }
+            LabelArray labels(shape);
+
+            // load the labels from the bounding box
+            z5::multiarray::readSubarray<NodeType>(labelDs, labels, roiBegin.begin());
+            const bool ignoreLabel = j[keys[2]];
+
+            // load the global edge indices for this block
+            std::vector<EdgeIndexType> edgeIndices;
+            loadEdgeIndices(blockPath, edgeIndices, 0);
+
+            // find which edges are 1d
+            find1DEdgesBlock(graph, labels, edgeIndices, ignoreLabel, out, edgeAxes);
+        }
+    }
+
+
+    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);
+        }
+    }
+
+
+}
+}

include/nifty/distributed/graph_extraction.hxx

@@ -402,7 +402,7 @@ namespace distributed {
                     lV = xtensor::read(labels, coord2.asStdArray());
                     // skip zero if we have an ignoreLabel
                     if(ignoreLabel && (lV == 0)) {
-                        return;
+                        continue;
                     }
                     if(lU != lV){
                         edges.insert(std::make_pair(std::min(lU, lV),

include/nifty/distributed/lifted_utils.hxx

@@ -68,12 +68,13 @@ namespace distributed {
     }
 
 
-    template<class NODE_LABELS>
+    template<class NODE_LABELS, class F>
     inline void findLiftedEdgesBfs(const Graph & graph,
                                    const uint64_t srcNode,
                                    const NODE_LABELS & nodeLabels,
                                    const unsigned graphDepth,
-                                   std::vector<EdgeType> & out) {
+                                   std::vector<EdgeType> & out,
+                                   F && addEdge) {
         // type to put on bfs queue, stores the node id and the graph
         // distance of this node from the start node
         typedef std::pair<uint64_t, unsigned> QueueElem;
@@ -84,6 +85,7 @@ namespace distributed {
         std::queue<QueueElem> queue;
 
         queue.emplace(std::make_pair(srcNode, 0));
+        const uint64_t srcLabel = nodeLabels[srcNode];
 
         while(queue.size()) {
             const QueueElem elem = queue.front();
@@ -116,23 +118,61 @@ namespace distributed {
             }
 
             // check if we make a lifted edge between the start node and this node:
-            // is this a lifted edge? i.e. graph depth > 1
+            // 1.) is this a lifted edge? i.e. graph depth > 1
             const bool isLiftedEdge = depth > 1;
-            // does the node have a label?
+            if(depth <= 1) {
+                continue;
+            }
+
+            // 2.) does the node have a label?
             const uint64_t label = nodeLabels[node];
-            // is the node's id bigger than srcNode? (otherwise edges would be redundant)
-            if(isLiftedEdge && label > 0 && srcNode < node) {
-                out.emplace_back(std::make_pair(srcNode, node));
+            if(label == 0) {
+                continue;
+            }
+
+            // 3.) is srcNode < node? (make sure not to duplicate lifted edges)
+            if(srcNode > node) {
+                continue;
+            }
+
+            // 4) additional lifted check dependend on the node labels
+            if(!addEdge(srcLabel, label)) {
+                continue;
             }
+
+            // all checks passed ? -> add the lifted edges
+            out.emplace_back(std::make_pair(srcNode, node));
         }
     }
 
 
+    inline bool addAll(const uint64_t labelA, const uint64_t labelB) {
+        return true;
+    }
+
+
+    inline bool addSame(const uint64_t labelA, const uint64_t labelB) {
+        return labelA == labelB;
+    }
+
+
+    inline bool addDifferent(const uint64_t labelA, const uint64_t labelB) {
+        return labelA != labelB;
+    }
+
+
     inline void computeLiftedNeighborhoodFromNodeLabels(const std::string & graphPath,
                                                         const std::string & nodeLabelPath,
                                                         const std::string & outputPath,
                                                         const unsigned graphDepth,
-                                                        const int numberOfThreads) {
+                                                        const int numberOfThreads,
+                                                        const std::string & mode="all") {
+        // modes can be
+        // "all": add all edges between nodes with a label
+        // "same": add edges only between nodes with the same label
+        // "different": add edges only between nodes with different labels
+        auto edgeChecker = (mode=="all") ? addAll : ((mode=="same") ? addSame : addDifferent);
+
         // load the graph
         const auto graph = Graph(graphPath, numberOfThreads);
 
@@ -160,7 +200,8 @@ namespace distributed {
             auto & threadData = perThreadData[tid];
             // do bfs for this node and find all relevant lifted edges
             findLiftedEdgesBfs(graph, nodeId,
-                               nodeLabels, graphDepth, threadData);
+                               nodeLabels, graphDepth, threadData,
+                               edgeChecker);
         });
 
         // merge the thread