diff --git a/gala/features/base.py b/gala/features/base.py
index 629c17c..1edd9c3 100644
--- a/gala/features/base.py
+++ b/gala/features/base.py
@@ -33,10 +33,6 @@ def update_node_cache(self, *args, **kwargs):
         pass
     def update_edge_cache(self, *args, **kwargs):
         pass
-    def pixelwise_update_node_cache(self, *args, **kwargs):
-        pass
-    def pixelwise_update_edge_cache(self, *args, **kwargs):
-        pass
     def compute_node_features(self, *args, **kwargs):
         return np.array([])
     def compute_edge_features(self, *args, **kwargs):
@@ -69,14 +65,6 @@ def update_edge_cache(self, g, e1, e2, dst, src):
         for i, child in enumerate(self.children):
             child.update_edge_cache(g, e1, e2, dst[i], src[i])
     
-    def pixelwise_update_node_cache(self, g, n, dst, idxs, remove=False):
-        for i, child in enumerate(self.children):
-            child.pixelwise_update_node_cache(g, n, dst[i], idxs, remove)
-
-    def pixelwise_update_edge_cache(self, g, n1, n2, dst, idxs, remove=False):
-        for i, child in enumerate(self.children):
-            child.pixelwise_update_edge_cache(g, n1, n2, dst[i], idxs, remove)
-
     def compute_node_features(self, g, n, cache=None):
         if cache is None: cache = g.node[n][self.default_cache]
         features = []
diff --git a/gala/features/contact.pyx b/gala/features/contact.pyx
index 0f39604..a693bb5 100644
--- a/gala/features/contact.pyx
+++ b/gala/features/contact.pyx
@@ -70,15 +70,6 @@ class Manager(base.Null):
     def update_edge_cache(self, g, e1, e2, dst, src):
         dst += src
 
-    def pixelwise_update_edge_cache(self, g, n1, n2, dst, idxs, remove=False):
-        if len(idxs) == 0: return
-        n1_idxs = np.array(list(g.extent(n1)))
-        n2_idxs = np.array(list(g.extent(n2)))
-        a = -1.0 if remove else 1.0
-        if self.oriented: ar = g.oriented_probabilities_r
-        else: ar = g.non_oriented_probabilities_r
-        dst += a * _compute_edge_cache(idxs, n1_idxs, n2_idxs, ar, self.thresholds)
-
 
 cdef _compute_contact_matrix(double[:,:,:] totals, double volume_ratio_1,
                              double volume_ratio_2):
diff --git a/gala/features/convex_hull.py b/gala/features/convex_hull.py
index 430f312..b085c87 100644
--- a/gala/features/convex_hull.py
+++ b/gala/features/convex_hull.py
@@ -84,12 +84,6 @@ def update_edge_cache(self, g, e1, e2, dst, src):
         vol, tri = self.convex_hull_vol(allind, g)
         dst = np.array([tri,vol])
 
-    def pixelwise_update_node_cache(self, g, n, dst, idxs, remove=False):
-        pass
-
-    def pixelwise_update_edge_cache(self, g, n1, n2, dst, idxs, remove=False):
-        pass
-
     def compute_node_features(self, g, n, cache=None):
         if cache is None: 
             cache = g.node[n][self.default_cache]
diff --git a/gala/features/histogram.pyx b/gala/features/histogram.pyx
index 6ea3454..757e7bf 100644
--- a/gala/features/histogram.pyx
+++ b/gala/features/histogram.pyx
@@ -124,29 +124,10 @@ class Manager(base.Null):
     def update_edge_cache(self, g, e1, e2, dst, src):
         dst += src
 
-    def pixelwise_update_node_cache(self, g, n, dst, idxs, remove=False):
-        if len(idxs) == 0: return
-        a = -1.0 if remove else 1.0
-        if self.oriented:
-            ar = g.max_probabilities_r
-        else:
-            ar = g.non_oriented_probabilities_r
-
-        dst += a * self.histogram(ar[idxs,:])
-
-    def pixelwise_update_edge_cache(self, g, n1, n2, dst, idxs, remove=False):
-        if len(idxs) == 0: return
-        a = -1.0 if remove else 1.0
-        if self.oriented:
-            ar = g.oriented_probabilities_r
-        else:
-            ar = g.non_oriented_probabilities_r
-
-        dst += a * self.histogram(ar[idxs,:])
-
     def JS_divergence(self, p, q):
         m = (p+q)/2
         return (self.KL_divergence(p, m) + self.KL_divergence(q, m))/2
+
     def KL_divergence(self, p, q):
         """Return the Kullback-Leibler Divergence between two histograms."""
         kl = []
diff --git a/gala/features/moments.pyx b/gala/features/moments.pyx
index 5f8cbab..ca61d69 100644
--- a/gala/features/moments.pyx
+++ b/gala/features/moments.pyx
@@ -57,24 +57,6 @@ class Manager(base.Null):
     def update_edge_cache(self, g, e1, e2, dst, src):
         dst += src
 
-    def pixelwise_update_node_cache(self, g, n, dst, idxs, remove=False):
-        if len(idxs) == 0: return
-        a = -1.0 if remove else 1.0
-        if self.oriented:
-            ar = g.max_probabilities_r
-        else:
-            ar = g.non_oriented_probabilities_r
-        dst += a * self.compute_moment_sums(ar, idxs)
-
-    def pixelwise_update_edge_cache(self, g, n1, n2, dst, idxs, remove=False):
-        if len(idxs) == 0: return
-        a = -1.0 if remove else 1.0
-        if self.oriented:
-            ar = g.max_probabilities_r
-        else:
-            ar = g.non_oriented_probabilities_r
-        dst += a * self.compute_moment_sums(ar, idxs)
-
     def compute_node_features(self, g, n, cache=None):
         if cache is None: 
             cache = g.node[n][self.default_cache]
diff --git a/gala/features/squiggliness.py b/gala/features/squiggliness.py
index 83d1d0b..c20e771 100644
--- a/gala/features/squiggliness.py
+++ b/gala/features/squiggliness.py
@@ -36,19 +36,3 @@ def update_edge_cache(self, g, e1, e2, dst, src):
             (dst[np.newaxis,self.ndim:], src[np.newaxis,self.ndim:]),
             axis=0).max(axis=0)
 
-    def pixelwise_update_edge_cache(self, g, n1, n2, dst, idxs, remove=False):
-        if remove:
-            pass
-            # dst = self.create_edge_cache(g, n1, n2)
-        if len(idxs) == 0: return
-        b = np.concatenate(
-            self.compute_bounding_box(idxs, g.watershed.shape))
-        self.update_edge_cache(g, (n1,n2), None, dst, b)
-
-    def compute_edge_features(self, g, n1, n2, cache=None):
-        if cache is None: 
-            cache = g[n1][n2][self.default_cache]
-        m, M = cache[:self.ndim], cache[self.ndim:]
-        plane_surface = np.sort(M-m)[1:].prod() * (3.0-g.pad_thickness)
-        return np.array([len(g[n1][n2]['boundary']) / plane_surface])
-