wip: experimenting with a smooth path feature

src/fastxs3d.cpp

@@ -155,8 +155,98 @@ auto projection(
     }
 }
 
+auto projection_with_frame(	
+	const py::array &labels,
+	const py::array_t<float> &point,
+	const py::array_t<float> &basis1,
+	const py::array_t<float> &basis2,
+	const py::array_t<float> &anisotropy
+) {
+	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];
+
+	float wx = anisotropy.at(0);
+	float wy = anisotropy.at(1);
+	float wz = anisotropy.at(2);
+
+	float minval = std::min(std::min(wx,wy), wz);
+	wx /= minval; wy /= minval; wz /= minval;
+	float maxval = std::max(std::max(std::abs(wx), std::abs(wy)), std::abs(wz));
+
+	const uint64_t distortion = static_cast<uint64_t>(ceil(maxval));
+
+	// rational approximation of sqrt(3) is 97/56
+	// result is more likely to be same across compilers
+	uint64_t psx = distortion * 2 * 97 * std::max(std::max(sx,sy), sz) / 56 + 1;
+	uint64_t pvoxels = psx * psx;
+
+	py::array arr; 
+
+	xs3d::Vec3 b1(basis1.at(0), basis1.at(1), basis1.at(2));
+	xs3d::Vec3 b2(basis2.at(0), basis2.at(1), basis2.at(2));
+
+	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),
+			b1, b2,
+			anisotropy.at(0), anisotropy.at(1), anisotropy.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);
+
+	    int64_t csx = bbox.sx();
+
+	    for (int64_t y = bbox.y_min; y < bbox.y_max; y++) {
+	        for (int64_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_with_frame", &projection_with_frame, "Project a cross section of a 3D image onto a 2D plane");
 	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

@@ -10,9 +10,7 @@
 #include <vector>
 #include <stdexcept>
 
-namespace {
-
-static uint8_t _dummy_contact = false;
+namespace xs3d {
 
 class Vec3 {
 public:
@@ -137,6 +135,14 @@ class Vec3 {
 	}
 };
 
+};
+
+namespace {
+
+using namespace xs3d;
+
+static uint8_t _dummy_contact = false;
+
 const Vec3 ihat = Vec3(1,0,0);
 const Vec3 jhat = Vec3(0,1,0);
 const Vec3 khat = Vec3(0,0,1);
@@ -763,11 +769,9 @@ 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,
+	const Vec3& basis_1, const Vec3& basis_2,
 	const float wx, const float wy, const float wz,
-	const bool positive_basis,
 	LABEL* out = NULL
 ) {
 
@@ -802,12 +806,9 @@ std::tuple<LABEL*, Bbox2d> cross_section_projection(
 	}
 
 	const Vec3 pos(px, py, pz);
-	Vec3 normal(nx, ny, nz);
-	normal /= normal.norm();
 
-	auto bases = create_orthonormal_basis(normal, positive_basis);
-	Vec3 basis1 = std::get<0>(bases) * anisotropy;
-	Vec3 basis2 = std::get<1>(bases) * anisotropy;
+	Vec3 basis1 = basis_1 / basis_1.norm() * anisotropy;
+	Vec3 basis2 = basis_2 / basis_2.norm() * anisotropy;
 
 	uint64_t plane_pos_x = psx / 2;
 	uint64_t plane_pos_y = psy / 2;
@@ -880,6 +881,37 @@ std::tuple<LABEL*, Bbox2d> cross_section_projection(
 	return std::tuple(out, bbx);
 }
 
+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,
+	const float wx, const float wy, const float wz,
+	const bool positive_basis,
+	LABEL* out = NULL
+) {
+	Vec3 normal(nx, ny, nz);
+	normal /= normal.norm();
+
+	Vec3 anisotropy(wx, wy, wz);
+	anisotropy /= anisotropy.min();
+	anisotropy = Vec3(1,1,1) / anisotropy;
+
+	auto bases = create_orthonormal_basis(normal, positive_basis);
+	Vec3 basis1 = std::get<0>(bases) * anisotropy;
+	Vec3 basis2 = std::get<1>(bases) * anisotropy;
+
+	return cross_section_projection<LABEL>(
+		labels,
+		sx, sy, sz,
+		px, py, pz,
+		basis1, basis2,
+		wx, wy, wz,
+		out
+	);
+}
+
 };
 
 #endif

xs3d/__init__.py

@@ -1,4 +1,4 @@
-from typing import Sequence, Optional
+from typing import Sequence, Optional, Tuple, Generator
 
 from .twod import cross_sectional_area_2d
 import fastxs3d
@@ -190,11 +190,138 @@ def slice(
 
   return fastxs3d.projection(labels, pos, normal, anisotropy, standardize_basis)
 
+def slice_path(
+  labels:np.ndarray,
+  path:Sequence[Sequence[int]],
+  anisotropy:Optional[Sequence[float]] = None,
+  smoothing:int = 1,
+) -> Generator[np.ndarray, None, None]:
+  """
+  Compute which voxels are intercepted by a section plane
+  and project them onto a plane.
+
+  Why is it a generator? Because paths can be artibrarily long
+  and this will avoid running out of memory (e.g. imagine a path
+  that passes through every voxel in the image, for a 512x512x512
+  uint64 volume, that could produce up to 550GB of image data).
+
+  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)
+  path: a sequence of points 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]
+  anisotropy: resolution of the x, y, and z axis
+    e.g. [4,4,40] for an electron microscope image with 
+    4nm XY resolution with a 40nm cutting plane in 
+    serial sectioning.
+  smoothing: number of verticies in the path to smooth the tangent
+    vectors with.
+  """
+  if anisotropy is None:
+    anisotropy = [ 1.0 ] * labels.ndim
 
+  path = np.array(path, dtype=np.float32)
 
+  if path.ndim != 2:
+    raise ValueError("pos must be a sequence of x,y,z points.")
 
+  if labels.ndim != 3:
+    raise ValueError(f"{labels.ndim} dimensions not supported")
 
+  anisotropy = np.array(anisotropy, dtype=np.float32)
+  labels = np.asfortranarray(labels)
 
-
+  if len(path) == 1:
+    yield slice(
+      labels,
+      pos=path[0],
+      normal=[0,0,1],
+      anisotropy=anisotropy,
+    )
+    return
+
+  # vectors aligned with the path
+  tangents = (path[1:] - path[:-1]).astype(np.float32)
+  tangents = np.concatenate([ tangents, [tangents[-1]] ])
+
+  # Running the filter in the forward and then backwards
+  # direction eliminates phase shift.
+  tangents = _moving_average(tangents, smoothing)
+  tangents = _moving_average(tangents[::-1], smoothing)[::-1]
+
+  basis1s = np.zeros([ len(tangents), 3 ], dtype=np.float32)
+  basis2s = np.zeros([ len(tangents), 3 ], dtype=np.float32)
+
+  if len(path) == 2:
+    basis1s[0] = np.cross(tangents[0], [0,1,0])
+    if np.isclose(np.linalg.norm(basis1s[0]), 0):
+      basis1s[0] = np.cross(tangents[0], [1,0,0])
+    basis2s[0] = np.cross(tangents[0], basis1s[0])
+
+    basis1s[1] = basis1s[0]
+    basis2s[1] = basis2s[0]
+  else:
+    basis1s[0] = np.cross(tangents[0], tangents[1])
+    basis2s[0] = np.cross(basis1s[0] , tangents[0])
+
+    for i in range(1, len(tangents)):
+      R = _rotation_matrix_from_vectors(tangents[i-1], tangents[i])
+      basis1s[i] = R @ basis1s[i-1].T
+      basis2s[i] = R @ basis2s[i-1].T
+
+  for i in range(len(basis1s)):
+    basis1s[i] /= np.linalg.norm(basis1s[i])
+    basis2s[i] /= np.linalg.norm(basis2s[i])
+
+  for pos, basis1, basis2 in zip(path, basis1s, basis2s):
+      yield fastxs3d.projection_with_frame(
+        labels, pos, 
+        basis1, basis2, 
+        anisotropy
+      )
+
+# From SO: https://stackoverflow.com/questions/14313510/how-to-calculate-rolling-moving-average-using-python-numpy-scipy
+def _moving_average(a:np.ndarray, n:int, mode:str = "symmetric") -> np.ndarray:
+  if n <= 0:
+    raise ValueError(f"Window size ({n}), must be >= 1.")
+  elif n == 1:
+    return a
+
+  if len(a) == 0:
+    return a
+
+  if a.ndim == 2:
+    a = np.pad(a, [[n, n],[0,0]], mode=mode)
+  else:
+    a = np.pad(a, [n, n], mode=mode)
 
+  ret = np.cumsum(a, dtype=float, axis=0)
+  ret = (ret[n:] - ret[:-n])[:-n]
+  ret /= float(n)
+  return ret
 
+def _rotation_matrix_from_vectors(vec1, vec2):
+  """
+  Find the rotation matrix that aligns vec1 to vec2
+  :param vec1: A 3d "source" vector
+  :param vec2: A 3d "destination" vector
+  :return mat: A transformation matrix (3x3) which when applied to vec1, aligns it with vec2.
+
+  Credit: help from chatgpt
+  """
+  a, b = (vec1 / np.linalg.norm(vec1)).reshape(3), (vec2 / np.linalg.norm(vec2)).reshape(3)
+  v = np.cross(a, b)
+  c = np.dot(a, b)
+  s = np.linalg.norm(v)
+  if s == 0:
+    return np.eye(3)
+
+  kmat = np.array([
+    [0, -v[2], v[1]], 
+    [v[2], 0, -v[0]], 
+    [-v[1], v[0], 0]
+  ])
+  return np.eye(3) + kmat + kmat @ kmat * ((1 - c) / (s * s))