feat: true 2d contact analysis w/ fix for area calculation (#142)

* feat: true 2d contact analysis w/ fix for area calculation

* perf: remove extraneous calculation

* perf: don't calculate for label==label

* fix: update cython code to use 4,8 correctly

* fix: use correct anisotropy

* fix: remove unused variable

* test: check 4 and 8 connected

* test: an additional orientation

* fix: ensure anisotropy is good for 2d or 3d

* test: check if surface area makes sense for 2d

* fix: ensure surface_area=False still works

Author: william-silversmith

automated_test.py

@@ -731,6 +731,15 @@ def test_contacts_surface_area():
   res = cc3d.contacts(labels_out, connectivity=6, surface_area=True)
   assert res[(1,2)] == 6
 
+  data = np.ones((3,3), dtype=np.uint32)
+  data[1,1] = 2
+  res = cc3d.contacts(data, connectivity=8, surface_area=True, anisotropy=(1.5, 2.3))
+  assert np.isclose(res[(1,2)], (1.5 * 2 + 2.3 * 2))
+
+  res = cc3d.contacts(data, connectivity=4, surface_area=True, anisotropy=(1.5, 2.3))
+  assert np.isclose(res[(1,2)], (1.5 * 2 + 2.3 * 2))
+
+
 def test_contacts_26():
   labels = np.zeros( (10, 10, 10), dtype=np.uint32 )
 
@@ -846,14 +855,25 @@ def test_contacts_2d():
   labels[7,2] = 6
   labels[8,2] = 7
 
+  labels[9,9] = 8
+  labels[8,8] = 9
+
+  labels[6,8] = 10
+  labels[5,9] = 11
+
   # not connected to anything else
-  labels[1,:] = 10
+  labels[1,:] = 12
 
-  res = cc3d.contacts(labels, connectivity=8)
+  res = cc3d.contacts(labels, connectivity=4)
   assert set(res.keys()) == set([ 
     (2,3), (4,5), (6,7)
   ])
 
+  res = cc3d.contacts(labels, connectivity=8)
+  assert set(res.keys()) == set([ 
+    (2,3), (4,5), (6,7), (8,9), (10,11),
+  ])
+
 def test_voxel_graph_2d():
   labels = np.ones((3,3), dtype=np.uint8)
   graph = cc3d.voxel_connectivity_graph(labels, connectivity=4)

cc3d/cc3d_graphs.hpp

@@ -279,6 +279,23 @@ inline void compute_neighborhood(
 	neighborhood[12] = (connectivity > 18) * (plus_x + plus_y + minus_z) * (plus_y && minus_z);
 }
 
+inline void compute_neighborhood(
+	int64_t *neighborhood, 
+	const int64_t x, const int64_t y,
+	const int64_t sx, const int64_t sy,
+	const int64_t connectivity = 8
+) {
+
+	const int64_t plus_x = (x < (static_cast<int64_t>(sx) - 1)); // +x
+	const int64_t minus_x = -1 * (x > 0); // -x
+	const int64_t minus_y = -static_cast<int64_t>(sx) * (y > 0); // -y
+
+	neighborhood[0] = minus_x;
+	neighborhood[1] = minus_y;
+	neighborhood[2] = (connectivity > 4) * (minus_x + minus_y);
+	neighborhood[3] = (connectivity > 4) * (plus_x + minus_y);
+}
+
 struct pair_hash {
 	inline std::size_t operator()(const std::pair<uint64_t,uint64_t> & v) const {
 		return v.first * 31 + v.second; // arbitrary hash fn
@@ -287,7 +304,66 @@ struct pair_hash {
 
 template <typename T>
 const std::unordered_map<std::pair<T,T>, float, pair_hash> 
-extract_region_graph(
+extract_region_graph_2d(
+	T* labels, 
+	const int64_t sx, const int64_t sy,
+	const float wx=1, const float wy=1,
+	const int64_t connectivity=8,
+	const bool surface_area=true
+) {
+
+	if (connectivity != 4 && connectivity != 8) {
+		throw std::runtime_error("Only 4 and 8 connectivities are supported.");
+	}
+
+	int64_t neighborhood[4];
+	float areas[4] = {
+		wy, wx, // edges
+		0, 0 // corners
+	};
+	if (!surface_area) {
+		std::fill(areas, areas+4, 1.0);
+	}
+
+	T cur = 0;
+	T label = 0;
+
+	std::unordered_map<std::pair<T,T>, float, pair_hash> edges;
+
+	for (int64_t y = 0; y < sy; y++) {
+		for (int64_t x = 0; x < sx; x++) {
+			int64_t loc = x + sx * y;
+			cur = labels[loc];
+
+			if (cur == 0) {
+				continue;
+			}
+
+			compute_neighborhood(neighborhood, x, y, sx, sy, connectivity);
+
+			for (int i = 0; i < connectivity / 2; i++) {
+				int64_t neighboridx = loc + neighborhood[i];
+				label = labels[neighboridx];
+
+				if (label == 0 || label == cur) {
+					continue;
+				}
+				else if (cur > label) {
+					edges[std::pair<T,T>(label, cur)] += areas[i];
+				}
+				else if (cur < label) {
+					edges[std::pair<T,T>(cur, label)] += areas[i];
+				}
+			}
+		}
+	}
+
+	return edges;
+}
+
+template <typename T>
+const std::unordered_map<std::pair<T,T>, float, pair_hash> 
+extract_region_graph_3d(
 	T* labels, 
 	const int64_t sx, const int64_t sy, const int64_t sz,
 	const float wx=1, const float wy=1, const float wz=1,
@@ -339,7 +415,7 @@ extract_region_graph(
 					int64_t neighboridx = loc + neighborhood[i];
 					label = labels[neighboridx];
 
-					if (label == 0) {
+					if (label == 0 || label == cur) {
 						continue;
 					}
 					else if (cur > label) {
@@ -356,6 +432,41 @@ extract_region_graph(
 	return edges;
 }
 
+
+template <typename T>
+const std::unordered_map<std::pair<T,T>, float, pair_hash> 
+extract_region_graph(
+	T* labels, 
+	const int64_t sx, const int64_t sy, const int64_t sz,
+	const float wx=1, const float wy=1, const float wz=1,
+	const int64_t connectivity=26,
+	const bool surface_area=true
+) {
+	if (connectivity == 4 || connectivity == 8) {
+		if (sz != 1) {
+			throw std::runtime_error("z thickness must be 1 for 2d region graph extraction.");
+		}
+
+		return extract_region_graph_2d<T>(
+			labels, 
+			sx, sy,
+			wx, wy,
+			connectivity, surface_area
+		);
+	}
+	else if (connectivity == 6 || connectivity == 18 || connectivity == 26) {
+		return extract_region_graph_3d<T>(
+			labels, 
+			sx, sy, sz,
+			wx, wy, wz,
+			connectivity, surface_area
+		);
+	}
+	else {
+		throw std::runtime_error("Only (2d) 4, 8, or (3d) 6, 18, and 26 connectivities are supported.");
+	}
+}
+
 template <typename T>
 std::map<T, std::vector<std::pair<size_t, size_t>>> 
 extract_runs(T* labels, const size_t voxels) {

cc3d/fastcc3d.pyx

@@ -1185,7 +1185,7 @@ def contacts(
   labels:np.ndarray, 
   int connectivity=26, 
   surface_area:bool = True, 
-  anisotropy:Tuple[int,int,int] = (1,1,1)
+  anisotropy:Union[Tuple[float,float,float], Tuple[float,float]] = (1,1,1)
 ) -> Dict[Tuple[int,int], float]:
   """
   Get the N-connected region adjacancy graph of a 3D image
@@ -1207,6 +1207,9 @@ def contacts(
   while len(labels.shape) < 3:
     labels = labels[..., np.newaxis ]
 
+  while len(anisotropy) < 3:
+    anisotropy = tuple(anisotropy) + (1,)
+
   return _contacts(labels, connectivity, surface_area, anisotropy)
 
 def _contacts(
@@ -1215,13 +1218,8 @@ def _contacts(
   surface_area=True,
   anisotropy=(1,1,1), 
 ):
-  if connectivity == 8 and labels.shape[2] == 1:
-    connectivity = 26
-  if connectivity == 4 and labels.shape[2] == 1:
-    connectivity = 6
-
-  if connectivity not in (6, 18, 26):
-    raise ValueError("Only 6, 18, and 26 connectivities are supported. Got: " + str(connectivity))
+  if connectivity not in (4, 8, 6, 18, 26):
+    raise ValueError(f"Only (2d) 4, 8, (3d) 6, 18, and 26 connectivities are supported. Got: {connectivity}")
 
   labels = np.asfortranarray(labels)