Planar index against fragmentation

src/simulation.rs

@@ -17,6 +17,22 @@ pub(crate) const COUPLING_NEIGHBORHOOD: [[i64; DIM]; 26] = [
     [1, 1, -1], [1, 1, 0], [1, 1, 1]
 ];
 
+/// Directions of adjacency plane
+pub(crate) type AdjacencyPlane = [[i64; DIM]; 4];
+
+/// Adjacency plane across X
+pub(crate) const ADJACENCY_X: AdjacencyPlane = [
+    [0, 0, 1], [0, 0, -1], [0, 1, 0], [0, -1, 0]
+];
+/// Adjacency plane across Y
+pub(crate) const ADJACENCY_Y: AdjacencyPlane = [
+    [0, 0, 1], [0, 0, -1], [1, 0, 0], [-1, 0, 0]
+];
+/// Adjacency plane across Z
+pub(crate) const ADJACENCY_Z: AdjacencyPlane = [
+    [0, 1, 0], [0, -1, 0], [1, 0, 0], [-1, 0, 0]
+];
+
 /// A stable implementation of the softmax function with multiplication in place
 pub(crate) fn softmax_(slice: &mut [f64]) {
     let max = slice.iter().fold(0f64, |max, x| if *x > max { *x } else { max });
@@ -210,8 +226,43 @@ impl Kernel {
         ).powi(2) / 2f64
     }
 
+    /// Count voxels in adjacency plane
+    fn plane_index(&self, voxel: &Voxel, voxel_value: usize, adjacency: &AdjacencyPlane) -> usize {
+        let mut result = 0;
+
+        for triplet in adjacency {
+            if voxel_value == self.get_index_value(
+                voxel.index[0] + triplet[0],
+                voxel.index[1] + triplet[1],
+                voxel.index[2] + triplet[2]
+            ) { result += 1 }
+        }
+
+        result
+    }
+    /// Calculate planar index
+    pub(crate) fn planar_index(&self, voxel: &Voxel, voxel_value: usize) -> [usize; DIM] {
+        let mut result = [0; DIM];
+
+        for (index, adjacency) in [&ADJACENCY_X, &ADJACENCY_Y, &ADJACENCY_Z].iter().enumerate() {
+            result[index] = self.plane_index(voxel, voxel_value, adjacency)
+        }
+        result.sort();
+
+        result
+    }
+
     /// Check whether local connectivity of a voxel allows a move
     pub(crate) fn check_change(&self, voxel: &Voxel, voxel_value: usize) -> bool {
+        // First screen planar index
+        let index = self.planar_index(voxel, voxel_value);
+        match index {
+            [0, 0, 0] | [0, 2, 2] | [2, 2, 4] | [4, 4, 4] => return false,
+            [0, 1, 1] => return true,
+            [1, 1, 2] | [1, 2, 3] | [2, 2, 2] | [2, 3, 3] | [3, 3, 4] => (),
+            _ => panic!("Unexpected planar index ({}, {}, {})", index[0], index[1], index[2])
+        }
+
         // The implementation relies on a graph connectivity check
         let mut graph_nodes = Vec::<[i64; DIM]>::with_capacity(5);
         for triplet in &TARGET_NEIGHBORHOOD {

src/tests.rs

@@ -405,3 +405,94 @@ fn connectivity_check() {
 
     assert!(!kernel.check_connectivity(&kernel.state.cells[&2]));
 }
+
+/// Test planar index
+#[test]
+fn planar_index() {
+    let mut kernel = Kernel::from_state(State {
+        side: 3, time: 0, cells: hashmap![
+            1usize => Cell {
+                type_id: EPI, preferred_volume: 2f64, division_volume: 2f64,
+                parent_id: 0usize, age: 0usize,
+                voxel_ids: hashset!{} // std::collections::HashSet::from([4i64, 10i64, 12i64, 14i64, 16i64, 22i64])
+            }
+        ],
+        action_probabilities: [0.001, 0.999, 0.999],
+        death_probability: 0f64,
+        growth_rate: [0.5, 0.5], division_distribution: vec![10f64],
+        volume_elasticity: [1.0, 1.5],
+        surface_tension: [1.0, 2.0, 3.0, 4.0, 5.0], modulation: [1.0, 1.0, 1.0],
+        rng: mt19937::MT19937::default()
+    });
+    let voxel = kernel.get_voxel_by_id(13i64);
+
+    // Case 1 (0)
+    kernel.lattice.copy_from_slice(&[0; 27]);
+    assert_eq!(kernel.planar_index(&voxel, 1usize), [0, 0, 0]);
+
+    // Case 2 (1)
+    kernel.lattice.copy_from_slice(&[0; 27]);
+    kernel.lattice[12] = 1;
+    assert_eq!(kernel.planar_index(&voxel, 1usize), [0, 1, 1]);
+
+    // Case 3 (2)
+    kernel.lattice.copy_from_slice(&[0; 27]);
+    kernel.lattice[10] = 1;
+    kernel.lattice[16] = 1;
+    assert_eq!(kernel.planar_index(&voxel, 1usize), [0, 2, 2]);
+
+    // Case 4 (2)
+    kernel.lattice.copy_from_slice(&[0; 27]);
+    kernel.lattice[12] = 1;
+    kernel.lattice[16] = 1;
+    assert_eq!(kernel.planar_index(&voxel, 1usize), [1, 1, 2]);
+
+    // Case 5 (3)
+    kernel.lattice.copy_from_slice(&[0; 27]);
+    kernel.lattice[10] = 1;
+    kernel.lattice[12] = 1;
+    kernel.lattice[16] = 1;
+    assert_eq!(kernel.planar_index(&voxel, 1usize), [1, 2, 3]);
+
+    // Case 6 (3)
+    kernel.lattice.copy_from_slice(&[0; 27]);
+    kernel.lattice[12] = 1;
+    kernel.lattice[16] = 1;
+    kernel.lattice[22] = 1;
+    assert_eq!(kernel.planar_index(&voxel, 1usize), [2, 2, 2]);
+
+    // Case 7 (4)
+    kernel.lattice.copy_from_slice(&[0; 27]);
+    kernel.lattice[10] = 1;
+    kernel.lattice[12] = 1;
+    kernel.lattice[16] = 1;
+    kernel.lattice[22] = 1;
+    assert_eq!(kernel.planar_index(&voxel, 1usize), [2, 3, 3]);
+
+    // Case 8 (4)
+    kernel.lattice.copy_from_slice(&[0; 27]);
+    kernel.lattice[10] = 1;
+    kernel.lattice[12] = 1;
+    kernel.lattice[14] = 1;
+    kernel.lattice[16] = 1;
+    assert_eq!(kernel.planar_index(&voxel, 1usize), [2, 2, 4]);
+
+    // Case 9 (5)
+    kernel.lattice.copy_from_slice(&[0; 27]);
+    kernel.lattice[10] = 1;
+    kernel.lattice[12] = 1;
+    kernel.lattice[14] = 1;
+    kernel.lattice[16] = 1;
+    kernel.lattice[22] = 1;
+    assert_eq!(kernel.planar_index(&voxel, 1usize), [3, 3, 4]);
+
+    // Case 10 (6)
+    kernel.lattice.copy_from_slice(&[0; 27]);
+    kernel.lattice[4] = 1;
+    kernel.lattice[10] = 1;
+    kernel.lattice[12] = 1;
+    kernel.lattice[14] = 1;
+    kernel.lattice[16] = 1;
+    kernel.lattice[22] = 1;
+    assert_eq!(kernel.planar_index(&voxel, 1usize), [4, 4, 4]);
+}
\ No newline at end of file