Merge pull request #1 from mino2357/kd-tree

Kd treeの実装

.gitignore

@@ -14,5 +14,6 @@ Cargo.lock
 *.pdb
 
 *.png
+*.log
 
 !images/*.png

Cargo.toml

@@ -1,5 +1,5 @@
 [package]
-name = "scrapbox"
+name = "kd-tree"
 version = "0.1.0"
 edition = "2021"
 

README.md

@@ -10,4 +10,26 @@ $$
 
 ## 計算例
 
-![sim](images/sim.png)
+![sim](images/sim.png)
+
+## kd-木を用いた近傍探索
+
+`test` で用いた例。アルゴリズムは参考文献[1]。
+
+![](images/debug-01.jpg)
+
+![](images/debug-02.jpg)
+
+結晶粒界のような構造が見られる。
+
+![](images/grain_boundary.png)
+
+## kd-木の動作テスト
+
+$[-1,1]^2 \in \mathbb{R}^2$ 上に乱数を10000000点（1000万点）ふって動く点から半径0.25の距離未満にある点を緑で表示している。あまりに数が多いので画像が緑と赤で覆われているが点群から出来ている。
+
+![](images/10e7.png)
+
+## 参考文献
+
+1. [コンピュータ・ジオメトリ―計算幾何学:アルゴリズムと応用](https://www.kindaikagaku.co.jp/book_list/detail/9784764903883/)

src/gen_grid2d.rs

@@ -0,0 +1,182 @@
+use super::kd_tree;
+
+#[derive(Debug, Clone, PartialEq, PartialOrd)]
+pub struct Grid2D {
+    pub x: f64,
+    pub y: f64,
+}
+
+impl Grid2D {
+    #[allow(dead_code)]
+    pub fn new(x_: f64, y_: f64) -> Self {
+        Grid2D { x: x_, y: y_ }
+    }
+
+    #[allow(dead_code)]
+    pub fn distance_square(&self, point: &Grid2D) -> f64 {
+        let dx2 = (self.x - point.x).powf(2.0);
+        let dy2 = (self.y - point.y).powf(2.0);
+        dx2 + dy2
+    }
+}
+
+#[derive(Debug, Clone)]
+pub struct Points2D {
+    pub points: Vec<Grid2D>,
+}
+
+impl Points2D {
+    #[allow(dead_code)]
+    pub fn new() -> Self {
+        Points2D { points: vec![] }
+    }
+
+    #[allow(dead_code)]
+    pub fn push(&mut self, x_r: f64, y_r: f64) {
+        self.points.push(Grid2D { x: x_r, y: y_r });
+    }
+
+    #[allow(dead_code)]
+    pub fn euler_step_by_near_points(
+        &mut self,
+        boundary: &Points2D,
+        tree: &kd_tree::KDTree,
+        radius: f64,
+    ) {
+        let dt = 1.0e-3;
+        for i in 0..self.points.len() {
+            let mut x = self.points[i].x
+                - dt * self
+                    .lennard_jones_potential_deriv_by_near_points(i, boundary, tree, radius)
+                    .x;
+            let mut y = self.points[i].y
+                - dt * self
+                    .lennard_jones_potential_deriv_by_near_points(i, boundary, tree, radius)
+                    .y;
+            if x * x + y * y < 1.0 {
+                self.points[i].x = x;
+                self.points[i].y = y;
+            } else {
+                loop {
+                    let tmp = x;
+                    x = -0.5 * y;
+                    y = 0.5 * tmp;
+                    if x * x + y * y < 1.0 {
+                        self.points[i].x = x;
+                        self.points[i].y = y;
+                        break;
+                    }
+                }
+            }
+        }
+    }
+
+    #[allow(dead_code)]
+    pub fn euler_step(&mut self, boundary: &Points2D) {
+        let dt = 1.0e-3;
+        for i in 0..self.points.len() {
+            let mut x = self.points[i].x - dt * self.lennard_jones_potential_deriv(i, boundary).x;
+            let mut y = self.points[i].y - dt * self.lennard_jones_potential_deriv(i, boundary).y;
+            if x * x + y * y < 1.0 {
+                self.points[i].x = x;
+                self.points[i].y = y;
+            } else {
+                loop {
+                    let tmp = x;
+                    x = -0.5 * y;
+                    y = 0.5 * tmp;
+                    if x * x + y * y < 1.0 {
+                        self.points[i].x = x;
+                        self.points[i].y = y;
+                        break;
+                    }
+                }
+            }
+        }
+    }
+
+    #[allow(dead_code)]
+    pub fn lennard_jones_potential_deriv_by_near_points(
+        &mut self,
+        index: usize,
+        boundary: &Points2D,
+        tree: &kd_tree::KDTree,
+        radius: f64,
+    ) -> Grid2D {
+        let epsilon: f64 = 1.0;
+        let sigma: f64 = 0.9 * (2.0_f64).powf(-1.0 / 6.0) / (self.points.len() as f64).powf(0.5);
+        let mut f_x = 0.0;
+        let mut f_y = 0.0;
+        let near = tree.neighbor_search(&self.points[index], radius);
+        for k in near.iter() {
+            if index != *k {
+                let dx = self.points[index].x - self.points[*k].x;
+                let dy = self.points[index].y - self.points[*k].y;
+                let r = (dx.powf(2.0) + dy.powf(2.0)).powf(0.5);
+                f_x += 4.0
+                    * epsilon
+                    * (12.0 * sigma.powf(12.0) / r.powf(13.0)
+                        - 6.0 * sigma.powf(6.0) / r.powf(7.0))
+                    * dx;
+                f_y += 4.0
+                    * epsilon
+                    * (12.0 * sigma.powf(12.0) / r.powf(13.0)
+                        - 6.0 * sigma.powf(6.0) / r.powf(7.0))
+                    * dy;
+            }
+        }
+        for i in 0..boundary.points.len() {
+            let dx = self.points[index].x - boundary.points[i].x;
+            let dy = self.points[index].y - boundary.points[i].y;
+            let r = (dx.powf(2.0) + dy.powf(2.0)).powf(0.5);
+            f_x += 4.0
+                * epsilon
+                * (12.0 * sigma.powf(12.0) / r.powf(13.0) - 6.0 * sigma.powf(6.0) / r.powf(7.0))
+                * dx;
+            f_y += 4.0
+                * epsilon
+                * (12.0 * sigma.powf(12.0) / r.powf(13.0) - 6.0 * sigma.powf(6.0) / r.powf(7.0))
+                * dy;
+        }
+        Grid2D { x: f_x, y: f_y }
+    }
+
+    #[allow(dead_code)]
+    pub fn lennard_jones_potential_deriv(&mut self, index: usize, boundary: &Points2D) -> Grid2D {
+        let epsilon: f64 = 1.0;
+        let sigma: f64 = 0.9 * (2.0_f64).powf(-1.0 / 6.0) / (self.points.len() as f64).powf(0.5);
+        let mut f_x = 0.0;
+        let mut f_y = 0.0;
+        for i in 0..self.points.len() {
+            if i != index {
+                let dx = self.points[index].x - self.points[i].x;
+                let dy = self.points[index].y - self.points[i].y;
+                let r = (dx.powf(2.0) + dy.powf(2.0)).powf(0.5);
+                f_x += 4.0
+                    * epsilon
+                    * (12.0 * sigma.powf(12.0) / r.powf(13.0)
+                        - 6.0 * sigma.powf(6.0) / r.powf(7.0))
+                    * dx;
+                f_y += 4.0
+                    * epsilon
+                    * (12.0 * sigma.powf(12.0) / r.powf(13.0)
+                        - 6.0 * sigma.powf(6.0) / r.powf(7.0))
+                    * dy;
+            }
+        }
+        for i in 0..boundary.points.len() {
+            let dx = self.points[index].x - boundary.points[i].x;
+            let dy = self.points[index].y - boundary.points[i].y;
+            let r = (dx.powf(2.0) + dy.powf(2.0)).powf(0.5);
+            f_x += 4.0
+                * epsilon
+                * (12.0 * sigma.powf(12.0) / r.powf(13.0) - 6.0 * sigma.powf(6.0) / r.powf(7.0))
+                * dx;
+            f_y += 4.0
+                * epsilon
+                * (12.0 * sigma.powf(12.0) / r.powf(13.0) - 6.0 * sigma.powf(6.0) / r.powf(7.0))
+                * dy;
+        }
+        Grid2D { x: f_x, y: f_y }
+    }
+}