與 Rust 勾心鬥角 · 包圍球

garfileo發表於2022-05-12
Rust

OFF 檔案記錄的是多面體資訊，將其內容解析為 Mesh 結構後，便可基於後者為多面體網格構造包圍球。

球體

用泛型的結構體定義球體：

struct Sphere<T> {
    n: usize,       // 維度
    center: Vec<T>, // 中心
    radius: T       // 半徑
}

然後為該結構定義 new 方法：

impl<T> Sphere<T> {
    fn new(n: usize) -> Sphere<T> {
        Sphere{n: n, center: vec![0.0; n], radius: 0.0}
    }
}

倘若編譯該方法，rustc 會有以下指責：

  • vec![...] 的第一個引數的型別本該是 T，不是浮點型；
  • 為 Sphere 的 radius 成員賦的值，其型別應該是 T，不是浮點型；
  • vec![...] 的第一個引數需要實現 std::Clone Trait。

前兩個指責，是希望我們為 T 定義 0，因為 rustc 不知道 T 型別的 0 值的形式。第三個指責是希望為 T 增加約束。要解決這些問題，我能想出的方案是

use std::marker;

struct Sphere<T> {
    n: usize,
    center: Vec<T>,
    radius: T
}

trait Zero {
    fn zero() -> Self;
}

impl Zero for f64 {
    fn zero() -> Self {
        0.0
    }
}

impl<T: Zero + marker::Copy> Sphere<T> {
    fn new(n: usize) -> Sphere<T> {
        Sphere{n: n, center: vec![T::zero(); n], radius: T::zero()}
    }
}

我沒按 rustc 的建議使用 std::clone::Clone Trait 作為 T 的約束，而是使用了 std::marker::Copy Trait。

基於以上程式碼定義的球體，能夠支援以下語句：

let sphere: Sphere<f64> = Sphere::new(3);

趁熱打鐵，再為球體實現 Display Trait 吧，現在已經駕輕就熟了，

use std::fmt;

impl<T: fmt::Display> fmt::Display for Sphere<T> {
    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
        let mut info = String::new();
        info += "球體：";
        info += format!("維度 {}；", self.n).as_str();
        info += format!("中心 (").as_str();
        for i in 0 .. self.n - 1 {
            info += format!("{}, ", self.center[i]).as_str();
        }
        info += format!("{})；", self.center[self.n - 1]).as_str();
        info += format!("半徑 {}.", self.radius).as_str();
        write!(f, "{}", info)
    }
}

網格的中心

Mesh 例項的中心即包圍球的中心。現在為 Mesh 結構增加 center 方法，用於計算 Mesh 例項的中心，以下是該過程的基本泛型框架：

impl<T: Zero + marker::Copy> Mesh<T> {
    fn center(&self) -> Vec<T> {
        let mut x = vec![T::zero(); self.n];
        // 計算 self 的中心，將結果存於 x
        return x;
    }
}

我敢肯定，rustc 會根據具體的網格中心計算程式碼繼續要求我為 T 增加型別約束，而且這個過程也會讓我有些焦慮。倘若我毫不焦慮，而且對 rustc 有所感激，認為它飽含聖光，指出了我的程式碼的疏漏，那我敢肯定，我被 rustc PUA 了。

網格的中心，可以取為網格頂點集合的均值點：

// 計算 self 的中心，將結果存於 x
for x_i in &self.points {
    for j in 0 .. self.n {
        x[j] += x_i[j] / self.points.len() as T;
    }
}

對於上述程式碼，rustc 認為：

  • 它不知該如何進行型別 T 的除法運算；
  • self.n 為 usize 型別，它無法使用 as 轉換為 T 型別，因為 as 只能用於基本型別的轉換。

對於第一個問題，為 T 增加 std::ops::Div<Output = T> 約束便可解決。對於第二個問題，一種可行的方案是，為 T 增加 std::convert::From<usize> 約束，然後將 self.points.len() as T 修改為 self.points.len().into()，以實現 self.points.len() 的型別從 usize 到 T 的轉換。於是，Mesh 的 center 方法的程式碼變為

impl<T: Zero
     + marker::Copy
     + std::ops::Div<Output = T>
     + std::convert::From<usize>> Mesh<T> {
    fn center(&self) -> Vec<T> {
        let mut x = vec![T::zero(); self.n];
        // 計算 self 的中心，將結果存於 x
        for x_i in &self.points {
            for j in 0 .. self.n {
                x[j] += x_i[j] / self.points.len().into();
            }
        }
        return x;
    }
}

然而，rustc 意猶未盡，繼續認為它不知道該怎麼用 += 處理 T 型別的值，於是我需要繼續為 T 增加約束 std::ops::AddAssign，結果 Mesh 的 center 方法的程式碼變成

impl<T: Zero
     + marker::Copy
     + std::ops::Div<Output = T>
     + std::convert::From<usize>
     + std::ops::AddAssign> Mesh<T> {
    fn center(&self) -> Vec<T> {
        let mut x = vec![T::zero(); self.n];
        // 計算 self 的中心，將結果存於 x
        for x_i in &self.points {
            for j in 0 .. self.n {
                x[j] += x_i[j] / self.points.len().into();
            }
        }
        return x;
    }
}

這樣便萬事大吉了嗎？當然不是，rustc 會繼續認為 x[j] += x_i[j] / ... 裡的 x_i[j] 無法移動，原因是它對應的型別 T 未實現 copy Trait，因此不得不繼續為 T 追加 std::marker::Copy 約束。現在，Mesh 的 center 方法的程式碼變為

impl<T: Zero
     + marker::Copy
     + std::ops::Div<Output = T>
     + std::convert::From<usize>
     + std::ops::AddAssign
     + std::marker::Copy> Mesh<T> {
    fn center(&self) -> Vec<T> {
        let mut x = vec![T::zero(); self.n];
        // 計算 self 的中心，將結果存於 x
        for x_i in &self.points {
            for j in 0 .. self.n {
                x[j] += x_i[j] / self.points.len().into();
            }
        }
        return x;
    }
}

然後，rustc 不再說什麼，這時我才有餘力看出程式碼裡存在一處效能問題需要解決，即

for x_i in &self.points {
    for j in 0 .. self.n {
        x[j] += x_i[j] / self.points.len().into();
    }
}

需要修改為

let n: T = self.points.len().into();
for x_i in &self.points {
    for j in 0 .. self.n {
        x[j] += x_i[j] / n;
    }
}

以下程式碼可用於測試 Mesh 的 center 方法是否真的能算出多面體的中心：

let dim = 3;
let mut mesh: Mesh<f64> = Mesh::new(dim);
mesh.load("foo.off");

let center: Vec<f64> = mesh.center();
let mut sphere: Sphere<f64> = Sphere::new(dim);
for i in 0 .. dim {
    sphere.center[i] = center[i];
}
println!("{}", sphere);

但是，rustc 編譯上述程式碼時，會很傲驕地說 f64: From<usize> 沒實現，也就是說 Rust 標準庫裡為 f64 型別實現了一大堆的的 From<...>，然而唯獨沒實現 From<usize，亦即 Mesh 的 center 方法裡的程式碼

let n: T = self.points.len().into();

無法通過編譯。於是，之前的一堆努力，崩潰於這最後一片無辜的雪花。為了挽回敗局，我只好在程式碼裡用了武當派的梯雲縱，左腳踩右腳，右腳踩左腳，扶搖直上……

impl<T: Zero
     + marker::Copy
     + std::ops::Div<Output = T>
     + std::convert::From<f64>
     + std::ops::AddAssign
     + std::marker::Copy> Mesh<T> {
    fn center(&self) -> Vec<T> {
        let mut x = vec![T::zero(); self.n];
        // 計算 self 的中心，將結果存於 x
        let n: T = (self.points.len() as f64).into();
        for x_i in &self.points {
            for j in 0 .. self.n {
                x[j] += x_i[j] / n;
            }
        }
        return x;
    }
}

網格的半徑

網格的半徑是網格頂點到網格中心的最大距離，為便於實現該過程，先定義一個泛型函式，用於計算兩點間的距離：

fn distance<T>(a: &Vec<T>, b: &Vec<T>) -> T {
    let na = a.len();
    let nb = b.len();
    assert_eq!(na, nb);
    let mut d: T = T::zero();
    for i in 0 .. na {
        let t = a[i] - b[i];
        d += t * t;
    }
    return d.sqrt();
}

經過 rustc 的一番調教，distance 函式變為

fn distance<T: Zero
            + std::ops::Sub<Output = T>
            + std::ops::Mul<Output = T>
            + std::ops::AddAssign
            + Copy>(a: &Vec<T>, b: &Vec<T>) -> T {
    let na = a.len();
    let nb = b.len();
    assert_eq!(na, nb);
    let mut d: T = T::zero();
    for i in 0 .. na {
        let t = a[i] - b[i];
        d += t * t;
    }
    return d.sqrt();
}

即便如此，該函式依然無法通過編譯，因為 rustc 認為它無法確定 T 型別的例項有 sqrt 方法。既然天不佑我，那就別怪我程式碼寫得醜：

trait Sqrt<T> {
    fn sqrt(self) -> T;
}

fn distance<T: Zero
            + std::ops::Sub<Output = T>
            + std::ops::Mul<Output = T>
            + std::ops::AddAssign
            + Copy
            + Sqrt<T>>(a: &Vec<T>, b: &Vec<T>) -> T {
    let na = a.len();
    let nb = b.len();
    assert_eq!(na, nb);
    let mut d: T = T::zero();
    for i in 0 .. na {
        let t = a[i] - b[i];
        d += t * t;
    }
    return d.sqrt();
}

若點的座標值是 f64 型別，只需為該型別實現 Sqrt Trait，

impl Sqrt<f64> for f64 {
    fn sqrt(self) -> f64 {
        self.sqrt()
    }
}

便可使用 distance 計算兩點距離，例如

let a: Vec<f64> = vec![0.0, 0.0, 0.0];
let b: Vec<f64> = vec![1.0, 1.0, 1.0];
println!("{}", distance(&a, &b));

結果為 1.7320508075688772。

有了 distance 函式，便可計算網格半徑：

impl <T: Zero
      + std::ops::AddAssign
      + std::marker::Copy
      + Sqrt<T>
      + std::ops::Sub<Output = T>
      + std::ops::Mul<Output = T>
      + std::cmp::PartialOrd> Mesh<T> {
    fn radius(&self, center: &Vec<T>) -> T {
        let mut r = T::zero();
        for x in &self.points {
            let d = distance(x, center);
            if r < d {
                r = d;
            }
        }
        return r;
    }
}

要寫出上述程式碼，自然少不了 rustc 對型別的 T 各種具體約束的循循善誘……

網格的包圍球

現在，將 Mesh 的 center 和 radius 方法合併為 bounding_sphere：

impl<T: Zero
     + marker::Copy
     + std::ops::Div<Output = T>
     + std::convert::From<f64>
     + std::ops::AddAssign
     + Sqrt<T>
     + std::ops::Sub<Output = T>
     + std::ops::Mul<Output = T>
     + std::cmp::PartialOrd> Mesh<T> {
    fn bounding_sphere(&self) -> Sphere<T> {
        let mut sphere: Sphere<T> = Sphere::new(self.n);
        // 計算包圍球中心
        let n: T = (self.points.len() as f64).into();
        for x_i in &self.points {
            for j in 0 .. self.n {
                sphere.center[j] += x_i[j] / n;
            }
        }
        // 計算包圍球半徑
        for x in &self.points {
            let d = distance(x, &sphere.center);
            if sphere.radius < d {
                sphere.radius = d;
            }
        }
        return sphere;
    }
}

以下為 Mesh 的 bounding_sphere 方法的呼叫示例：

let dim = 3;
let mut mesh: Mesh<f64> = Mesh::new(dim);
mesh.load("foo.off");
let sphere: Sphere<f64> = mesh.bounding_sphere();
println!("{}", sphere);

Rust 泛型之我見

最好別用泛型。

最好別用泛型。

最好別用泛型。

小結

use std::{fmt, ops, convert, marker, cmp};
use std::path::Path;
use std::fs::File;
use std::io::{BufRead, BufReader};
use std::str::FromStr;
use std::num::ParseFloatError;
use std::ops::Index;

trait Zero {
    fn zero() -> Self;
}

impl Zero for f64 {
    fn zero() -> Self {
        0.0
    }
}

trait Length {
    fn len(&self) -> usize;
}

impl<T> Length for Vec<T> {
    fn len(&self) -> usize {
        return self.len();
    }
}

struct Mesh<T> {
    n: usize, // 維度
    points: Vec<Vec<T>>,  // 點表
    facets: Vec<Vec<usize>> // 面表
}

impl<T: FromStr<Err = ParseFloatError>> Mesh<T> {
    fn new(n: usize) -> Mesh<T> {
        return Mesh {n: n, points: Vec::new(), facets: Vec::new()};
    }
    fn load(&mut self, path: &str) {
        let path = Path::new(path);
        let file = File::open(path).unwrap();
        let buf = BufReader::new(file);

        let mut lines_iter = buf.lines().map(|l| l.unwrap());
        assert_eq!(lines_iter.next(), Some(String::from("OFF")));
        let second_line = lines_iter.next().unwrap();
        let mut split = second_line.split_whitespace();
        let n_of_points: usize = split.next().unwrap().parse().unwrap();
        let n_of_facets: usize = split.next().unwrap().parse().unwrap();

        for _i in 0 .. n_of_points {
            let line = lines_iter.next().unwrap();
            let mut p: Vec<T> = Vec::new();
            for x in line.split_whitespace() {
                p.push(x.parse().unwrap());
            }
            self.points.push(p);
        }
        for _i in 0 .. n_of_facets {
            let line = lines_iter.next().unwrap();
            let mut f: Vec<usize> = Vec::new();
            let mut split = line.split_whitespace();
            let n:usize = split.next().unwrap().parse().unwrap();
            assert_eq!(n, self.n);
            for x in split {
                f.push(x.parse().unwrap());
            }
            assert_eq!(n, f.len());
            self.facets.push(f);        
        }
    }
}

struct Prefix<T> {
    status: bool,
    body: fn(&T) -> String
}

impl<T> Prefix<T> {
    fn new() -> Prefix<T> {
        Prefix{status: false, body: |_| "".to_string()}
    }
}

fn matrix_fmt<T: Length + Index<usize>>(v: &Vec<T>,
                                        prefix: Prefix<T>) -> String
where <T as Index<usize>>::Output: fmt::Display,
      <T as Index<usize>>::Output: Sized {
    let mut s = String::new();
    for x in v {
        let n = x.len();
        if prefix.status {
            s += (prefix.body)(x).as_str();
        }
        for i in 0 .. n {
            if i == n - 1 {
                s += format!("{}\n", x[i]).as_str();
            } else {
                s += format!("{} ", x[i]).as_str();
            }
        }
    }
    return s;
}

impl<T: fmt::Display> fmt::Display for Mesh<T> {
    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
        let mut info = String::new();
        info += format!("OFF\n").as_str();
        info += format!("{0} {1} 0\n", self.points.len(), self.facets.len()).as_str();
        info += matrix_fmt(&self.points,  Prefix::new()).as_str();
        info += matrix_fmt(&self.facets, Prefix{status: true,
                                                body: |x| format!("{} ", x.len())}).as_str();
        write!(f, "{}", info)
    }
}

trait Sqrt<T> {
    fn sqrt(self) -> T;
}

impl Sqrt<f64> for f64 {
    fn sqrt(self) -> f64 {
        self.sqrt()
    }
}

fn distance<T: Zero
            + ops::Sub<Output = T>
            + ops::Mul<Output = T>
            + ops::AddAssign
            + marker::Copy
            + Sqrt<T>>(a: &Vec<T>, b: &Vec<T>) -> T {
    let na = a.len();
    let nb = b.len();
    assert_eq!(na, nb);
    let mut d: T = T::zero();
    for i in 0 .. na {
        let t = a[i] - b[i];
        d += t * t;
    }
    return d.sqrt();
}

struct Sphere<T> {
    n: usize,
    center: Vec<T>,
    radius: T
}

impl<T: Zero + marker::Copy> Sphere<T> {
    fn new(n: usize) -> Sphere<T> {
        Sphere{n: n, center: vec![T::zero(); n], radius: T::zero()}
    }
}

impl<T: fmt::Display> fmt::Display for Sphere<T> {
    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
        let mut info = String::new();
        info += "球體：";
        info += format!("維度 {}；", self.n).as_str();
        info += format!("中心 (").as_str();
        for i in 0 .. self.n - 1 {
            info += format!("{}, ", self.center[i]).as_str();
        }
        info += format!("{})；", self.center[self.n - 1]).as_str();
        info += format!("半徑 {}.", self.radius).as_str();
        write!(f, "{}", info)
    }
}

impl<T: Zero
     + ops::Div<Output = T>
     + convert::From<f64>
     + ops::AddAssign
     + marker::Copy
     + Sqrt<T>
     + ops::Sub<Output = T>
     + ops::Mul<Output = T>
     + cmp::PartialOrd> Mesh<T> {
    fn bounding_sphere(&self) -> Sphere<T> {
        let mut sphere: Sphere<T> = Sphere::new(self.n);
        // 計算包圍球中心
        let n: T = (self.points.len() as f64).into();
        for x_i in &self.points {
            for j in 0 .. self.n {
                sphere.center[j] += x_i[j] / n;
            }
        }
        // 計算包圍球半徑
        for x in &self.points {
            let d = distance(x, &sphere.center);
            if sphere.radius < d {
                sphere.radius = d;
            }
        }
        return sphere;
    }
}

fn main() {
    let dim = 3;
    let mut mesh: Mesh<f64> = Mesh::new(dim);
    mesh.load("foo.off");
    let sphere: Sphere<f64> = mesh.bounding_sphere();
    println!("{}", sphere);
}

贈品

Sphere 的 uv_to_xyz 方法可將球面上的經緯度座標轉化為三維座標，該方法僅適用於三維球體，以下為其程式碼實現：

use std::f64::consts::PI;

trait Math<T> {
    fn sqrt(self) -> T;
    fn sin(self) -> T;
    fn cos(self) -> T;
}
impl Math<f64> for f64 {
    fn sqrt(self) -> f64 {
        self.sqrt()
    }
    fn sin(self) -> f64 {
        self.sin()
    }
    fn cos(self) -> f64 {
        self.cos()
    }
}

impl<T: Zero
     + marker::Copy
     + ops::Div<Output = T>
     + convert::From<f64>
     + ops::Mul<Output = T>
     + Math<T>
     + ops::AddAssign> Sphere<T> {
    // 僅適用於三維球體
    fn uv_to_xyz(&self, u: T, v: T) -> Vec<T> {
        assert_eq!(self.n, 3);
        let mut xyz = vec![T::zero(); self.n];
        let u_r = (u / 180.0_f64.into()) * PI.into();
        let v_r = (v / 180.0_f64.into()) * PI.into();
        xyz[0] = self.radius * v_r.cos() * u_r.sin();
        xyz[1] = self.radius * v_r.sin();
        xyz[2] = self.radius * v_r.cos() * u_r.cos();
        // 平移
        for i in 0 .. self.n {
            xyz[i] += self.center[i];
        }
        return xyz;
    }
}

以下為測試程式碼，計算球面上東經 90 度，北緯 45 度的點的三維座標：

let xyz: Vec<f64> = sphere.uv_to_xyz(90.0, 45.0);
println!("({0}, {1}, {2})", xyz[0], xyz[1], xyz[2]);

倘若是西經或南緯，角度值為負值。

相關文章