Overhaul cluster algorithm (3d extension needed) (#89)

docs/examples/clusters.py

@@ -22,6 +22,7 @@
     parser.add_argument("--save", action="store_true")
     parser.add_argument("--check", action="store_true")
     parser.add_argument("--plot", action="store_true")
+    parser.add_argument("--show", action="store_true")
     args = parser.parse_args()
 
     if args.save:
@@ -44,6 +45,7 @@
         import matplotlib.pyplot as plt
         import matplotlib as mpl
         import matplotlib.cm as cm
+        import prrng
         from mpl_toolkits.axes_grid1 import make_axes_locatable
 
         # color-scheme: modify such that the background is white
@@ -52,12 +54,30 @@
         cmap[0, :3] = 1.0
         cmap = mpl.colors.ListedColormap(cmap)
 
+        # reshuffle for better visualisation
+        assert np.all(np.diff(np.unique(clusters)) == 1)
+        assert np.all(np.diff(np.unique(clusters_periodic)) == 1)
+        assert np.unique(clusters).size >= np.unique(clusters_periodic).size
+        rng = prrng.pcg32()
+        lab = np.unique(clusters)
+        lab = lab[lab != 0]
+        new = np.copy(lab).astype(np.int64)
+        rng.shuffle(new)
+        rename = np.array(([0] + list(lab), [0] + list(new))).T
+        clusters = GooseEYE.labels_rename(clusters, rename)
+        lmap = GooseEYE.labels_map(clusters_periodic, clusters)
+        unq, unq_idx, unq_cnt = np.unique(lmap[:, 1], return_inverse=True, return_counts=True)
+        assert np.all(np.in1d(np.unique(clusters_periodic), lmap[:, 0]))
+        assert np.all(np.equal(np.sort(lmap[:, 0]), np.unique(lmap[:, 0])))
+        assert np.all(np.equal(np.sort(lmap[:, 1]), np.unique(lmap[:, 1])))
+        clusters_periodic = GooseEYE.labels_rename(clusters_periodic, lmap)
+
         try:
             plt.style.use(["goose", "goose-latex"])
         except OSError:
             pass
 
-        fig, axes = plt.subplots(figsize=(18, 6), nrows=1, ncols=3)
+        fig, axes = plt.subplots(figsize=(8 * 4, 6), nrows=1, ncols=4)
 
         ax = axes[0]
         im = ax.imshow(img, clim=(0, 1), cmap=mpl.colors.ListedColormap(cm.gray([0, 255])))
@@ -82,24 +102,33 @@
         ax.set_xlabel(r"$x$")
         ax.set_ylabel(r"$y$")
         ax.set_title(r"clusters")
-        div = make_axes_locatable(ax)
-        cax = div.append_axes("right", size="5%", pad=0.1)
-        cbar = plt.colorbar(im, cax=cax)
-        cbar.set_ticks([])
 
         ax = axes[2]
-        im = ax.imshow(clusters_periodic, clim=(0, np.max(clusters) + 1), cmap=cmap)
+        im = ax.imshow(clusters_periodic, clim=(0, np.max(clusters_periodic) + 1), cmap=cmap)
         ax.xaxis.set_ticks([0, 500])
         ax.yaxis.set_ticks([0, 500])
         ax.set_xlim([0, 500])
         ax.set_ylim([0, 500])
         ax.set_xlabel(r"$x$")
         ax.set_ylabel(r"$y$")
         ax.set_title(r"clusters (periodic)")
-        div = make_axes_locatable(ax)
-        cax = div.append_axes("right", size="5%", pad=0.1)
-        cbar = plt.colorbar(im, cax=cax)
-        cbar.set_ticks([])
 
-        fig.savefig(root / "clusters.svg")
+        ax = axes[3]
+        im = ax.imshow(
+            np.where(clusters_periodic != clusters, clusters_periodic, 0),
+            clim=(0, np.max(clusters_periodic) + 1),
+            cmap=cmap,
+        )
+        ax.xaxis.set_ticks([0, 500])
+        ax.yaxis.set_ticks([0, 500])
+        ax.set_xlim([0, 500])
+        ax.set_ylim([0, 500])
+        ax.set_xlabel(r"$x$")
+        ax.set_ylabel(r"$y$")
+        ax.set_title(r"difference")
+
+        if args.show:
+            plt.show()
+        else:
+            fig.savefig(root / "clusters.svg")
         plt.close(fig)

environment.yaml

@@ -3,7 +3,6 @@ channels:
 dependencies:
 - catch2
 - cmake
-- docopt
 - doxygen
 - enstat >=0.9.7
 - h5py
@@ -16,6 +15,7 @@ dependencies:
 - python
 - python-prrng
 - scikit-build
+- scipy
 - setuptools_scm
 - xsimd
 - xtensor

include/GooseEYE/Ensemble.hpp

@@ -14,7 +14,7 @@ namespace GooseEYE {
 inline Ensemble::Ensemble(const std::vector<size_t>& roi, bool periodic, bool variance)
     : m_periodic(periodic), m_variance(variance), m_shape_orig(roi)
 {
-    GOOSEEYE_ASSERT(m_shape_orig.size() <= MAX_DIM);
+    GOOSEEYE_ASSERT(m_shape_orig.size() <= MAX_DIM, std::out_of_range);
 
     m_first = xt::atleast_3d(xt::zeros<double>(m_shape_orig));
     m_second = zeros_like(m_first);
@@ -70,7 +70,7 @@ inline array_type::array<double> Ensemble::norm() const
 
 inline array_type::array<double> Ensemble::distance(size_t axis) const
 {
-    GOOSEEYE_ASSERT(axis < m_shape_orig.size());
+    GOOSEEYE_ASSERT(axis < m_shape_orig.size(), std::out_of_range);
     axis = detail::atleast_3d_axis(m_shape_orig.size(), axis);
 
     array_type::tensor<double, 3> dist = xt::empty<double>(m_shape);
@@ -113,7 +113,7 @@ inline array_type::array<double> Ensemble::distance() const
 
 inline array_type::array<double> Ensemble::distance(const std::vector<double>& h) const
 {
-    GOOSEEYE_ASSERT(m_shape_orig.size() == h.size());
+    GOOSEEYE_ASSERT(m_shape_orig.size() == h.size(), std::out_of_range);
 
     array_type::array<double> ret = xt::zeros<double>(m_shape_orig);
 
@@ -126,7 +126,7 @@ inline array_type::array<double> Ensemble::distance(const std::vector<double>& h
 
 inline array_type::array<double> Ensemble::distance(const std::vector<double>& h, size_t axis) const
 {
-    GOOSEEYE_ASSERT(m_shape_orig.size() == h.size());
+    GOOSEEYE_ASSERT(m_shape_orig.size() == h.size(), std::out_of_range);
     return this->distance(axis) * h[axis];
 }
 

include/GooseEYE/Ensemble_C2.hpp

@@ -20,13 +20,13 @@ inline void Ensemble::C2(const T& f, const T& g, const M& fmask, const M& gmask)
     static_assert(std::is_integral<value_type>::value, "Integral image required.");
     static_assert(std::is_integral<mask_type>::value, "Integral mask required.");
 
-    GOOSEEYE_ASSERT(xt::has_shape(f, g.shape()));
-    GOOSEEYE_ASSERT(xt::has_shape(f, fmask.shape()));
-    GOOSEEYE_ASSERT(xt::has_shape(f, gmask.shape()));
-    GOOSEEYE_ASSERT(f.dimension() == m_shape_orig.size());
-    GOOSEEYE_ASSERT(xt::all(xt::equal(fmask, 0) || xt::equal(fmask, 1)));
-    GOOSEEYE_ASSERT(xt::all(xt::equal(gmask, 0) || xt::equal(gmask, 1)));
-    GOOSEEYE_ASSERT(m_stat == Type::C2 || m_stat == Type::Unset);
+    GOOSEEYE_ASSERT(xt::has_shape(f, g.shape()), std::out_of_range);
+    GOOSEEYE_ASSERT(xt::has_shape(f, fmask.shape()), std::out_of_range);
+    GOOSEEYE_ASSERT(xt::has_shape(f, gmask.shape()), std::out_of_range);
+    GOOSEEYE_ASSERT(f.dimension() == m_shape_orig.size(), std::out_of_range);
+    GOOSEEYE_ASSERT(xt::all(xt::equal(fmask, 0) || xt::equal(fmask, 1)), std::out_of_range);
+    GOOSEEYE_ASSERT(xt::all(xt::equal(gmask, 0) || xt::equal(gmask, 1)), std::out_of_range);
+    GOOSEEYE_ASSERT(m_stat == Type::C2 || m_stat == Type::Unset, std::out_of_range);
 
     // lock statistic
     m_stat = Type::C2;

include/GooseEYE/Ensemble_L.hpp

@@ -18,8 +18,8 @@ inline void Ensemble::L(const T& f, path_mode mode)
 
     static_assert(std::is_integral<value_type>::value, "Integral image required.");
 
-    GOOSEEYE_ASSERT(f.dimension() == m_shape_orig.size());
-    GOOSEEYE_ASSERT(m_stat == Type::L || m_stat == Type::Unset);
+    GOOSEEYE_ASSERT(f.dimension() == m_shape_orig.size(), std::out_of_range);
+    GOOSEEYE_ASSERT(m_stat == Type::L || m_stat == Type::Unset, std::out_of_range);
 
     // lock statistics
     m_stat = Type::L;

include/GooseEYE/Ensemble_S2.hpp

@@ -19,13 +19,13 @@ inline void Ensemble::S2(const T& f, const T& g, const M& fmask, const M& gmask)
 
     static_assert(std::is_integral<mask_type>::value, "Integral mask required.");
 
-    GOOSEEYE_ASSERT(xt::has_shape(f, g.shape()));
-    GOOSEEYE_ASSERT(xt::has_shape(f, fmask.shape()));
-    GOOSEEYE_ASSERT(xt::has_shape(f, gmask.shape()));
-    GOOSEEYE_ASSERT(f.dimension() == m_shape_orig.size());
-    GOOSEEYE_ASSERT(xt::all(xt::equal(fmask, 0) || xt::equal(fmask, 1)));
-    GOOSEEYE_ASSERT(xt::all(xt::equal(gmask, 0) || xt::equal(gmask, 1)));
-    GOOSEEYE_ASSERT(m_stat == Type::S2 || m_stat == Type::Unset);
+    GOOSEEYE_ASSERT(xt::has_shape(f, g.shape()), std::out_of_range);
+    GOOSEEYE_ASSERT(xt::has_shape(f, fmask.shape()), std::out_of_range);
+    GOOSEEYE_ASSERT(xt::has_shape(f, gmask.shape()), std::out_of_range);
+    GOOSEEYE_ASSERT(f.dimension() == m_shape_orig.size(), std::out_of_range);
+    GOOSEEYE_ASSERT(xt::all(xt::equal(fmask, 0) || xt::equal(fmask, 1)), std::out_of_range);
+    GOOSEEYE_ASSERT(xt::all(xt::equal(gmask, 0) || xt::equal(gmask, 1)), std::out_of_range);
+    GOOSEEYE_ASSERT(m_stat == Type::S2 || m_stat == Type::Unset, std::out_of_range);
 
     // lock statistic
     m_stat = Type::S2;

include/GooseEYE/Ensemble_W2.hpp

@@ -19,11 +19,11 @@ inline void Ensemble::W2(const T& f, const T& g, const M& gmask)
 
     static_assert(std::is_integral<mask_type>::value, "Integral mask required.");
 
-    GOOSEEYE_ASSERT(xt::has_shape(f, g.shape()));
-    GOOSEEYE_ASSERT(xt::has_shape(f, gmask.shape()));
-    GOOSEEYE_ASSERT(f.dimension() == m_shape_orig.size());
-    GOOSEEYE_ASSERT(xt::all(xt::equal(gmask, 0) || xt::equal(gmask, 1)));
-    GOOSEEYE_ASSERT(m_stat == Type::W2 || m_stat == Type::Unset);
+    GOOSEEYE_ASSERT(xt::has_shape(f, g.shape()), std::out_of_range);
+    GOOSEEYE_ASSERT(xt::has_shape(f, gmask.shape()), std::out_of_range);
+    GOOSEEYE_ASSERT(f.dimension() == m_shape_orig.size(), std::out_of_range);
+    GOOSEEYE_ASSERT(xt::all(xt::equal(gmask, 0) || xt::equal(gmask, 1)), std::out_of_range);
+    GOOSEEYE_ASSERT(m_stat == Type::W2 || m_stat == Type::Unset, std::out_of_range);
 
     // lock statistic
     m_stat = Type::W2;

include/GooseEYE/Ensemble_W2c.hpp

@@ -22,12 +22,12 @@ Ensemble::W2c(const C& clusters, const C& centers, const T& f, const M& fmask, p
     static_assert(std::is_integral<cluster_type>::value, "Integral clusters required.");
     static_assert(std::is_integral<mask_type>::value, "Integral mask required.");
 
-    GOOSEEYE_ASSERT(f.shape() == clusters.shape());
-    GOOSEEYE_ASSERT(f.shape() == centers.shape());
-    GOOSEEYE_ASSERT(xt::has_shape(f, fmask.shape()));
-    GOOSEEYE_ASSERT(f.dimension() == m_shape_orig.size());
-    GOOSEEYE_ASSERT(xt::all(xt::equal(fmask, 0) || xt::equal(fmask, 1)));
-    GOOSEEYE_ASSERT(m_stat == Type::W2c || m_stat == Type::Unset);
+    GOOSEEYE_ASSERT(f.shape() == clusters.shape(), std::out_of_range);
+    GOOSEEYE_ASSERT(f.shape() == centers.shape(), std::out_of_range);
+    GOOSEEYE_ASSERT(xt::has_shape(f, fmask.shape()), std::out_of_range);
+    GOOSEEYE_ASSERT(f.dimension() == m_shape_orig.size(), std::out_of_range);
+    GOOSEEYE_ASSERT(xt::all(xt::equal(fmask, 0) || xt::equal(fmask, 1)), std::out_of_range);
+    GOOSEEYE_ASSERT(m_stat == Type::W2c || m_stat == Type::Unset, std::out_of_range);
 
     // lock statistic
     m_stat = Type::W2c;

include/GooseEYE/Ensemble_heightheight.hpp

@@ -18,10 +18,10 @@ inline void Ensemble::heightheight(const T& f, const M& fmask)
 
     static_assert(std::is_integral<mask_type>::value, "Integral mask required.");
 
-    GOOSEEYE_ASSERT(xt::has_shape(f, fmask.shape()));
-    GOOSEEYE_ASSERT(f.dimension() == m_shape_orig.size());
-    GOOSEEYE_ASSERT(xt::all(xt::equal(fmask, 0) || xt::equal(fmask, 1)));
-    GOOSEEYE_ASSERT(m_stat == Type::heightheight || m_stat == Type::Unset);
+    GOOSEEYE_ASSERT(xt::has_shape(f, fmask.shape()), std::out_of_range);
+    GOOSEEYE_ASSERT(f.dimension() == m_shape_orig.size(), std::out_of_range);
+    GOOSEEYE_ASSERT(xt::all(xt::equal(fmask, 0) || xt::equal(fmask, 1)), std::out_of_range);
+    GOOSEEYE_ASSERT(m_stat == Type::heightheight || m_stat == Type::Unset, std::out_of_range);
 
     // lock statistic
     m_stat = Type::heightheight;

include/GooseEYE/Ensemble_mean.hpp

@@ -14,8 +14,8 @@ namespace GooseEYE {
 template <class T>
 void Ensemble::mean(const T& f)
 {
-    GOOSEEYE_ASSERT(m_shape == std::vector<size_t>(MAX_DIM, 1));
-    GOOSEEYE_ASSERT(m_stat == Type::mean || m_stat == Type::Unset);
+    GOOSEEYE_ASSERT(m_shape == std::vector<size_t>(MAX_DIM, 1), std::out_of_range);
+    GOOSEEYE_ASSERT(m_stat == Type::mean || m_stat == Type::Unset, std::out_of_range);
 
     m_stat = Type::mean;
 
@@ -29,10 +29,10 @@ void Ensemble::mean(const T& f, const M& fmask)
 {
     static_assert(std::is_integral<typename M::value_type>::value, "Integral mask required.");
 
-    GOOSEEYE_ASSERT(xt::has_shape(f, fmask.shape()));
-    GOOSEEYE_ASSERT(xt::all(xt::equal(fmask, 0) || xt::equal(fmask, 1)));
-    GOOSEEYE_ASSERT(m_shape == std::vector<size_t>(MAX_DIM, 1));
-    GOOSEEYE_ASSERT(m_stat == Type::mean || m_stat == Type::Unset);
+    GOOSEEYE_ASSERT(xt::has_shape(f, fmask.shape()), std::out_of_range);
+    GOOSEEYE_ASSERT(xt::all(xt::equal(fmask, 0) || xt::equal(fmask, 1)), std::out_of_range);
+    GOOSEEYE_ASSERT(m_shape == std::vector<size_t>(MAX_DIM, 1), std::out_of_range);
+    GOOSEEYE_ASSERT(m_stat == Type::mean || m_stat == Type::Unset, std::out_of_range);
 
     m_stat = Type::mean;