diff --git a/setup.py b/setup.py
index a90e566..61a24c1 100644
--- a/setup.py
+++ b/setup.py
@@ -1,4 +1,6 @@
 import sys
+import tempfile
+from pathlib import Path
 
 import numpy
 from Cython.Build import cythonize
@@ -24,4 +26,26 @@
     extra_compile_args=["-std=c++11", "-w"],
 )
 
-setup(ext_modules=cythonize([evaluate]))
+# Pre-compile default agglomeration module (OneMinus<MeanAffinity>, PriorityQueue)
+# so it doesn't need JIT compilation at runtime.
+_generated_dir = Path(tempfile.mkdtemp())
+(_generated_dir / "ScoringFunction.h").write_text(
+    "typedef OneMinus<MeanAffinity<RegionGraphType, ScoreValue>> ScoringFunctionType;"
+)
+(_generated_dir / "Queue.h").write_text(
+    "template<typename T, typename S> using QueueType = PriorityQueue<T, S>;"
+)
+
+agglomerate = Extension(
+    name="waterz._agglomerate_default",
+    sources=[
+        "src/waterz/agglomerate.pyx",
+        "src/waterz/frontend_agglomerate.cpp",
+    ],
+    include_dirs=include_dirs + [str(_generated_dir)],
+    language="c++",
+    extra_link_args=["-std=c++11"],
+    extra_compile_args=["-std=c++11", "-w", "-O3"],
+)
+
+setup(ext_modules=cythonize([evaluate, agglomerate]))
diff --git a/src/waterz/_agglomerate.py b/src/waterz/_agglomerate.py
index 63d0f9d..9b0f7bd 100644
--- a/src/waterz/_agglomerate.py
+++ b/src/waterz/_agglomerate.py
@@ -15,15 +15,28 @@
 
 
 def agglomerate(
-    affs: NDArray[np.float32],
-    thresholds: Sequence[float],
+    affs: NDArray[np.float32] | None = None,
+    thresholds: Sequence[float] | None = None,
     gt: NDArray[np.uint32] | None = None,
     fragments: NDArray[np.uint64] | None = None,
+    semantic: NDArray[np.uint8] | None = None,
+    segconstraint: NDArray[np.uint64] | None = None,
+    input_rag=None,
+    input_rag_metadata=None,
     aff_threshold_low: float = 0.0001,
     aff_threshold_high: float = 0.9999,
     return_merge_history: bool = False,
     return_region_graph: bool = False,
+    return_region_graph_metadata: bool = False,
     scoring_function: str = "OneMinus<MeanAffinity<RegionGraphType, ScoreValue>>",
+    semantic_aff_threshold: float = 0.5,
+    semantic_size_threshold: int = 100_000,
+    semantic_signal_ratio: float = 0.6,
+    semantic_taint_labels: list[int] = [],
+    semantic_taint_threshold: float = 0.0,
+    size_heuristic_aff_threshold: float = 1.0,
+    size_heuristic_small_threshold: int = 1_000_000,
+    size_heuristic_large_threshold: int = 10_000_000,
     discretize_queue: int = 0,
     force_rebuild: bool = False,
 ) -> Iterator[tuple | NDArray[np.uint64]]:
@@ -142,48 +155,88 @@ def agglomerate(
             affs, range(100,10000,100), gt, return_merge_history = True):
             # ...
     """
-    import witty
-
-    with TemporaryDirectory() as tmpdir:
-        # supply #include <ScoringFunction.h> in frontend_agglomerate.h
-        tmp_path = Path(tmpdir)
-        scoredef = f"typedef {scoring_function} ScoringFunctionType;"
-        (tmp_path / "ScoringFunction.h").write_text(scoredef)
-
-        # supply #include <Queue.h> in frontend_agglomerate.h
-        queue_src = "template<typename T, typename S> using QueueType = " + (
-            "PriorityQueue<T, S>;"
-            if discretize_queue == 0
-            else f"BinQueue<T, S, {discretize_queue}>;"
-        )
-        (tmp_path / "Queue.h").write_text(queue_src)
-
-        # compile module
-        module = witty.compile_cython(
-            (HERE / "agglomerate.pyx").read_text(),
-            source_files=[str(HERE / "frontend_agglomerate.cpp")],
-            extra_link_args=["-std=c++11"],
-            extra_compile_args=["-std=c++11", "-w"],
-            include_dirs=[
+    _DEFAULT_SCORING = "OneMinus<MeanAffinity<RegionGraphType, ScoreValue>>"
+    _use_precompiled = (scoring_function == _DEFAULT_SCORING and discretize_queue == 0)
+
+    if _use_precompiled:
+        from waterz import _agglomerate_default as module
+    else:
+        import witty
+
+        with TemporaryDirectory() as tmpdir:
+            tmp_path = Path(tmpdir)
+            scoredef = f"typedef {scoring_function} ScoringFunctionType;"
+            (tmp_path / "ScoringFunction.h").write_text(scoredef)
+
+            queue_src = "template<typename T, typename S> using QueueType = " + (
+                "PriorityQueue<T, S>;"
+                if discretize_queue == 0
+                else f"BinQueue<T, S, {discretize_queue}>;"
+            )
+            (tmp_path / "Queue.h").write_text(queue_src)
+
+            _include_dirs = [
                 str(HERE),
                 tmpdir,
                 str(HERE / "backend"),
                 np.get_include(),
                 "/opt/homebrew/include",
-            ],
-            language="c++",
-            quiet=True,
-            force_rebuild=force_rebuild,
-        )
+            ]
+            _compile_args = ["-std=c++11", "-w"]
+
+            def _build_frontend(cache_dir: Path) -> list[str]:
+                import subprocess
+
+                obj_path = cache_dir / "frontend_agglomerate.o"
+                cpp_path = HERE / "frontend_agglomerate.cpp"
+                if not obj_path.exists() or obj_path.stat().st_mtime < cpp_path.stat().st_mtime:
+                    cmd = [
+                        "c++", *_compile_args,
+                        *[f"-I{d}" for d in _include_dirs],
+                        "-fPIC", "-c", str(cpp_path), "-o", str(obj_path),
+                    ]
+                    subprocess.check_call(cmd)
+                return [str(obj_path)]
+
+            module = witty.compile_cython(
+                (HERE / "agglomerate.pyx").read_text(),
+                source_files=[str(HERE / "frontend_agglomerate.cpp")],
+                build_extra_objects=_build_frontend,
+                extra_link_args=["-std=c++11"],
+                extra_compile_args=_compile_args,
+                include_dirs=_include_dirs,
+                language="c++",
+                quiet=True,
+                force_rebuild=force_rebuild,
+            )
 
     # call compiled function
+    if input_rag is not None or input_rag_metadata is not None:
+        return module.agglomerate_rag(
+            rag=input_rag,
+            rag_metadata=input_rag_metadata,
+            thresholds=thresholds,
+            fragments=fragments,
+        )
+
     return module.agglomerate(
-        affs,
-        thresholds,
-        gt,
-        fragments,
-        aff_threshold_low,
-        aff_threshold_high,
-        return_merge_history,
-        return_region_graph,
+        affs=affs,
+        thresholds=thresholds,
+        gt=gt,
+        fragments=fragments,
+        semantic=semantic,
+        segconstraint=segconstraint,
+        aff_threshold_low=aff_threshold_low,
+        aff_threshold_high=aff_threshold_high,
+        semantic_aff_threshold=semantic_aff_threshold,
+        semantic_size_threshold=semantic_size_threshold,
+        semantic_signal_ratio=semantic_signal_ratio,
+        semantic_taint_labels=semantic_taint_labels,
+        semantic_taint_threshold=semantic_taint_threshold,
+        size_heuristic_aff_threshold=size_heuristic_aff_threshold,
+        size_heuristic_small_threshold=size_heuristic_small_threshold,
+        size_heuristic_large_threshold=size_heuristic_large_threshold,
+        return_merge_history=return_merge_history,
+        return_region_graph=return_region_graph,
+        return_region_graph_metadata=return_region_graph_metadata,
     )
diff --git a/src/waterz/agglomerate.pyx b/src/waterz/agglomerate.pyx
index 07a3583..218a8a3 100644
--- a/src/waterz/agglomerate.pyx
+++ b/src/waterz/agglomerate.pyx
@@ -1,22 +1,62 @@
 from libcpp.vector cimport vector
-from libc.stdint cimport uint64_t, uint32_t
+from libc.stdint cimport uint64_t, uint32_t, uint8_t
 from libcpp cimport bool
+from libcpp.string cimport string
+
 import numpy as np
 cimport numpy as np
 
-def agglomerate(
-        affs,
+def agglomerate_rag(
+        rag,
+        rag_metadata,
         thresholds,
-        gt=None,
-        fragments=None,
-        aff_threshold_low=0.0001,
-        aff_threshold_high=0.9999,
-        return_merge_history=False,
-        return_region_graph=False):
+        fragments,
+        ):
+
+    if fragments is not None and not fragments.flags['C_CONTIGUOUS']:
+        print("Creating memory-contiguous fragments arrray (avoid this by passing C_CONTIGUOUS arrays)")
+        fragments = np.ascontiguousarray(fragments)
+
+    cdef WaterzState state = __initialize_with_rag(rag, rag_metadata, fragments)
+
+    thresholds.sort()
+    for threshold in thresholds:
+        merge_history = mergeUntil(state, threshold)
+        yield merge_history, fragments
+
+    free(state)
+
+
+def agglomerate(
+        affs: np.array | None,
+        thresholds: np.array | None,
+        gt: np.array | None,
+        fragments: np.array | None,
+        semantic: np.array | None,
+        segconstraint: np.array | None,
+
+        aff_threshold_low: float,
+        aff_threshold_high: float,
+
+        semantic_aff_threshold: float,
+        semantic_size_threshold: int,
+        semantic_signal_ratio: float,
+
+        semantic_taint_labels: list,
+        semantic_taint_threshold: float,
+
+        size_heuristic_aff_threshold: float,
+        size_heuristic_small_threshold: int,
+        size_heuristic_large_threshold: int,
+
+        return_merge_history: bool,
+        return_region_graph: bool,
+        return_region_graph_metadata: bool,
+    ):
 
     # the C++ part assumes contiguous memory, make sure we have it (and do 
     # nothing, if we do)
-    if not affs.flags['C_CONTIGUOUS']:
+    if affs is not None and not affs.flags['C_CONTIGUOUS']:
         print("Creating memory-contiguous affinity arrray (avoid this by passing C_CONTIGUOUS arrays)")
         affs = np.ascontiguousarray(affs)
     if gt is not None and not gt.flags['C_CONTIGUOUS']:
@@ -25,6 +65,12 @@ def agglomerate(
     if fragments is not None and not fragments.flags['C_CONTIGUOUS']:
         print("Creating memory-contiguous fragments arrray (avoid this by passing C_CONTIGUOUS arrays)")
         fragments = np.ascontiguousarray(fragments)
+    if semantic is not None and not semantic.flags['C_CONTIGUOUS']:
+        print("Creating memory-contiguous semantic arrray (avoid this by passing C_CONTIGUOUS arrays)")
+        semantic = np.ascontiguousarray(semantic)
+    if segconstraint is not None and not segconstraint.flags['C_CONTIGUOUS']:
+        print("Creating memory-contiguous segconstraint arrray (avoid this by passing C_CONTIGUOUS arrays)")
+        segconstraint = np.ascontiguousarray(segconstraint)
 
     print("Preparing segmentation volume...")
 
@@ -36,7 +82,24 @@ def agglomerate(
         segmentation = fragments
         find_fragments = False
 
-    cdef WaterzState state = __initialize(affs, segmentation, gt, aff_threshold_low, aff_threshold_high, find_fragments)
+    cdef WaterzState state = __initialize(
+        affs=affs,
+        segmentation=segmentation,
+        gt=gt,
+        semantic=semantic,
+        segconstraint=segconstraint,
+        aff_threshold_low=aff_threshold_low,
+        aff_threshold_high=aff_threshold_high,
+        semantic_aff_threshold=semantic_aff_threshold,
+        semantic_size_threshold=semantic_size_threshold,
+        semantic_signal_ratio=semantic_signal_ratio,
+        semantic_taint_labels=semantic_taint_labels,
+        semantic_taint_threshold=semantic_taint_threshold,
+        size_heuristic_aff_threshold=size_heuristic_aff_threshold,
+        size_heuristic_small_threshold=size_heuristic_small_threshold,
+        size_heuristic_large_threshold=size_heuristic_large_threshold,
+        find_fragments=find_fragments,
+        )
 
     thresholds.sort()
     for threshold in thresholds:
@@ -63,6 +126,10 @@ def agglomerate(
 
             result += (getRegionGraph(state),)
 
+        if return_region_graph_metadata:
+
+            result += (getRegionGraphMeta(state),)
+
         if len(result) == 1:
             yield result[0]
         else:
@@ -72,29 +139,95 @@ def agglomerate(
 
 def __initialize(
         np.ndarray[np.float32_t, ndim=4] affs,
-        np.ndarray[uint64_t, ndim=3]     segmentation,
-        np.ndarray[uint32_t, ndim=3]     gt = None,
-        aff_threshold_low  = 0.0001,
-        aff_threshold_high = 0.9999,
-        find_fragments = True):
+        np.ndarray[uint64_t, ndim=3] segmentation,
+        np.ndarray[uint32_t, ndim=3] gt,
+        np.ndarray[uint8_t, ndim=3] semantic,
+        np.ndarray[uint64_t, ndim=3] segconstraint,
+        aff_threshold_low: float,
+        aff_threshold_high: float,
+        semantic_aff_threshold: float,
+        semantic_size_threshold: int,
+        semantic_signal_ratio: float,
+        semantic_taint_labels: list,
+        semantic_taint_threshold: float,
+        size_heuristic_aff_threshold: float,
+        size_heuristic_small_threshold: int,
+        size_heuristic_large_threshold: int,
+        find_fragments: bool,
+    ):
 
     cdef float*    aff_data
     cdef uint64_t* segmentation_data
     cdef uint32_t* gt_data = NULL
+    cdef uint8_t* semantic_data = NULL
+    cdef uint64_t* segconstraint_data = NULL
+    cdef vector[uint8_t] taint_labels_vec
+    for label in semantic_taint_labels:
+        taint_labels_vec.push_back(<uint8_t>label)
 
     aff_data = &affs[0,0,0,0]
     segmentation_data = &segmentation[0,0,0]
     if gt is not None:
         gt_data = &gt[0,0,0]
+    if semantic is not None:
+        semantic_data = &semantic[0,0,0]
+    if segconstraint is not None:
+        segconstraint_data = &segconstraint[0,0,0]
+
+    # return initialize(
+    #     affs.shape[1], affs.shape[2], affs.shape[3],
+    #     aff_data,
+    #     segmentation_data,
+    #     gt_data,
+    #     aff_threshold_low,
+    #     aff_threshold_high,
+    #     find_fragments
 
     return initialize(
-        affs.shape[1], affs.shape[2], affs.shape[3],
-        aff_data,
+        width=affs.shape[1],
+        height=affs.shape[2],
+        depth=affs.shape[3],
+        affinity_data=aff_data,
+        segmentation_data=segmentation_data,
+        groundtruth_data=gt_data,
+        semantic_data=semantic_data,
+        segconstraint_data=segconstraint_data,
+
+        affThresholdLow=aff_threshold_low,
+        affThresholdHigh=aff_threshold_high,
+
+        semantic_aff_threshold=semantic_aff_threshold,
+        semantic_size_threshold=semantic_size_threshold,
+        semantic_signal_ratio=semantic_signal_ratio,
+
+        semantic_taint_labels=taint_labels_vec,
+        semantic_taint_threshold=semantic_taint_threshold,
+
+        size_heuristic_aff_threshold=size_heuristic_aff_threshold,
+        size_heuristic_small_threshold=size_heuristic_small_threshold,
+        size_heuristic_large_threshold=size_heuristic_large_threshold,
+
+        findFragments=find_fragments,
+    )
+
+def __initialize_with_rag(
+        rag,
+        rag_metadata,
+        np.ndarray[uint64_t, ndim=3]     segmentation = None):
+
+    cdef uint64_t* segmentation_data = NULL
+    shape = (0, 0, 0)
+
+    if segmentation is not None:
+        segmentation_data = &segmentation[0,0,0]
+        shape = (segmentation.shape[0], segmentation.shape[1], segmentation.shape[2])
+
+
+    return initialize_with_rag(
+        rag,
+        rag_metadata,
         segmentation_data,
-        gt_data,
-        aff_threshold_low,
-        aff_threshold_high,
-        find_fragments)
+        shape[0], shape[1], shape[2])
 
 cdef extern from "frontend_agglomerate.h":
 
@@ -126,14 +259,34 @@ cdef extern from "frontend_agglomerate.h":
             const float*    affinity_data,
             uint64_t*       segmentation_data,
             const uint32_t* groundtruth_data,
+            const uint8_t*  semantic_data,
+            const uint64_t* segconstraint_data,
             float           affThresholdLow,
             float           affThresholdHigh,
+            float           semantic_aff_threshold,
+            uint64_t        semantic_size_threshold,
+            float           semantic_signal_ratio,
+            const vector[uint8_t]& semantic_taint_labels,
+            float           semantic_taint_threshold,
+            float           size_heuristic_aff_threshold,
+            uint64_t        size_heuristic_small_threshold,
+            uint64_t        size_heuristic_large_threshold,
             bool            findFragments);
 
+    WaterzState initialize_with_rag(
+            const vector[ScoredEdge]& rag,
+            const vector[double]& rag_metadata,
+            uint64_t*       segmentation_data,
+            size_t          width,
+            size_t          height,
+            size_t          depth,);
+
     vector[Merge] mergeUntil(
             WaterzState& state,
             float        threshold)
 
     vector[ScoredEdge] getRegionGraph(WaterzState& state)
 
+    vector[double] getRegionGraphMeta(WaterzState& state)
+
     void free(WaterzState& state)
diff --git a/src/waterz/backend/ConstraintProvider.hpp b/src/waterz/backend/ConstraintProvider.hpp
new file mode 100644
index 0000000..a2cbd09
--- /dev/null
+++ b/src/waterz/backend/ConstraintProvider.hpp
@@ -0,0 +1,13 @@
+#ifndef WATERZ_CONSTRAINT_PROVIDER_H__
+#define WATERZ_CONSTRAINT_PROVIDER_H__
+
+/**
+ * Base class for statistics providers with fallback implementations.
+ */
+class ConstraintProvider {
+public:
+	virtual inline bool notifyNodeMerge(uint64_t from, uint64_t to) = 0;
+	virtual inline bool isConstrained(uint64_t from, uint64_t to, float score) const = 0;
+};
+
+#endif // WATERZ_CONSTRAINT_PROVIDER_H__
diff --git a/src/waterz/backend/DefaultDict.hpp b/src/waterz/backend/DefaultDict.hpp
new file mode 100644
index 0000000..4faa08e
--- /dev/null
+++ b/src/waterz/backend/DefaultDict.hpp
@@ -0,0 +1,70 @@
+#ifndef WATERZ_DEFAULTDICT_H__
+#define WATERZ_DEFAULTDICT_H__
+
+#include <unordered_map>
+
+template<typename T>
+uint64_t getDictSum(const T& map) {
+    uint64_t sum = 0;
+    for (const auto& pair : map) {
+        sum += pair.second;
+    }
+    return sum;
+}
+
+template<typename T>
+auto getDictMaxKey(const T& map) -> decltype(map.begin()->first) {
+    if (map.empty()) {
+        throw std::runtime_error("Map is empty");
+    }
+    auto max_key = map.begin()->first;
+    auto max_value = map.begin()->second;
+    for (const auto& pair : map) {
+        if (pair.second > max_value) {
+            max_key = pair.first;
+            max_value = pair.second;
+        }
+    }
+    return max_key;
+}
+
+template<typename K, typename V>
+class DefaultDict {
+private:
+    std::unordered_map<K, V> container;
+    V default_value;
+
+public:
+    DefaultDict() : default_value() {}
+    DefaultDict(const V& default_val) : default_value(default_val) {}
+
+    V& operator[](const K& key) {
+        auto result = container.emplace(key, default_value);
+        return result.first->second;
+    }
+
+    V operator[](const K& key) const {
+        auto it = container.find(key);
+        if (it != container.end()) {
+            return it->second;
+        }
+        return default_value;
+    }
+
+    // V at(const K& key) const {
+    //     auto it = container.find(key);
+    //     if (it != container.end()) {
+    //         return it->second;
+    //     }
+    //     return default_value;
+    // }
+
+    void erase(const K& key) {
+        container.erase(key);
+    }
+
+    const std::unordered_map<K, V>& getContainer() const { return container; }
+};
+
+
+#endif // WATERZ_DEFAULTDICT_H__
diff --git a/src/waterz/backend/IterativeRegionMerging.hpp b/src/waterz/backend/IterativeRegionMerging.hpp
index 9689714..513a0e3 100644
--- a/src/waterz/backend/IterativeRegionMerging.hpp
+++ b/src/waterz/backend/IterativeRegionMerging.hpp
@@ -4,12 +4,14 @@
 #include <iostream>
 #include <vector>
 #include <map>
+#include <unordered_map>
 #include <queue>
 #include <cassert>
 #include <limits>
 
 #include "RegionGraph.hpp"
 #include "PriorityQueue.hpp"
+#include "ConstraintProvider.hpp"
 
 template <typename NodeIdType, typename ScoreType, template <typename T, typename S> class QueueType = PriorityQueue>
 class IterativeRegionMerging {
@@ -37,6 +39,7 @@ class IterativeRegionMerging {
 	std::size_t mergeUntil(
 			EdgeScoringFunction& edgeScoringFunction,
 			StatisticsProviderType& statisticsProvider,
+			std::vector<ConstraintProvider*>& constraints,
 			ScoreType threshold,
 			Visitor& visitor) {
 
@@ -100,12 +103,23 @@ class IterativeRegionMerging {
 				continue;
 			}
 
-			NodeIdType newRegion = mergeRegions(next, statisticsProvider);
+			NodeIdType a = _regionGraph.edge(next).u;
+			NodeIdType b = _regionGraph.edge(next).v;
+
+			bool is_constrained = false;
+			for (auto constraint : constraints) {
+				is_constrained |= constraint->isConstrained(a, b, score);
+			}
+
+			if (is_constrained) {
+				continue;  // skip merging
+			}
+
+			NodeIdType newRegion = mergeRegions(next, statisticsProvider, constraints);
 			merged++;
 
 			visitor.onMerge(
-					_regionGraph.edge(next).u,
-					_regionGraph.edge(next).v,
+					a, b,
 					newRegion,
 					score);
 		}
@@ -163,6 +177,31 @@ class IterativeRegionMerging {
 		return edges;
 	}
 
+	template <typename EdgeMetadata, typename EdgeScoringFunction, typename StatisticsProviderType>
+	std::vector<EdgeMetadata> extractRegionGraphMeta(EdgeScoringFunction& edgeScoringFunction, StatisticsProviderType& statisticsProvider) {
+
+		std::vector<EdgeMetadata> ret;
+
+		for (EdgeIdType e = 0; e < _regionGraph.numEdges(); e++) {
+
+			if (_deleted[e])
+				continue;
+
+			ScoreType score;
+			if (_stale[e])
+				score = scoreEdge(e, edgeScoringFunction);
+			else
+				score = _edgeScores[e];
+
+			if (score < _mergedUntil)
+				continue;
+
+			ret.emplace_back(statisticsProvider.getEdgeMetadata(e));
+		}
+
+		return ret;
+	}
+
 private:
 
 	/**
@@ -171,7 +210,8 @@ class IterativeRegionMerging {
 	template <typename StatisticsProviderType>
 	NodeIdType mergeRegions(
 			EdgeIdType e,
-			StatisticsProviderType& statisticsProvider) {
+			StatisticsProviderType& statisticsProvider,
+			std::vector<ConstraintProvider*>& constraints) {
 
 		NodeIdType a = _regionGraph.edge(e).u;
 		NodeIdType b = _regionGraph.edge(e).v;
@@ -179,6 +219,10 @@ class IterativeRegionMerging {
 		// assign new node a = a + b
 		bool nodeStatisticsChanged = statisticsProvider.notifyNodeMerge(b, a);
 
+		for (auto constraint : constraints) {
+			nodeStatisticsChanged |= constraint->notifyNodeMerge(b, a);
+		}
+
 		// set path
 		_rootPaths[b] = a;
 
@@ -190,7 +234,10 @@ class IterativeRegionMerging {
 		}
 
 		// ...and update incident edges of b
-		std::vector<EdgeIdType> neighborEdges = _regionGraph.incEdges(b);
+		// Take b's incident list via swap to avoid copying the vector.
+		// b's list is cleared, so use reassignEdge/removeEdgeSkipNode
+		// which don't touch b's incEdges.
+		std::vector<EdgeIdType> neighborEdges = _regionGraph.takeIncEdges(b);
 		for (EdgeIdType neighborEdge : neighborEdges) {
 
 			if (neighborEdge == e)
@@ -209,8 +256,7 @@ class IterativeRegionMerging {
 
 				// We encountered an exclusive neighbor of b.
 
-				_regionGraph.moveEdge(neighborEdge, a, neighbor);
-				assert(_regionGraph.findEdge(a, neighbor) == neighborEdge);
+				_regionGraph.reassignEdge(neighborEdge, b, a);
 
 				if (nodeStatisticsChanged)
 					_stale[neighborEdge] = true;
@@ -223,32 +269,30 @@ class IterativeRegionMerging {
 				// * mark the cheaper one as stale (if it isn't already)
 				// * delete the more expensive one
 				//
-				// This ensures that the stale edge bubbles up early enough 
-				// to consider it's real score (which is assumed to be 
+				// This ensures that the stale edge bubbles up early enough
+				// to consider it's real score (which is assumed to be
 				// larger than the minium of the two original scores).
 
 				if (_edgeScores[neighborEdge] > _edgeScores[aNeighborEdge]) {
 
-					// We got lucky, we can reuse the edge that is attached to a 
+					// We got lucky, we can reuse the edge that is attached to a
 					// already
 
 					bool edgeStatisticChanged = statisticsProvider.notifyEdgeMerge(neighborEdge, aNeighborEdge);
 
-					_regionGraph.removeEdge(neighborEdge);
+					_regionGraph.removeEdgeSkipNode(neighborEdge, b);
 					_deleted[neighborEdge] = true;
 					if (edgeStatisticChanged)
 						_stale[aNeighborEdge] = true;
 
 				} else {
 
-					// Bummer. The new edge should be the one pointing from 
+					// Bummer. The new edge should be the one pointing from
 					// a to neighbor.
 
 					bool edgeStatisticChanged = statisticsProvider.notifyEdgeMerge(aNeighborEdge, neighborEdge);
 
-					_regionGraph.removeEdge(aNeighborEdge);
-					_regionGraph.moveEdge(neighborEdge, a, neighbor);
-					assert(_regionGraph.findEdge(a, neighbor) == neighborEdge);
+					_regionGraph.replaceEdge(aNeighborEdge, neighborEdge, a, neighbor, b);
 
 					if (edgeStatisticChanged)
 						_stale[neighborEdge] = true;
@@ -324,7 +368,7 @@ class IterativeRegionMerging {
 	// root nodes are not in the map
 	//
 	// paths will be compressed when read
-	std::map<NodeIdType, NodeIdType> _rootPaths;
+	std::unordered_map<NodeIdType, NodeIdType> _rootPaths;
 
 	// current state of merging
 	ScoreType _mergedUntil;
diff --git a/src/waterz/backend/MeanAffinityProvider.hpp b/src/waterz/backend/MeanAffinityProvider.hpp
index 5ff45a4..a8a3b3a 100644
--- a/src/waterz/backend/MeanAffinityProvider.hpp
+++ b/src/waterz/backend/MeanAffinityProvider.hpp
@@ -27,6 +27,12 @@ class MeanAffinityProvider : public StatisticsProvider {
 		_numValues[e]++;
 	}
 
+	template <typename EdgeIdType, typename ScoreType>
+	inline void addEdge(EdgeIdType e, ScoreType affinity, double value) {
+		_meanAffinities[e] = affinity;
+		_numValues[e] = value;
+	}
+
 	inline bool notifyEdgeMerge(EdgeIdType from, EdgeIdType to) {
 
 		size_t fromN = _numValues[from];
@@ -46,6 +52,10 @@ class MeanAffinityProvider : public StatisticsProvider {
 		return _meanAffinities[e];
 	}
 
+	inline double getEdgeMetadata(EdgeIdType e) const {
+		return _numValues[e];
+	}
+
 private:
 
 	typename RegionGraphType::template EdgeMap<size_t> _numValues;
diff --git a/src/waterz/backend/RegionGraph.hpp b/src/waterz/backend/RegionGraph.hpp
index 4dda0ae..ac7a540 100644
--- a/src/waterz/backend/RegionGraph.hpp
+++ b/src/waterz/backend/RegionGraph.hpp
@@ -3,6 +3,7 @@
 
 #include <algorithm>
 #include <vector>
+#include <unordered_map>
 #include <limits>
 #include <cassert>
 
@@ -164,7 +165,8 @@ class RegionGraph {
 
 	RegionGraph(ID numNodes = 0) :
 		_numNodes(numNodes),
-		_incEdges(numNodes) {}
+		_incEdges(numNodes),
+		_adjMap(numNodes) {}
 
 	ID numNodes() const { return _numNodes; }
 
@@ -175,6 +177,7 @@ class RegionGraph {
 		NodeIdType id = _numNodes;
 		_numNodes++;
 		_incEdges.emplace_back();
+		_adjMap.emplace_back();
 
 		for (RegionGraphNodeMapBase<ID>* map : _nodeMaps)
 			map->onNewNode(id);
@@ -189,6 +192,8 @@ class RegionGraph {
 
 		_incEdges[u].push_back(id);
 		_incEdges[v].push_back(id);
+		_adjMap[u][v] = id;
+		_adjMap[v][u] = id;
 
 		for (RegionGraphEdgeMapBase<ID>* map : _edgeMaps)
 			map->onNewEdge(id);
@@ -198,8 +203,12 @@ class RegionGraph {
 
 	void removeEdge(EdgeIdType e) {
 
-		removeIncEdge(_edges[e].u, e);
-		removeIncEdge(_edges[e].v, e);
+		NodeIdType u = _edges[e].u;
+		NodeIdType v = _edges[e].v;
+		removeIncEdge(u, e);
+		removeIncEdge(v, e);
+		_adjMap[u].erase(v);
+		_adjMap[v].erase(u);
 	}
 
 	void moveEdge(EdgeIdType e, NodeIdType u, NodeIdType v) {
@@ -263,32 +272,105 @@ class RegionGraph {
 
 	inline const std::vector<EdgeIdType>& incEdges(ID node) const { return _incEdges[node]; }
 
+	inline std::vector<EdgeIdType> takeIncEdges(ID node) {
+		std::vector<EdgeIdType> result;
+		result.swap(_incEdges[node]);
+		return result;
+	}
+
 	inline NodeIdType getOpposite(NodeIdType n, EdgeIdType e) const {
 
 		return (_edges[e].u == n ? _edges[e].v : _edges[e].u);
 	}
 
 	/**
-	 * Find the edge connecting u and v. Returns NoEdge, if there is none.
+	 * Fast edge reassignment: move edge e from oldNode to newNode.
+	 * Caller must ensure one endpoint of e is oldNode.
+	 * Does NOT touch oldNode's incEdges (caller handles that).
 	 */
-	inline EdgeIdType findEdge(NodeIdType u, NodeIdType v) {
+	void reassignEdge(EdgeIdType e, NodeIdType oldNode, NodeIdType newNode) {
+
+		NodeIdType other = getOpposite(oldNode, e);
+
+		// Update adjacency maps
+		_adjMap[oldNode].erase(other);
+		_adjMap[other].erase(oldNode);
+		_adjMap[newNode][other] = e;
+		_adjMap[other][newNode] = e;
+
+		// Update edge endpoints
+		if (_edges[e].u == oldNode)
+			_edges[e].u = newNode;
+		else
+			_edges[e].v = oldNode == _edges[e].v ? newNode : _edges[e].v;
+
+		// Ensure u < v
+		if (_edges[e].u > _edges[e].v)
+			std::swap(_edges[e].u, _edges[e].v);
 
-		return findEdge(u, v, (_incEdges[u].size() < _incEdges[v].size() ? _incEdges[u] : _incEdges[v]));
+		// Add to newNode's incident list (don't remove from oldNode - caller owns that)
+		_incEdges[newNode].push_back(e);
 	}
 
 	/**
-	 * Same as findEdge(u, v), but restricted to edges in pool.
+	 * Remove edge from graph but don't touch incEdges of the given skipNode.
+	 * Used when the caller already owns/cleared that node's incident list.
 	 */
-	inline EdgeIdType findEdge(NodeIdType u, NodeIdType v, const std::vector<EdgeIdType>& pool) {
+	void removeEdgeSkipNode(EdgeIdType e, NodeIdType skipNode) {
+
+		NodeIdType u = _edges[e].u;
+		NodeIdType v = _edges[e].v;
+		NodeIdType other = (u == skipNode) ? v : u;
+		removeIncEdge(other, e);
+		_adjMap[u].erase(v);
+		_adjMap[v].erase(u);
+	}
 
-		NodeIdType min = std::min(u, v);
-		NodeIdType max = std::max(u, v);
+	/**
+	 * Replace oldEdge (between survivor and neighbor) with newEdge
+	 * (being reassigned from oldNode to survivor). Fuses removeEdge +
+	 * reassignEdge to avoid redundant incEdge and adjMap operations.
+	 * Caller must have already taken oldNode's incEdges.
+	 */
+	void replaceEdge(EdgeIdType oldEdge, EdgeIdType newEdge,
+			NodeIdType survivor, NodeIdType neighbor, NodeIdType oldNode) {
+
+		// Remove oldEdge from survivor's and neighbor's incEdges
+		removeIncEdge(survivor, oldEdge);
+		removeIncEdge(neighbor, oldEdge);
+
+		// Remove newEdge from neighbor's incEdges (oldNode's already taken)
+		removeIncEdge(neighbor, newEdge);
+
+		// Update newEdge endpoints: oldNode -> survivor
+		if (_edges[newEdge].u == oldNode)
+			_edges[newEdge].u = survivor;
+		else
+			_edges[newEdge].v = survivor;
+		if (_edges[newEdge].u > _edges[newEdge].v)
+			std::swap(_edges[newEdge].u, _edges[newEdge].v);
+
+		// Add newEdge to survivor's and neighbor's incEdges
+		_incEdges[survivor].push_back(newEdge);
+		_incEdges[neighbor].push_back(newEdge);
+
+		// Update adjMaps: survivor↔neighbor now points to newEdge
+		_adjMap[survivor][neighbor] = newEdge;
+		_adjMap[neighbor][survivor] = newEdge;
+
+		// Clean up oldNode's adjMap entries
+		_adjMap[oldNode].erase(neighbor);
+		_adjMap[neighbor].erase(oldNode);
+	}
 
-		for (EdgeIdType e : pool)
-			if (std::min(_edges[e].u, _edges[e].v) == min &&
-				std::max(_edges[e].u, _edges[e].v) == max)
-				return e;
+	/**
+	 * Find the edge connecting u and v. Returns NoEdge, if there is none.
+	 */
+	inline EdgeIdType findEdge(NodeIdType u, NodeIdType v) {
 
+		auto it = _adjMap[u].find(v);
+		if (it != _adjMap[u].end())
+			return it->second;
 		return NoEdge;
 	}
 
@@ -323,15 +405,27 @@ class RegionGraph {
 
 	inline void moveEdgeNodeV(EdgeIdType e, NodeIdType v) {
 
-		removeIncEdge(_edges[e].v, e);
+		NodeIdType oldV = _edges[e].v;
+		NodeIdType otherNode = _edges[e].u;
+		removeIncEdge(oldV, e);
+		_adjMap[oldV].erase(otherNode);
+		_adjMap[otherNode].erase(oldV);
 		_incEdges[v].push_back(e);
+		_adjMap[v][otherNode] = e;
+		_adjMap[otherNode][v] = e;
 		_edges[e].v = v;
 	}
 
 	inline void moveEdgeNodeU(EdgeIdType e, NodeIdType u) {
 
-		removeIncEdge(_edges[e].u, e);
+		NodeIdType oldU = _edges[e].u;
+		NodeIdType otherNode = _edges[e].v;
+		removeIncEdge(oldU, e);
+		_adjMap[oldU].erase(otherNode);
+		_adjMap[otherNode].erase(oldU);
 		_incEdges[u].push_back(e);
+		_adjMap[u][otherNode] = e;
+		_adjMap[otherNode][u] = e;
 		_edges[e].u = u;
 	}
 
@@ -339,7 +433,8 @@ class RegionGraph {
 
 		auto it = std::find(_incEdges[n].begin(), _incEdges[n].end(), e);
 		assert(it != _incEdges[n].end());
-		_incEdges[n].erase(it);
+		std::swap(*it, _incEdges[n].back());
+		_incEdges[n].pop_back();
 		assert(std::find(_incEdges[n].begin(), _incEdges[n].end(), e) == _incEdges[n].end());
 	}
 
@@ -349,6 +444,9 @@ class RegionGraph {
 
 	std::vector<std::vector<EdgeIdType>> _incEdges;
 
+	// per-node adjacency map: neighbor -> edge id (O(1) findEdge)
+	std::vector<std::unordered_map<NodeIdType, EdgeIdType>> _adjMap;
+
 	std::vector<RegionGraphNodeMapBase<ID>*> _nodeMaps;
 	std::vector<RegionGraphEdgeMapBase<ID>*> _edgeMaps;
 };
diff --git a/src/waterz/backend/SegConstraintProvider.hpp b/src/waterz/backend/SegConstraintProvider.hpp
new file mode 100644
index 0000000..d28a724
--- /dev/null
+++ b/src/waterz/backend/SegConstraintProvider.hpp
@@ -0,0 +1,52 @@
+#include "ConstraintProvider.hpp"
+#include "DefaultDict.hpp"
+
+#include <unordered_map>
+
+template <typename RegionGraphType, typename SegType>
+class SegConstraintProvider : public ConstraintProvider {
+
+    typedef typename RegionGraphType::NodeIdType NodeIdType;
+
+private:
+    std::unordered_map<SegType, SegType> _constraint;
+
+public:
+	SegConstraintProvider(
+            const SegType *segconstraint_data,
+            const SegType *seg_data,
+            size_t num_voxels) {
+
+        // Collect all mapped constraint, accounting for when an object
+        // is straddled across multiple constraint
+        DefaultDict<SegType, DefaultDict<SegType, size_t>> constraint{
+            DefaultDict<SegType, size_t>{0}
+        };
+		for (std::size_t i = 0; i < num_voxels; i++) {
+            constraint[seg_data[i]][segconstraint_data[i]] += 1;
+        }
+
+        // Now we get a single dominant constraint per object
+        for (const auto& pair : constraint.getContainer()) {
+            const auto& segid = pair.first;
+            const auto& mapped_constraint_ids = pair.second;
+            _constraint[segid] = getDictMaxKey(mapped_constraint_ids.getContainer());
+        }
+    }
+
+    inline bool notifyNodeMerge(NodeIdType from, NodeIdType to) {
+        if (_constraint.at(to) == 0) {
+            _constraint[to] = _constraint.at(from);
+        }
+        _constraint.erase(from);
+        return true;
+    }
+
+    inline bool isConstrained(NodeIdType from, NodeIdType to, float score) const {
+        if (_constraint.at(from) == 0 || _constraint.at(to) == 0)
+            return false;
+        if (_constraint.at(from) == _constraint.at(to))
+            return false;
+        return true;
+    }
+};
diff --git a/src/waterz/backend/SemanticConstraintProvider.hpp b/src/waterz/backend/SemanticConstraintProvider.hpp
new file mode 100644
index 0000000..885e69e
--- /dev/null
+++ b/src/waterz/backend/SemanticConstraintProvider.hpp
@@ -0,0 +1,74 @@
+#include "ConstraintProvider.hpp"
+#include "DefaultDict.hpp"
+
+#include <unordered_map>
+
+template <typename RegionGraphType, typename SemValue, typename SegType>
+class SemanticConstraintProvider : public ConstraintProvider {
+
+	typedef SemValue ValueType;
+	typedef typename RegionGraphType::EdgeIdType EdgeIdType;
+    typedef typename RegionGraphType::NodeIdType NodeIdType;
+
+private:
+    DefaultDict<SegType, DefaultDict<SemValue, size_t>> _semantic{
+        DefaultDict<SemValue, size_t>{0}
+    };
+    float semantic_aff_threshold;
+    uint64_t semantic_size_threshold;
+    float semantic_signal_ratio;
+
+public:
+	SemanticConstraintProvider(
+        const SemValue *semantic_data,
+        const SegType *seg_data,
+        size_t num_voxels,
+        float semantic_aff_threshold,
+        uint64_t semantic_size_threshold,
+        float semantic_signal_ratio
+    ) :
+    semantic_aff_threshold(semantic_aff_threshold),
+    semantic_size_threshold(semantic_size_threshold),
+    semantic_signal_ratio(semantic_signal_ratio)
+    {
+		for (std::size_t i = 0; i < num_voxels; i++) {
+            _semantic[seg_data[i]][semantic_data[i]] += 1;
+        }
+    }
+
+    inline bool notifyNodeMerge(NodeIdType from, NodeIdType to) {
+        for (auto k : _semantic[from].getContainer()) {
+            // _semantic[to][k] += _semantic[from][k];
+            _semantic[to][k.first] += k.second;
+        }
+        _semantic.erase(from);
+        return true;
+    }
+
+    inline bool isConstrained(NodeIdType from, NodeIdType to, float score) const {
+
+        if (score < semantic_aff_threshold)
+            return false;
+
+        auto max_sem1_label = getDictMaxKey(_semantic[from].getContainer());
+        auto max_sem1 = _semantic[from][max_sem1_label];
+        auto total_sem1 = getDictSum(_semantic[from].getContainer());
+        auto max_sem2_label = getDictMaxKey(_semantic[to].getContainer());
+        auto max_sem2 = _semantic[to][max_sem2_label];
+        auto total_sem2 = getDictSum(_semantic[to].getContainer());
+
+        if (total_sem1 < semantic_size_threshold)
+            return false;
+        if (total_sem2 < semantic_size_threshold)
+            return false;
+        if (max_sem1 < semantic_signal_ratio * total_sem1)
+            return false;
+        if (max_sem2 < semantic_signal_ratio * total_sem2)
+            return false;
+        if (max_sem1_label == 0 || max_sem2_label == 0)
+            return false;
+        if (max_sem1_label == max_sem2_label)
+            return false;
+        return true;
+    }
+};
diff --git a/src/waterz/backend/SemanticTaintConstraintProvider.hpp b/src/waterz/backend/SemanticTaintConstraintProvider.hpp
new file mode 100644
index 0000000..2bcabb7
--- /dev/null
+++ b/src/waterz/backend/SemanticTaintConstraintProvider.hpp
@@ -0,0 +1,57 @@
+#ifndef WATERZ_SEMANTIC_TAINT_CONSTRAINT_PROVIDER_H__
+#define WATERZ_SEMANTIC_TAINT_CONSTRAINT_PROVIDER_H__
+
+#include "ConstraintProvider.hpp"
+#include "DefaultDict.hpp"
+
+#include <unordered_set>
+#include <vector>
+
+template <typename RegionGraphType, typename SemValue, typename SegType>
+class SemanticTaintConstraintProvider : public ConstraintProvider {
+
+    typedef typename RegionGraphType::NodeIdType NodeIdType;
+
+private:
+    DefaultDict<SegType, size_t> _taint_counts{0};
+    DefaultDict<SegType, size_t> _total_counts{0};
+    float _threshold;
+
+public:
+    SemanticTaintConstraintProvider(
+        const SemValue* semantic_data,
+        const SegType* seg_data,
+        size_t num_voxels,
+        const std::vector<SemValue>& taint_labels,
+        float threshold
+    ) : _threshold(threshold) {
+        std::unordered_set<SemValue> taint_set(taint_labels.begin(), taint_labels.end());
+        for (size_t i = 0; i < num_voxels; i++) {
+            _total_counts[seg_data[i]] += 1;
+            if (taint_set.count(semantic_data[i])) {
+                _taint_counts[seg_data[i]] += 1;
+            }
+        }
+    }
+
+    inline bool notifyNodeMerge(NodeIdType from, NodeIdType to) override {
+        _taint_counts[to] += _taint_counts[from];
+        _taint_counts.erase(from);
+        _total_counts[to] += _total_counts[from];
+        _total_counts.erase(from);
+        return true;
+    }
+
+    inline bool isConstrained(NodeIdType from, NodeIdType to, float score) const override {
+        size_t total_from = _total_counts[from];
+        size_t total_to = _total_counts[to];
+        bool from_tainted = total_from > 0 && (float)_taint_counts[from] / total_from > _threshold;
+        bool to_tainted = total_to > 0 && (float)_taint_counts[to] / total_to > _threshold;
+        // tainted segments can only merge with other tainted segments
+        if (from_tainted != to_tainted)
+            return true;
+        return false;
+    }
+};
+
+#endif // WATERZ_SEMANTIC_TAINT_CONSTRAINT_PROVIDER_H__
diff --git a/src/waterz/backend/SizeHeuristicConstraintProvider.hpp b/src/waterz/backend/SizeHeuristicConstraintProvider.hpp
new file mode 100644
index 0000000..2fd6011
--- /dev/null
+++ b/src/waterz/backend/SizeHeuristicConstraintProvider.hpp
@@ -0,0 +1,49 @@
+#include "StatisticsProvider.hpp"
+#include "DefaultDict.hpp"
+
+template <typename RegionGraphType, typename SegType>
+class SizeHeuristicConstraintProvider : public ConstraintProvider {
+
+    typedef typename RegionGraphType::NodeIdType NodeIdType;
+
+private:
+    DefaultDict<SegType, size_t> _size{0};
+    float size_heuristic_aff_threshold;
+    size_t size_heuristic_small_threshold;
+    size_t size_heuristic_large_threshold;
+
+public:
+	SizeHeuristicConstraintProvider(
+            const std::vector<size_t>& precomputed_sizes,
+            float size_heuristic_aff_threshold,
+            size_t size_heuristic_small_threshold,
+            size_t size_heuristic_large_threshold
+    ):
+    size_heuristic_aff_threshold(size_heuristic_aff_threshold),
+    size_heuristic_small_threshold(size_heuristic_small_threshold),
+    size_heuristic_large_threshold(size_heuristic_large_threshold)
+    {
+		for (std::size_t i = 0; i < precomputed_sizes.size(); i++) {
+            if (precomputed_sizes[i] > 0)
+                _size[static_cast<SegType>(i)] = precomputed_sizes[i];
+        }
+    }
+
+    inline bool notifyNodeMerge(NodeIdType from, NodeIdType to) override {
+        _size[to] += _size[from];
+        _size.erase(from);
+        return true;  // statistics changed
+    }
+
+    inline bool isConstrained(NodeIdType from, NodeIdType to, float score) const override {
+        if (score < size_heuristic_aff_threshold)
+            return false;
+        if (_size[from] < size_heuristic_small_threshold)
+            return false;
+        if (_size[to] < size_heuristic_small_threshold)
+            return false;
+        if ((_size[from] + _size[to]) < size_heuristic_large_threshold)
+            return false;
+        return true;
+    }
+};
diff --git a/src/waterz/backend/StatisticsProvider.hpp b/src/waterz/backend/StatisticsProvider.hpp
index fe7ec82..96a040f 100644
--- a/src/waterz/backend/StatisticsProvider.hpp
+++ b/src/waterz/backend/StatisticsProvider.hpp
@@ -24,6 +24,9 @@ class StatisticsProvider {
 	template <typename NodeIdType>
 	inline void addVoxel(NodeIdType n, std::size_t x, std::size_t y, std::size_t z) {}
 
+	template <typename EdgeIdType, typename ScoreType>
+	inline void addAffinity(EdgeIdType e, ScoreType affinity, double metadata) {}
+
 	/**
 	 * Callback for node merges: 'from' will be merged into 'to'. Return true, 
 	 * if this changed the statistics of this provider.
@@ -37,6 +40,9 @@ class StatisticsProvider {
 	 */
 	template<typename EdgeIdType>
 	inline bool notifyEdgeMerge(EdgeIdType from, EdgeIdType to) { return false; }
+
+	template<typename EdgeIdType>
+	inline double getEdgeMetadata(EdgeIdType e) const { return 0; }
 };
 
 #endif // WATERZ_STATISTICS_PROVIDER_H__
diff --git a/src/waterz/backend/region_graph.hpp b/src/waterz/backend/region_graph.hpp
index 4af7fae..656db6e 100644
--- a/src/waterz/backend/region_graph.hpp
+++ b/src/waterz/backend/region_graph.hpp
@@ -4,7 +4,6 @@
 
 #include <cstddef>
 #include <iostream>
-#include <map>
 
 /**
  * Extract the region graph from a segmentation. Edges are annotated with the 
@@ -40,45 +39,80 @@ get_region_graph(
 	std::ptrdiff_t ydim = aff.shape()[2];
 	std::ptrdiff_t xdim = aff.shape()[3];
 
-	// list of affinities between pairs of regions
-	std::vector<std::map<ID, std::vector<F>>> affinities(max_segid+1);
+	// Use raw pointers for direct memory access instead of boost operator[][]
+	const ID* seg_data = seg.data();
+	const F* aff_data = aff.data();
+	const std::size_t slice_size = ydim * xdim;
+	const std::size_t aff_channel_size = zdim * slice_size;
 
-	EdgeIdType e;
-	std::size_t p[3];
-	for (p[0] = 0; p[0] < zdim; ++p[0])
-		for (p[1] = 0; p[1] < ydim; ++p[1])
-			for (p[2] = 0; p[2] < xdim; ++p[2]) {
+	for (std::ptrdiff_t z = 0; z < zdim; ++z)
+		for (std::ptrdiff_t y = 0; y < ydim; ++y)
+			for (std::ptrdiff_t x = 0; x < xdim; ++x) {
 
-				ID id1 = seg[p[0]][p[1]][p[2]];
-				statisticsProvider.addVoxel(id1, p[2], p[1], p[0]);
+				std::size_t idx = z * slice_size + y * xdim + x;
+				ID id1 = seg_data[idx];
+				statisticsProvider.addVoxel(id1, x, y, z);
 
-				for (int d = 0; d < 3; d++) {
-
-					if (p[d] == 0)
-						continue;
-
-					ID id2 = seg[p[0]-(d==0)][p[1]-(d==1)][p[2]-(d==2)];
+				// d=0: z-affinity, neighbor at z-1
+				if (z > 0) {
+					ID id2 = seg_data[idx - slice_size];
+					if (id1 != id2) {
+						EdgeIdType e = rg.findEdge(id1, id2);
+						if (e == RegionGraphType::NoEdge) {
+							e = rg.addEdge(id1, id2);
+							statisticsProvider.notifyNewEdge(e);
+						}
+						statisticsProvider.addAffinity(e, aff_data[idx]);
+					}
+				}
 
+				// d=1: y-affinity, neighbor at y-1
+				if (y > 0) {
+					ID id2 = seg_data[idx - xdim];
 					if (id1 != id2) {
+						EdgeIdType e = rg.findEdge(id1, id2);
+						if (e == RegionGraphType::NoEdge) {
+							e = rg.addEdge(id1, id2);
+							statisticsProvider.notifyNewEdge(e);
+						}
+						statisticsProvider.addAffinity(e, aff_data[aff_channel_size + idx]);
+					}
+				}
 
-						auto mm = std::minmax(id1, id2);
-						affinities[mm.first][mm.second].push_back(aff[d][p[0]][p[1]][p[2]]);
+				// d=2: x-affinity, neighbor at x-1
+				if (x > 0) {
+					ID id2 = seg_data[idx - 1];
+					if (id1 != id2) {
+						EdgeIdType e = rg.findEdge(id1, id2);
+						if (e == RegionGraphType::NoEdge) {
+							e = rg.addEdge(id1, id2);
+							statisticsProvider.notifyNewEdge(e);
+						}
+						statisticsProvider.addAffinity(e, aff_data[2 * aff_channel_size + idx]);
 					}
 				}
 			}
 
-	for (ID id1 = 1; id1 <= max_segid; ++id1) {
-		for (const auto& p: affinities[id1]) {
+	std::cout << "Region graph number of edges: " << rg.edges().size() << std::endl;
+}
 
-			// p.first is ID
-			// p.second is list of affiliated edges
-			EdgeIdType e = rg.addEdge(id1, p.first);
-			statisticsProvider.notifyNewEdge(e);
+template<typename ID, typename StatisticsProviderType>
+inline
+void
+initialize_with_region_graph(
+		StatisticsProviderType& statisticsProvider,
+		RegionGraph<ID>& rg,
+		const std::vector<ScoredEdge>& edges,
+		const std::vector<double>& edges_metadata) {
 
-			for (F affinity : p.second)
-				statisticsProvider.addAffinity(e, affinity);
-        }
-    }
+	typedef RegionGraph<ID> RegionGraphType;
+	typedef typename RegionGraphType::EdgeIdType EdgeIdType;
 
-	std::cout << "Region graph number of edges: " << rg.edges().size() << std::endl;
+	for (int i = 0; i < edges.size(); ++i) {
+		const auto& edge = edges[i];
+		double size = edges_metadata[i];
+		auto aff = 1.0 - edge.score;
+		EdgeIdType e = rg.addEdge(edge.u, edge.v);
+		statisticsProvider.addEdge(e, aff, size);
+	}
 }
diff --git a/src/waterz/frontend_agglomerate.cpp b/src/waterz/frontend_agglomerate.cpp
index 100bcd6..ebb3205 100644
--- a/src/waterz/frontend_agglomerate.cpp
+++ b/src/waterz/frontend_agglomerate.cpp
@@ -10,9 +10,90 @@
 #include "backend/basic_watershed.hpp"
 #include "backend/region_graph.hpp"
 
+#include "backend/SemanticConstraintProvider.hpp"
+#include "backend/SemanticTaintConstraintProvider.hpp"
+#include "backend/SegConstraintProvider.hpp"
+#include "backend/SizeHeuristicConstraintProvider.hpp"
+
 std::map<int, WaterzContext*> WaterzContext::_contexts;
 int WaterzContext::_nextId = 0;
 
+std::vector<ScoredEdge> getRegionGraph(WaterzState& state);
+std::vector<double> getRegionGraphMeta(WaterzState& state);
+
+WaterzState
+initialize_with_rag(
+		const std::vector<ScoredEdge>& rag,
+		const std::vector<double>& rag_metadata,
+		SegID* segmentation_data,
+		std::size_t     width,
+		std::size_t     height,
+		std::size_t     depth) {
+
+	std::size_t maxId = 0;
+	for (const auto& edge : rag) {
+		auto max_uv = std::max(edge.u, edge.v);
+		maxId = std::max(max_uv, maxId);
+	}
+
+	std::size_t numNodes = maxId + 1;
+	std::cout << "creating region graph for " << numNodes << " nodes" << std::endl;
+
+	std::shared_ptr<RegionGraphType> regionGraph(
+			new RegionGraphType(numNodes)
+	);
+
+	std::cout << "creating statistics provider" << std::endl;
+	std::shared_ptr<StatisticsProviderType> statisticsProvider(
+			new StatisticsProviderType(*regionGraph)
+	);
+
+	std::cout << "initializing region graph..." << std::endl;
+
+	initialize_with_region_graph<SegID>(
+			*statisticsProvider,
+			*regionGraph,
+			rag,
+			rag_metadata);
+
+	std::shared_ptr<ScoringFunctionType> scoringFunction(
+			new ScoringFunctionType(*regionGraph, *statisticsProvider)
+	);
+
+	std::shared_ptr<RegionMergingType> regionMerging(
+			new RegionMergingType(*regionGraph)
+	);
+
+	
+	std::shared_ptr<vector<ConstraintProvider*>> dummy_constraints(
+		new vector<ConstraintProvider*>()
+	);
+
+	WaterzContext* context = WaterzContext::createNew();
+	context->regionGraph        = regionGraph;
+	context->regionMerging      = regionMerging;
+	context->scoringFunction    = scoringFunction;
+	context->statisticsProvider = statisticsProvider;
+	// context->segmentation       = segmentation_data;
+	context->constraints 		= dummy_constraints;
+
+	if (segmentation_data != NULL) {
+		// wrap data (no copy)
+		volume_ref_ptr<SegID> segmentation(
+				new volume_ref<SegID>(
+						segmentation_data,
+						boost::extents[width][height][depth]
+				)
+		);
+		context->segmentation = segmentation;
+	}
+
+	WaterzState initial_state;
+	initial_state.context = context->id;
+
+	return initial_state;
+}
+
 WaterzState
 initialize(
 		std::size_t     width,
@@ -21,8 +102,18 @@ initialize(
 		const AffValue* affinity_data,
 		SegID*          segmentation_data,
 		const GtID*     ground_truth_data,
+		const SemValue* semantic_data,
+		const SegID*    segconstraint_data,
 		AffValue        affThresholdLow,
 		AffValue        affThresholdHigh,
+		AffValue        semantic_aff_threshold,
+		size_t          semantic_size_threshold,
+		AffValue        semantic_signal_ratio,
+		const std::vector<SemValue>& semantic_taint_labels,
+		AffValue        semantic_taint_threshold,
+		AffValue        size_heuristic_aff_threshold,
+		size_t          size_heuristic_small_threshold,
+		size_t          size_heuristic_large_threshold,
 		bool            findFragments) {
 
 	std::size_t num_voxels = width*height*depth;
@@ -44,15 +135,10 @@ initialize(
 	counts_t<std::size_t> sizes;
 
 	if (findFragments) {
-
 		std::cout << "performing initial watershed segmentation..." << std::endl;
-
 		watershed(affinities, affThresholdLow, affThresholdHigh, *segmentation, sizes);
-
 	} else {
-
 		std::cout << "counting regions and sizes..." << std::endl;
-
 		std::size_t maxId = *std::max_element(segmentation_data, segmentation_data + num_voxels);
 		sizes.resize(maxId + 1);
 		for (std::size_t i = 0; i < num_voxels; i++)
@@ -88,12 +174,58 @@ initialize(
 			new RegionMergingType(*regionGraph)
 	);
 
+	std::shared_ptr<vector<ConstraintProvider*>> constraints(
+		new vector<ConstraintProvider*>()
+	);
+
+	if (semantic_data != NULL) {
+		std::cout << "getting semantic information..." << std::endl;
+		constraints->push_back(
+			new SemanticConstraintProvider<RegionGraphType, SemValue, SegID>(
+				semantic_data, segmentation_data, num_voxels,
+				semantic_aff_threshold,
+				semantic_size_threshold,
+				semantic_signal_ratio
+			)
+		);
+	}
+
+	if (semantic_data != NULL && !semantic_taint_labels.empty()) {
+		std::cout << "getting semantic taint constraint information..." << std::endl;
+		constraints->push_back(
+			new SemanticTaintConstraintProvider<RegionGraphType, SemValue, SegID>(
+				semantic_data, segmentation_data, num_voxels,
+				semantic_taint_labels,
+				semantic_taint_threshold
+			)
+		);
+	}
+
+	if (segconstraint_data != NULL) {
+		std::cout << "getting seg constraint information..." << std::endl;
+		constraints->push_back(
+			new SegConstraintProvider<RegionGraphType, SegID>(
+				segconstraint_data, segmentation_data, num_voxels
+			)
+		);
+	}
+
+	constraints->push_back(
+		new SizeHeuristicConstraintProvider<RegionGraphType, SegID>(
+			sizes,
+			size_heuristic_aff_threshold,
+			size_heuristic_small_threshold,
+			size_heuristic_large_threshold
+		)
+	);
+
 	WaterzContext* context = WaterzContext::createNew();
 	context->regionGraph        = regionGraph;
 	context->regionMerging      = regionMerging;
 	context->scoringFunction    = scoringFunction;
 	context->statisticsProvider = statisticsProvider;
 	context->segmentation       = segmentation;
+	context->constraints 		= constraints;
 
 	WaterzState initial_state;
 	initial_state.context = context->id;
@@ -129,10 +261,12 @@ mergeUntil(
 	std::size_t merged = context->regionMerging->mergeUntil(
 			*context->scoringFunction,
 			*context->statisticsProvider,
+			*context->constraints,
 			threshold,
-			mergeHistoryVisitor);
+			mergeHistoryVisitor
+		);
 
-	if (merged) {
+	if (merged && context->segmentation) {
 
 		std::cout << "extracting segmentation" << std::endl;
 
@@ -164,6 +298,16 @@ getRegionGraph(WaterzState& state) {
 	return regionMerging->extractRegionGraph<ScoredEdge>(*scoringFunction);
 }
 
+std::vector<double>
+getRegionGraphMeta(WaterzState& state) {
+
+	WaterzContext* context = WaterzContext::get(state.context);
+	std::shared_ptr<RegionMergingType> regionMerging = context->regionMerging;
+	std::shared_ptr<ScoringFunctionType> scoringFunction = context->scoringFunction;
+
+	return regionMerging->extractRegionGraphMeta<double>(*scoringFunction, *(context->statisticsProvider));
+}
+
 void
 free(WaterzState& state) {
 
diff --git a/src/waterz/frontend_agglomerate.h b/src/waterz/frontend_agglomerate.h
index 7724ed7..9551a30 100644
--- a/src/waterz/frontend_agglomerate.h
+++ b/src/waterz/frontend_agglomerate.h
@@ -11,10 +11,13 @@
 #include "backend/PriorityQueue.hpp"
 #include "backend/HistogramQuantileProvider.hpp"
 #include "backend/VectorQuantileProvider.hpp"
+#include "backend/ConstraintProvider.hpp"
 #include "evaluate.hpp"
 
+
 typedef uint64_t SegID;
 typedef uint32_t GtID;
+typedef uint8_t SemValue;
 typedef float AffValue;
 typedef float ScoreValue;
 typedef RegionGraph<SegID> RegionGraphType;
@@ -44,6 +47,8 @@ struct Merge {
 
 struct ScoredEdge {
 
+	ScoredEdge() {};
+
 	ScoredEdge(SegID u_, SegID v_, ScoreValue score_) :
 		u(u_),
 		v(v_),
@@ -99,6 +104,7 @@ class WaterzContext {
 	std::shared_ptr<RegionMergingType> regionMerging;
 	std::shared_ptr<ScoringFunctionType> scoringFunction;
 	std::shared_ptr<StatisticsProviderType> statisticsProvider;
+	std::shared_ptr<vector<ConstraintProvider*>> constraints;
 	volume_ref_ptr<SegID> segmentation;
 	volume_const_ref_ptr<GtID> groundtruth;
 
@@ -147,16 +153,36 @@ WaterzState initialize(
 		size_t          depth,
 		const AffValue* affinity_data,
 		SegID*          segmentation_data,
-		const GtID*     groundtruth_data = NULL,
-		AffValue        affThresholdLow  = 0.0001,
-		AffValue        affThresholdHigh = 0.9999,
-		bool            findFragments = true);
+		const GtID*     groundtruth_data,
+		const SemValue* semantic_data,
+		const SegID*    segconstraint_data,
+		AffValue        affThresholdLow,
+		AffValue        affThresholdHigh,
+		AffValue        semantic_aff_threshold,
+		size_t          semantic_size_threshold,
+		AffValue        semantic_signal_ratio,
+		const std::vector<SemValue>& semantic_taint_labels,
+		AffValue        semantic_taint_threshold,
+		AffValue        size_heuristic_aff_threshold,
+		size_t          size_heuristic_small_threshold,
+		size_t          size_heuristic_large_threshold,
+		bool            findFragments
+	);
+
+WaterzState initialize_with_rag(
+		const std::vector<ScoredEdge>& rag,
+		const std::vector<double>& rag_metadata,
+		SegID*          segmentation_data,
+		std::size_t     width,
+		std::size_t     height,
+		std::size_t     depth);
 
 std::vector<Merge> mergeUntil(
 		WaterzState& state,
 		float        threshold);
 
 std::vector<ScoredEdge> getRegionGraph(WaterzState& state);
+std::vector<double> getRegionGraphMeta(WaterzState& state);
 
 void free(WaterzState& state);