feat: slice operator (#6)

* feat: projection operator

* fix: ensure other types are handled right

* fix: projections need to handle sqrt(3) in any direction

* docs: better help string

* fix: ensure no off-by-one error and use vector<bool>

* wip: trying to reduce the area of the projection

* fix: ensure that for normal=0,0,1, an xy plane is correcly rendered

Previously, at least one axis was reflected.

* fix: warning on Bbox2d print statement

* fix: cutout and fortran orientation of array

* fix: change basis calculation to produce XY plane in std orientation

Previously was reflecting the X axis when given a positive z normal.

* redesign: change order of stack for easier debugging

* fix: only modify bbox if valid coordinate

* redesign: change name of projection to slice

* fix: use rational appx of sqrt(3) for reliable behavior

floating point multiplication could be interpreted differently
on various compilers, architectures

* fix: use a more accurate rational approximation

* docs: show how to use new slice operator

* fix: replace sqrt(3) with rational appx in another location

* test: add the most basic test of the slice operator

* perf: remove diagnonal moves

Not necessary to explore in an 8-connected fashion since 4 will
fully flood fill the plane since there are no obstacles.

README.md

@@ -39,15 +39,18 @@ area, contact_warning = xs3d.cross_sectional_area(
 	return_contact=True
 )
 
-# Returns the cross section as a float32 image
+# Returns the cross section as a float32 3d image
 # where each voxel represents its contribution
 # to the cross sectional area
 image = xs3d.cross_section(
 	binary_image, vertex, 
 	normal, resolution, 
 )
 
-
+# Get a slice of a 3d image in any orientation.
+# Note: result may be reflected or transposed
+# compared with what you might expect.
+image2d = xs3d.slice(labels, vertex, normal)
 ```
 
 When using skeletons (one dimensional stick figure representations) to create electrophysiological compartment simulations of neurons, some additional information is required for accuracy. The caliber of the neurite changes over the length of the cell.

automated_test.py

@@ -304,7 +304,10 @@ def angle(theta):
 		assert image.dtype == np.float32
 		assert np.isclose(image.sum(), area)
 
-
+def test_slice():
+	labels = np.arange(9, dtype=np.uint8).reshape([3,3,1], order="F")
+	slc = xs3d.slice(labels, [0,0,0], [0,0,1])
+	assert np.all(slc == labels[:,:,0])
 
 
 

src/fastxs3d.cpp

@@ -68,8 +68,82 @@ auto area(
 	return std::tuple(area, contact);
 }
 
+auto projection(	
+	const py::array &labels,
+	const py::array_t<float> &point,
+	const py::array_t<float> &normal
+) {
+	const uint64_t sx = labels.shape()[0];
+	const uint64_t sy = labels.ndim() < 2
+		? 1 
+		: labels.shape()[1];
+	const uint64_t sz = labels.ndim() < 3 
+		? 1 
+		: labels.shape()[2];
+
+	// rational approximation of sqrt(3) is 97/56
+	// result is more likely to be same across compilers
+	uint64_t psx = 2 * 97 * std::max(std::max(sx,sy), sz) / 56 + 1;
+	uint64_t pvoxels = psx * psx;
+
+	py::array arr; 
+
+	auto projectionfn = [&](auto dtype) {
+		arr = py::array_t<decltype(dtype), py::array::f_style>({ psx, psx });
+		auto out = reinterpret_cast<decltype(dtype)*>(arr.request().ptr);
+		auto data = reinterpret_cast<decltype(dtype)*>(labels.request().ptr);
+		std::fill(out, out + pvoxels, 0);
+
+		std::tuple<decltype(dtype)*, xs3d::Bbox2d> tup = xs3d::cross_section_projection<decltype(dtype)>(
+			data,
+			sx, sy, sz,
+			point.at(0), point.at(1), point.at(2),
+			normal.at(0), normal.at(1), normal.at(2),
+			out
+		);
+
+		xs3d::Bbox2d bbox = std::get<1>(tup);
+		bbox.x_max++;
+		bbox.y_max++;
+
+		auto cutout = py::array_t<decltype(dtype), py::array::f_style>({ bbox.sx(), bbox.sy() });
+	    auto cutout_ptr = reinterpret_cast<decltype(dtype)*>(cutout.request().ptr);
+
+	    ssize_t csx = bbox.sx();
+
+	    for (ssize_t y = bbox.y_min; y < bbox.y_max; y++) {
+	        for (ssize_t x = bbox.x_min; x < bbox.x_max; x++) {
+	            cutout_ptr[
+	            	(x - bbox.x_min) + csx * (y - bbox.y_min)
+	            ] = out[x + psx * y];
+	        }
+	    }
+
+		return cutout.view(py::str(labels.dtype()));
+	};
+
+	int data_width = labels.dtype().itemsize();
+
+    if (data_width == 1) {
+    	return projectionfn(uint8_t{});
+    }
+    else if (data_width == 2) {
+    	return projectionfn(uint16_t{});
+    }
+    else if (data_width == 4) {
+    	return projectionfn(uint32_t{});
+    }
+    else if (data_width == 8) {
+    	return projectionfn(uint64_t{});
+    }
+    else {
+    	throw new std::runtime_error("should never happen");
+    }
+}
+
 PYBIND11_MODULE(fastxs3d, m) {
 	m.doc() = "Finding cross sectional area in 3D voxelized images."; 
+	m.def("projection", &projection, "Project a cross section of a 3D image onto a 2D plane");
 	m.def("section", &section, "Return a binary image that highlights the voxels contributing area to a cross section.");
 	m.def("area", &area, "Find the cross sectional area for a given binary image, point, and normal vector.");
 }

src/xs3d.hpp

@@ -423,7 +423,9 @@ float cross_sectional_area_helper(
 ) {
 	std::vector<bool> ccl(sx * sy * sz);
 
-	uint64_t plane_size = 2 * sqrt(3) * std::max(std::max(sx,sy), sz) + 1;
+	// rational approximation of sqrt(3) is 97/56
+	// more reliable behavior across compilers/architectures
+	uint64_t plane_size = 2 * 97 * std::max(std::max(sx,sy), sz) / 56 + 1;
 
 	// maximum possible size of plane
 	uint64_t psx = plane_size;
@@ -544,6 +546,27 @@ float cross_sectional_area_helper(
 
 namespace xs3d {
 
+struct Bbox2d {
+	int64_t x_min, x_max;
+	int64_t y_min, y_max;
+	Bbox2d() : x_min(0), x_max(0), y_min(0), y_max(0) {};
+	Bbox2d(int64_t x_min, int64_t x_max, int64_t y_min, int64_t y_max) 
+		: x_min(x_min), x_max(x_max), y_min(y_min), y_max(y_max) {};
+
+	int64_t sx() const {
+		return x_max - x_min;
+	}
+	int64_t sy() const {
+		return y_max - y_min;
+	}
+	int64_t pixels() const {
+		return sx() * sy();
+	}
+	void print() const {
+		printf("Bbox2d(%llu, %llu, %llu, %llu)\n", x_min, x_max, y_min, y_max);
+	}
+};
+
 float cross_sectional_area(
 	const uint8_t* binimg,
 	const uint64_t sx, const uint64_t sy, const uint64_t sz,
@@ -633,6 +656,121 @@ float* cross_section(
 	return plane_visualization;
 }
 
+template <typename LABEL>
+std::tuple<LABEL*, Bbox2d> cross_section_projection(
+	const LABEL* labels,
+	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,
+	LABEL* out = NULL
+) {
+	// maximum possible size of plane
+	// rational approximation of sqrt(3) is 97/56
+	const uint64_t psx = 2 * 97 * std::max(std::max(sx,sy), sz) / 56 + 1;
+	const uint64_t psy = psx;
+
+	Bbox2d bbx;
+
+	std::vector<bool> visited(psx * psy);
+
+	if (out == NULL) {
+		out = new LABEL[psx * psy]();
+	}
+
+	if (px < 0 || px >= sx) {
+		return std::tuple(out, bbx);
+	}
+	else if (py < 0 || py >= sy) {
+		return std::tuple(out, bbx);
+	}
+	else if (pz < 0 || pz >= sz) {
+		return std::tuple(out, bbx);
+	}
+
+	const Vec3 pos(px, py, pz);
+	Vec3 normal(nx, ny, nz);
+	normal /= normal.norm();
+
+	Vec3 basis1 = jhat.cross(normal);
+	if (basis1.is_null()) {
+		basis1 = normal.cross(ihat);
+	}
+	basis1 /= basis1.norm();
+
+	Vec3 basis2 = normal.cross(basis1);
+	basis2 /= basis2.norm();
+
+	uint64_t plane_pos_x = psx / 2;
+	uint64_t plane_pos_y = psy / 2;
+
+	bbx.x_min = plane_pos_x;
+	bbx.x_max = plane_pos_x;
+	bbx.y_min = plane_pos_y;
+	bbx.y_max = plane_pos_y;
+
+	uint64_t ploc = plane_pos_x + psx * plane_pos_y;
+
+	std::stack<uint64_t> stack;
+	stack.push(ploc);
+
+	while (!stack.empty()) {
+		ploc = stack.top();
+		stack.pop();
+
+		if (visited[ploc]) {
+			continue;
+		}
+
+		visited[ploc] = true;
+
+		uint64_t y = ploc / psx;
+		uint64_t x = ploc - y * psx;
+
+		float dx = static_cast<float>(x) - static_cast<float>(plane_pos_x);
+		float dy = static_cast<float>(y) - static_cast<float>(plane_pos_y);
+
+		Vec3 cur = pos + basis1 * dx + basis2 * dy;
+
+		if (cur.x < 0 || cur.y < 0 || cur.z < 0) {
+			continue;
+		}
+		else if (cur.x >= sx || cur.y >= sy || cur.z >= sz) {
+			continue;
+		}
+
+		bbx.x_min = std::min(bbx.x_min, static_cast<int64_t>(x));
+		bbx.x_max = std::max(bbx.x_max, static_cast<int64_t>(x));
+		bbx.y_min = std::min(bbx.y_min, static_cast<int64_t>(y));
+		bbx.y_max = std::max(bbx.y_max, static_cast<int64_t>(y));
+
+		uint64_t loc = static_cast<uint64_t>(cur.x) + sx * (
+			static_cast<uint64_t>(cur.y) + sy * static_cast<uint64_t>(cur.z)
+		);
+
+		out[ploc] = labels[loc];
+
+		uint64_t up = ploc - psx; 
+		uint64_t down = ploc + psx;
+		uint64_t left = ploc - 1;
+		uint64_t right = ploc + 1;
+
+		if (x > 0 && !visited[left]) {
+			stack.push(left);
+		}
+		if (x < psx - 1 && !visited[right]) {
+			stack.push(right);
+		}
+		if (y > 0 && !visited[up]) {
+			stack.push(up);
+		}
+		if (y < psy - 1 && !visited[down]) {
+			stack.push(down);
+		}
+	}
+
+	return std::tuple(out, bbx);
+}
 
 };
 

xs3d/__init__.py

@@ -118,3 +118,46 @@ def cross_section(
   else:
     raise ValueError("dimensions not supported")
 
+def slice(
+  labels:np.ndarray,
+  pos:Sequence[int],
+  normal:Sequence[float],
+) -> np.ndarray:
+  """
+  Compute which voxels are intercepted by a section plane
+  and project them onto a plane.
+
+  NB: The orientation of this projection is not guaranteed. 
+  The axes can be reflected and transposed compared to what
+  you might expect.
+
+  labels: a binary 2d or 3d numpy image (e.g. a bool datatype)
+  pos: the point in the image from which to extract the section
+    must be an integer (it's an index into the image).
+    e.g. [5,10,2]
+  normal: a vector normal to the section plane, does not
+    need to be a unit vector. e.g. [sqrt(2)/2. sqrt(2)/2, 0]
+
+  Returns: ndarray
+  """
+  pos = np.array(pos, dtype=np.float32)
+  normal = np.array(normal, dtype=np.float32)
+
+  if np.all(normal == 0):
+    raise ValueError("normal vector must not be a null vector (all zeros).")
+
+  labels = np.asfortranarray(labels)
+
+  if labels.ndim != 3:
+    raise ValueError(f"{labels.ndim} dimensions not supported")
+
+  return fastxs3d.projection(labels, pos, normal)
+
+
+
+
+
+
+
+
+