perf: persistent buffer with colored pens

src/aux.hpp

@@ -190,14 +190,13 @@ std::tuple<float*, uint8_t> cross_section_slow(
  * connected component, so pre-filtering must be performed to 
  * ensure a match.
  */
-float cross_sectional_area_slow(
+std::tuple<float, uint8_t> cross_sectional_area_slow(
 	const uint8_t* binimg,
 	const uint64_t sx, const uint64_t sy, const uint64_t sz,
 
 	const float px, const float py, const float pz,
 	const float nx, const float ny, const float nz,
-	const float wx, const float wy, const float wz,
-	uint8_t &contact = _dummy_contact
+	const float wx, const float wy, const float wz
 ) {
 
 	const Vec3 pos(px, py, pz);
@@ -208,7 +207,7 @@ float cross_sectional_area_slow(
 		|| rpos.y < 0 || rpos.y >= sy 
 		|| rpos.z < 0 || rpos.z >= sz
 	) {
-		return 0.0;
+		return std::make_tuple(0.0, 0);
 	}
 
 	const Vec3 anisotropy(wx, wy, wz);
@@ -233,7 +232,7 @@ float cross_sectional_area_slow(
 			: 1.0 / projections[i];
 	}
 
-	contact = 0;
+	uint8_t contact = 0;
 
 	for (uint64_t z = 0; z < sz; z++) {
 		for (uint64_t y = 0; y < sy; y++) {
@@ -263,7 +262,7 @@ float cross_sectional_area_slow(
 		}
 	}
 
-	return area;
+	return std::make_tuple(area, contact);
 }
 
 };

src/fastxs3d.cpp

@@ -84,7 +84,8 @@ auto calculate_area(
 	const py::array_t<float> &point,
 	const py::array_t<float> &normal,
 	const py::array_t<float> &anisotropy,
-	const bool slow_method
+	const bool slow_method,
+	const bool use_persistent_data
 ) {
 	const uint64_t sx = binimg.shape()[0];
 	const uint64_t sy = binimg.ndim() < 2
@@ -94,31 +95,25 @@ auto calculate_area(
 		? 1 
 		: binimg.shape()[2];
 
-	uint8_t contact = false;
-	float area = 0;
-
 	if (slow_method) {
-		area = xs3d::cross_sectional_area_slow(
+		return xs3d::cross_sectional_area_slow(
 			binimg.data(),
 			sx, sy, sz,
 			point.at(0), point.at(1), point.at(2),
 			normal.at(0), normal.at(1), normal.at(2),
-			anisotropy.at(0), anisotropy.at(1), anisotropy.at(2),
-			contact
+			anisotropy.at(0), anisotropy.at(1), anisotropy.at(2)
 		);
 	}
 	else {
-		area = xs3d::cross_sectional_area(
+		return xs3d::cross_sectional_area(
 			binimg.data(),
 			sx, sy, sz,
 			point.at(0), point.at(1), point.at(2),
 			normal.at(0), normal.at(1), normal.at(2),
 			anisotropy.at(0), anisotropy.at(1), anisotropy.at(2),
-			contact
+			use_persistent_data
 		);
 	}
-
-	return std::tuple(area, contact);
 }
 
 auto projection(	
@@ -225,4 +220,7 @@ PYBIND11_MODULE(fastxs3d, m) {
 	m.def("projection", &projection, "Project a cross section of a 3D image onto a 2D plane");
 	m.def("section", &section, "Return a floating point image that shows the voxels contributing area to a cross section.");
 	m.def("area", &calculate_area, "Find the cross sectional area for a given binary image, point, and normal vector.");
+
+	m.def("set_shape", &xs3d::set_shape, "Accelerate the area function across many evaluation points by saving some attributes of the input shape upfront. Call clear_shape when you are done.");
+	m.def("clear_shape", &xs3d::clear_shape, "Delete the data that was persisted by set_shape.");
 }