From 0c2e7fc8d4fae32acaeaf1c187b0b20eb4d33eac Mon Sep 17 00:00:00 2001
From: William Silversmith <william.silversmith@gmail.com>
Date: Tue, 25 Jun 2024 22:49:26 -0400
Subject: [PATCH 1/3] wip: experimenting with a smooth path feature

---
 src/fastxs3d.cpp |  90 +++++++++++++++++++++++++++++++++++++
 src/xs3d.hpp     |  54 +++++++++++++++++-----
 xs3d/__init__.py | 115 ++++++++++++++++++++++++++++++++++++++++++++++-
 3 files changed, 246 insertions(+), 13 deletions(-)

diff --git a/src/fastxs3d.cpp b/src/fastxs3d.cpp
index bfadda2..cced849 100644
--- a/src/fastxs3d.cpp
+++ b/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.");
diff --git a/src/xs3d.hpp b/src/xs3d.hpp
index c7af68f..7dd98bc 100644
--- a/src/xs3d.hpp
+++ b/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
diff --git a/xs3d/__init__.py b/xs3d/__init__.py
index a2cbead..c00493e 100644
--- a/xs3d/__init__.py
+++ b/xs3d/__init__.py
@@ -1,4 +1,4 @@
-from typing import Sequence, Optional
+from typing import Sequence, Optional, Tuple
 
 from .twod import cross_sectional_area_2d
 import fastxs3d
@@ -190,11 +190,122 @@ 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,
+) -> 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)
+  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.
+
+  Returns: ndarray
+  """
+  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)
+
+  # 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 = (tangents[1:] - tangents[:-1]).astype(np.float32)
+  basis1s = np.concatenate([ basis1s, [basis1s[-1]] ])
+
+  basis2s = []
+
+  basis1 = basis1s[0]
+  if np.all(basis1 == 0):
+    basis1 = np.cross(tangents[0], [1,0,0])
+    if np.all(basis1 == 0):
+      basis1 = np.cross(tangents[0], [0,1,0])
+
+  basis1s[0] = basis1
+
+  for i in range(1, len(basis1s)):
+    if np.all(basis1s[i] == 0):
+      basis1s[i] = basis1s[i-1]
+
+  basis2 = np.cross(tangents[0], basis1)
+  if np.all(basis2 == 0):
+    basis2 = np.cross(tangents[0], [1,0,0])
+    if np.all(basis2 == 0):
+      basis2 = np.cross(tangents[0], [0,1,0])
+
+  basis2s.append(basis2)
+
+  for tangent, delta in zip(tangents[1:], basis1s[1:]):
+    basis1 = delta
+    basis2 = np.cross(tangent, basis1)
+    if np.all(basis2 == 0):
+      basis2 = basis2s[-1]
+    basis2s.append(basis2)
+
+  for i in range(len(basis1s)):
+    basis1s[i] /= np.linalg.norm(basis1s[i])
+    basis2s[i] /= np.linalg.norm(basis2s[i])
+
+  slices = []
+  from tqdm import tqdm
+  for pos, basis1, basis2 in tqdm(zip(path, basis1s, basis2s)):
+    slices.append(
+      fastxs3d.projection_with_frame(
+        labels, pos, 
+        basis1, basis2, 
+        anisotropy
+      )
+    )
+  return slices
+
+# 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
 

From 39c0f2f64fdfa964647e7847c2be2d30abbfd6a6 Mon Sep 17 00:00:00 2001
From: William Silversmith <william.silversmith@gmail.com>
Date: Mon, 8 Jul 2024 18:24:34 -0400
Subject: [PATCH 2/3] wip: experimenting with rotating the basis vectors

---
 xs3d/__init__.py | 60 ++++++++++++++++++++++++++----------------------
 1 file changed, 32 insertions(+), 28 deletions(-)

diff --git a/xs3d/__init__.py b/xs3d/__init__.py
index c00493e..d7d3e8a 100644
--- a/xs3d/__init__.py
+++ b/xs3d/__init__.py
@@ -229,6 +229,29 @@ def slice_path(
   if labels.ndim != 3:
     raise ValueError(f"{labels.ndim} dimensions not supported")
 
+  # def degrees(v1,v2):
+  #   cos = np.dot(v1, v2) / np.linalg.norm(v1) / np.linalg.norm(v2)
+  #   cos = np.clip(cos, -1, 1)
+  #   return np.arccos(cos) / np.pi / 2.0 * 360.0
+
+  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.
+    """
+    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]])
+    rotation_matrix = np.eye(3) + kmat + kmat @ kmat * ((1 - c) / (s ** 2))
+    return rotation_matrix
+
   anisotropy = np.array(anisotropy, dtype=np.float32)
   labels = np.asfortranarray(labels)
 
@@ -241,37 +264,18 @@ def slice_path(
   tangents = _moving_average(tangents, smoothing)
   tangents = _moving_average(tangents[::-1], smoothing)[::-1]
 
-  basis1s = (tangents[1:] - tangents[:-1]).astype(np.float32)
-  basis1s = np.concatenate([ basis1s, [basis1s[-1]] ])
-
-  basis2s = []
-
-  basis1 = basis1s[0]
-  if np.all(basis1 == 0):
-    basis1 = np.cross(tangents[0], [1,0,0])
-    if np.all(basis1 == 0):
-      basis1 = np.cross(tangents[0], [0,1,0])
-
-  basis1s[0] = basis1
-
-  for i in range(1, len(basis1s)):
-    if np.all(basis1s[i] == 0):
-      basis1s[i] = basis1s[i-1]
+  # need to consider single tangent too
 
-  basis2 = np.cross(tangents[0], basis1)
-  if np.all(basis2 == 0):
-    basis2 = np.cross(tangents[0], [1,0,0])
-    if np.all(basis2 == 0):
-      basis2 = np.cross(tangents[0], [0,1,0])
+  basis1s = np.zeros([ len(tangents), 3 ], dtype=np.float32)
+  basis2s = np.zeros([ len(tangents), 3 ], dtype=np.float32)
 
-  basis2s.append(basis2)
+  basis1s[0] = np.cross(tangents[0], tangents[1])
+  basis2s[0] = np.cross(basis1s[0] , tangents[0])
 
-  for tangent, delta in zip(tangents[1:], basis1s[1:]):
-    basis1 = delta
-    basis2 = np.cross(tangent, basis1)
-    if np.all(basis2 == 0):
-      basis2 = basis2s[-1]
-    basis2s.append(basis2)
+  for i in range(1, len(tangents)):
+    R = rotation_matrix_from_vectors(tangents[i-1], tangents[i])
+    basis1s[i] = np.matmul(R, basis1s[i-1].T)
+    basis2s[i] = np.matmul(R, basis2s[i-1].T)
 
   for i in range(len(basis1s)):
     basis1s[i] /= np.linalg.norm(basis1s[i])

From 7ad93093b65e11980ed33c7c7da93ad8b6330a81 Mon Sep 17 00:00:00 2001
From: William Silversmith <william.silversmith@gmail.com>
Date: Tue, 9 Jul 2024 20:47:42 -0400
Subject: [PATCH 3/3] fix: cover edge cases, add type annotation, generator

---
 xs3d/__init__.py | 98 +++++++++++++++++++++++++++---------------------
 1 file changed, 55 insertions(+), 43 deletions(-)

diff --git a/xs3d/__init__.py b/xs3d/__init__.py
index d7d3e8a..91e050b 100644
--- a/xs3d/__init__.py
+++ b/xs3d/__init__.py
@@ -1,4 +1,4 @@
-from typing import Sequence, Optional, Tuple
+from typing import Sequence, Optional, Tuple, Generator
 
 from .twod import cross_sectional_area_2d
 import fastxs3d
@@ -195,11 +195,16 @@ def slice_path(
   path:Sequence[Sequence[int]],
   anisotropy:Optional[Sequence[float]] = None,
   smoothing:int = 1,
-) -> np.ndarray:
+) -> 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.
@@ -212,11 +217,8 @@ def slice_path(
     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.
-
-  Returns: ndarray
   """
   if anisotropy is None:
     anisotropy = [ 1.0 ] * labels.ndim
@@ -229,32 +231,18 @@ def slice_path(
   if labels.ndim != 3:
     raise ValueError(f"{labels.ndim} dimensions not supported")
 
-  # def degrees(v1,v2):
-  #   cos = np.dot(v1, v2) / np.linalg.norm(v1) / np.linalg.norm(v2)
-  #   cos = np.clip(cos, -1, 1)
-  #   return np.arccos(cos) / np.pi / 2.0 * 360.0
-
-  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.
-    """
-    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]])
-    rotation_matrix = np.eye(3) + kmat + kmat @ kmat * ((1 - c) / (s ** 2))
-    return rotation_matrix
-
   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]] ])
@@ -264,34 +252,36 @@ def rotation_matrix_from_vectors(vec1, vec2):
   tangents = _moving_average(tangents, smoothing)
   tangents = _moving_average(tangents[::-1], smoothing)[::-1]
 
-  # need to consider single tangent too
-
   basis1s = np.zeros([ len(tangents), 3 ], dtype=np.float32)
   basis2s = np.zeros([ len(tangents), 3 ], dtype=np.float32)
 
-  basis1s[0] = np.cross(tangents[0], tangents[1])
-  basis2s[0] = np.cross(basis1s[0] , tangents[0])
+  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] = np.matmul(R, basis1s[i-1].T)
-    basis2s[i] = np.matmul(R, basis2s[i-1].T)
+    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])
 
-  slices = []
-  from tqdm import tqdm
-  for pos, basis1, basis2 in tqdm(zip(path, basis1s, basis2s)):
-    slices.append(
-      fastxs3d.projection_with_frame(
+  for pos, basis1, basis2 in zip(path, basis1s, basis2s):
+      yield fastxs3d.projection_with_frame(
         labels, pos, 
         basis1, basis2, 
         anisotropy
       )
-    )
-  return slices
 
 # 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:
@@ -313,3 +303,25 @@ def _moving_average(a:np.ndarray, n:int, mode:str = "symmetric") -> np.ndarray:
   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))