Merge pull request #114 from dyollb/feature/topology_simplification

add topological simplification of voxel set = "spherize"

Author: dyollb

Core/Graph.h

@@ -0,0 +1,153 @@
+/// \file FastMarching.h
+///
+/// \author Bryn Lloyd
+/// \copyright 2016/2019, IT'IS Foundation
+
+#pragma once
+
+#include <array>
+#include <cmath>
+
+namespace XCore
+{
+
+	/**	Uniform n-dim grid
+	*/
+	template<typename TIndex, int D>
+	class CUniformGridGraph
+	{
+		template <typename T> static constexpr T Pow(T num, unsigned pow) { return pow == 0 ? 1 : num * Pow(num, pow - 1); }
+
+		template<int N> struct type{};
+
+		template<typename T, int SIZE>
+		class StackVector
+		{
+		public:
+			using value_type = T;
+
+			const T* begin() const { return &_data[0]; }
+			const T* end() const { return begin() + _size; }
+
+			T* begin() { return &_data[0]; }
+			T* end() { return begin() + _size; }
+
+			inline T& operator[](size_t idx) { return _data[idx]; }
+			inline const T& operator[](size_t idx) const { return _data[idx]; }
+	
+			size_t size() const { return _size; }
+
+			void push_back(T v)
+			{
+				_data[_size] = v;
+				_size++;
+			}
+
+		private:
+			unsigned _size = 0;
+			std::array<T, SIZE> _data;
+		};
+
+	public:
+		using idx_type = TIndex;
+		using idx_value = typename idx_type::IndexValueType;
+		using idx_list_type = StackVector<idx_type, Pow(3, D) - 1>;
+
+		CUniformGridGraph(const std::array<idx_value, D>& dims, const std::array<double, D>& spacing)
+			: m_Dims(dims)
+		{
+			static_assert(std::is_signed<idx_value>::value, "idx_value should be signed");
+			static_assert(D > 0, "Dimension should be positive");
+
+			// total number of nodes in grid
+			m_Size = 1;
+			for (int k=0; k<D; ++k)
+			{
+				m_Strides[k] = m_Size;
+				m_Size *= dims[k];
+			}
+			// squared spacing
+			std::transform(spacing.begin(), spacing.end(), m_Spacing2.begin(), 
+				[](double s){ return s*s; });
+		}
+
+		size_t Size() const
+		{
+			return m_Size;
+		}
+
+		const std::array<idx_value, D>& Shape() const
+		{
+			return m_Dims;
+		}
+
+		const std::array<size_t, D>& Strides() const
+		{
+			return m_Strides;
+		}
+
+		double Length(const idx_type& a, const idx_type& b) const
+		{
+			return std::sqrt(LengthD<D-1>(type<D-1>(), a, b));
+		}
+
+		idx_list_type Neighbors(const idx_type& idx) const
+		{
+			idx_list_type neighbors;
+			NeighborsD<D-1>(idx, false, neighbors);
+			return neighbors;
+		}
+
+	private:
+		template<int N>
+		inline double LengthD(type<N>, const idx_type& a, const idx_type& b) const
+		{
+			return (a[N] - b[N])*(a[N] - b[N])*m_Spacing2[N] + LengthD<N - 1>(a, b);
+		}
+
+		inline double LengthD(type<-1>, const idx_type& a, const idx_type& b) const
+		{
+			return 0;
+		}
+
+		template<int N>
+		inline void NeighborsD(const idx_type& idx, bool collect, idx_list_type& collected) const
+		{
+			NeighborsD(type<N>(), idx, collect, collected);
+		}
+
+		template<int N>
+		inline void NeighborsD(type<N>, const idx_type& idx, bool collect, idx_list_type& collected) const
+		{
+			// make modifiable copy
+			auto ref{ idx };
+
+			// center 
+			NeighborsD<N - 1>(ref, collect, collected);
+			// bottom
+			if (idx[N] > 0)
+			{
+				ref[N] = idx[N] - 1;
+				NeighborsD<N - 1>(ref, true, collected);
+			}
+			// top
+			if (idx[N] + 1 < m_Dims[N])
+			{
+				ref[N] = idx[N] + 1;
+				NeighborsD<N - 1>(ref, true, collected);
+			}
+		}
+
+		void NeighborsD(type<-1>, const idx_type& idx, bool collect, idx_list_type& collected) const
+		{
+			if (collect)
+				collected.push_back(idx);
+		}
+
+		std::array<idx_value, D> m_Dims;
+		std::array<double, D> m_Spacing2;
+		size_t m_Size;
+		std::array<size_t, D> m_Strides;
+	};
+
+}

Core/TopologyInvariants.h

@@ -0,0 +1,214 @@
+/// \file TopologyInvariants.h
+///
+/// \author Bryn Lloyd
+/// \copyright 2020, IT'IS Foundation
+
+#pragma once
+
+#include <array>
+#include <deque>
+
+namespace topology {
+
+template<typename T>
+inline int CountIf(const T& v)
+{
+	static_assert(std::is_same<T,bool>::value, "CountIf argument must be a boolean.");
+	return (v != 0) ? 1 : 0;
+}
+
+/** \brief True if Euler characteristic of voxel set does not change
+*/
+template<typename TNeighborhood, typename TLabel>
+bool EulerInvariant(const TNeighborhood& neighbors, const TLabel label)
+{
+	// For a face-connected set, a vertex, edge, or facet is included 
+	// only if all voxels it is incident to are in the set, and for a 
+	// vertex-connected set it is included if any incident voxels are in the set
+	
+	// Here we count only the changes at the center voxel.
+	// ==> If we mark the center as BG, then P, F, E and V are '0'.
+	static const int c = 27 / 2;
+
+//fprintf(stderr, "neighbors[c]: %d vs %d", static_cast<int>(neighbors[c]), static_cast<int>(label));
+	assert(neighbors[c] == label);
+
+	// Voxels (parallelepipeds) - count center
+	const int P = 1;
+
+	// Faces - count faces around center voxel
+	const int F = CountIf(neighbors[c - 1] == label) // -x
+			+ CountIf(neighbors[c + 1] == label) // +x
+			+ CountIf(neighbors[c - 3] == label) // -y
+			+ CountIf(neighbors[c + 3] == label) // +y
+			+ CountIf(neighbors[c - 9] == label) // -z
+			+ CountIf(neighbors[c + 9] == label); // +z
+
+	// Edges - count edges around center voxel
+	int E = 0;
+	{
+		static const short offsets[12][2] = {
+				{-1, -3},
+				{+1, -3},
+				{-1, +3},
+				{+1, +3},
+				{-1, -9},
+				{+1, -9},
+				{-1, +9},
+				{+1, +9},
+				{-3, -9},
+				{+3, -9},
+				{-3, +9},
+				{+3, +9}};
+		for (int i = 0; i < 12; i++)
+		{
+			short o = offsets[i][0];
+			short u = offsets[i][1];
+
+			E += CountIf(
+					neighbors[c + o] == label &&
+					neighbors[c + u] == label &&
+					neighbors[c + o + u] == label);
+		}
+	}
+
+	// Vertices - count verts around center voxel
+	int V = 0;
+	{
+		static const short offsets[8][3] = {
+				{-1, -3, -9},
+				{+1, -3, -9},
+				{-1, +3, -9},
+				{+1, +3, -9},
+				{-1, -3, +9},
+				{+1, -3, +9},
+				{-1, +3, +9},
+				{+1, +3, +9}};
+		for (int i = 0; i < 8; i++)
+		{
+			short x = offsets[i][0];
+			short y = offsets[i][1];
+			short z = offsets[i][2];
+
+			V += CountIf(
+					neighbors[c + x] == label &&
+					neighbors[c + y] == label &&
+					neighbors[c + z] == label &&
+					neighbors[c + x + y] == label &&
+					neighbors[c + x + z] == label &&
+					neighbors[c + y + z] == label &&
+					neighbors[c + x + y + z] == label);
+		}
+	}
+
+	return (V - E + F - P) == 0;
+}
+
+/** \brief Returns the number of connected components for selected label
+*/
+template<typename TNeighborhood, typename TLabel>
+unsigned ConnectedComponents(const TNeighborhood& neighbors, const TLabel label)
+{
+	std::deque<unsigned> queue;
+	std::array<bool, 27> visited;
+	visited.fill(false);
+
+	static const short ijk_lut[27][3] = {
+			{0, 0, 0},
+			{1, 0, 0},
+			{2, 0, 0},
+			{0, 1, 0},
+			{1, 1, 0},
+			{2, 1, 0},
+			{0, 2, 0},
+			{1, 2, 0},
+			{2, 2, 0},
+			//
+			{0, 0, 1},
+			{1, 0, 1},
+			{2, 0, 1},
+			{0, 1, 1},
+			{1, 1, 1},
+			{2, 1, 1},
+			{0, 2, 1},
+			{1, 2, 1},
+			{2, 2, 1},
+			//
+			{0, 0, 2},
+			{1, 0, 2},
+			{2, 0, 2},
+			{0, 1, 2},
+			{1, 1, 2},
+			{2, 1, 2},
+			{0, 2, 2},
+			{1, 2, 2},
+			{2, 2, 2},
+	};
+
+	unsigned new_label = 0;
+	for (unsigned n = 0; n < 27; ++n)
+	{
+		// skip other labels and visited voxels
+		if (neighbors[n] != label || visited[n])
+			continue;
+
+		queue.push_back(n);
+		visited[n] = true;
+		new_label++;
+
+		while (!queue.empty())
+		{
+			auto id = queue.front();
+			queue.pop_front();
+			auto ijk = ijk_lut[id];
+
+			if (ijk[0] < 2 && neighbors[id] == neighbors[id + 1] && !visited[id + 1]) // +x
+			{
+				queue.push_back(id + 1);
+				visited[id + 1] = true;
+			}
+			if (ijk[0] > 0 && neighbors[id] == neighbors[id - 1] && !visited[id - 1]) // -x
+			{
+				queue.push_back(id - 1);
+				visited[id - 1] = true;
+			}
+			if (ijk[1] < 2 && neighbors[id] == neighbors[id + 3] && !visited[id + 3]) // +y
+			{
+				queue.push_back(id + 3);
+				visited[id + 3] = true;
+			}
+			if (ijk[1] > 0 && neighbors[id] == neighbors[id - 3] && !visited[id - 3]) // -y
+			{
+				queue.push_back(id - 3);
+				visited[id - 3] = true;
+			}
+			if (ijk[2] < 2 && neighbors[id] == neighbors[id + 9] && !visited[id + 9]) // +z
+			{
+				queue.push_back(id + 9);
+				visited[id + 9] = true;
+			}
+			if (ijk[2] > 0 && neighbors[id] == neighbors[id - 9] && !visited[id - 9]) // -z
+			{
+				queue.push_back(id - 9);
+				visited[id - 9] = true;
+			}
+		}
+	}
+
+	return new_label;
+}
+
+/** \brief True if number of connected components of voxel set does not change
+*/
+template<typename TNeighborhood, typename TLabel>
+bool CCInvariant(TNeighborhood neighbors, const TLabel label)
+{
+	neighbors[27 / 2] = 1;
+	unsigned cc_before = ConnectedComponents(neighbors, label);
+
+	neighbors[27 / 2] = 0;
+	unsigned cc_after = ConnectedComponents(neighbors, label);
+	return (cc_before == cc_after);
+}
+
+} // namespace topology