feat: add return_contact feature to cross_section

src/fastxs3d.cpp

@@ -30,7 +30,7 @@ auto section(
     float* data = static_cast<float*>(arr.request().ptr);
     std::fill(data, data + voxels, 0.0f);
 
-	xs3d::cross_section(
+	std::tuple<float*, uint8_t> tup = xs3d::cross_section(
 		binimg.data(),
 		sx, sy, sz,
 		point.at(0), point.at(1), point.at(2),
@@ -39,7 +39,7 @@ auto section(
 		data
 	);
 
-	return arr;
+	return std::make_tuple(arr, std::get<1>(tup));
 }
 
 auto area(

src/xs3d.hpp

@@ -665,7 +665,7 @@ float cross_sectional_area(
 	);
 }
 
-float* cross_section(
+std::tuple<float*, uint8_t> cross_section(
 	const uint8_t* binimg,
 	const uint64_t sx, const uint64_t sy, const uint64_t sz,
 
@@ -678,30 +678,31 @@ float* cross_section(
 		plane_visualization = new float[sx * sy * sz]();
 	}
 
+	uint8_t contact = 0;
+
 	if (px < 0 || px >= sx) {
-		return plane_visualization;
+		return std::make_tuple(plane_visualization, contact);
 	}
 	else if (py < 0 || py >= sy) {
-		return plane_visualization;
+		return std::make_tuple(plane_visualization, contact);
 	}
 	else if (pz < 0 || pz >= sz) {
-		return plane_visualization;
+		return std::make_tuple(plane_visualization, contact);
 	}
 
 	uint64_t loc = static_cast<uint64_t>(px) + sx * (
 		static_cast<uint64_t>(py) + sy * static_cast<uint64_t>(pz)
 	);
 
 	if (!binimg[loc]) {
-		return plane_visualization;
+		return std::make_tuple(plane_visualization, contact);
 	}
 
 	const Vec3 pos(px, py, pz);
 	const Vec3 anisotropy(wx, wy, wz);
 	Vec3 normal(nx, ny, nz);
 	normal /= normal.norm();
 
-	uint8_t contact = 0;
 
 	cross_sectional_area_helper(
 		binimg, 
@@ -710,7 +711,7 @@ float* cross_section(
 		contact, plane_visualization
 	);
 
-	return plane_visualization;
+	return std::make_tuple(plane_visualization, contact);
 }
 
 std::tuple<Vec3, Vec3> create_orthonormal_basis(

xs3d/__init__.py

@@ -76,6 +76,7 @@ def cross_section(
   pos:Sequence[int],
   normal:Sequence[float],
   anisotropy:Optional[Sequence[float]] = None,
+  return_contact:bool = False,
 ) -> np.ndarray:
   """
   Compute which voxels are intercepted by a section plane
@@ -91,6 +92,17 @@ def cross_section(
     e.g. [4,4,40] for an electron microscope image with 
     4nm XY resolution with a 40nm cutting plane in 
     serial sectioning.
+  return_contact: if true, return a tuple of (area, contact)
+    where area is the usual output and contact is non-zero if
+    the section plane has contacted the edge of the image
+    indicating the area may be an underestimate if you are
+    working with a cutout of a larger image.
+
+    Contact is an 8-bit bitfield that represents which image faces
+    have been touched. The bits are organized as follows.
+
+    0: 0 X     2: 0 Y     4: 0 Z      6: Unused
+    1: Max X   3: Max Y   5: Max Z    7: Unused
 
   Returns: float32 volume where each voxel's value is its
     contribution to the cross sectional area
@@ -113,10 +125,15 @@ def cross_section(
 
   binimg = np.asfortranarray(binimg)
 
-  if binimg.ndim == 3:
-    return fastxs3d.section(binimg, pos, normal, anisotropy)
-  else:
+  if binimg.ndim != 3:
     raise ValueError("dimensions not supported")
+
+  section, contact = fastxs3d.section(binimg, pos, normal, anisotropy)
+
+  if return_contact:
+    return (section, contact)
+
+  return section
 
 def slice(
   labels:np.ndarray,