//! # startin
//!
//! A Delaunay triangulator where the input are 2.5D points, the DT is computed in 2D but the elevation of the vertices are kept.
//! This is used mostly for the modelling of terrains.
//! Constructing a 2D Delaunay triangulation is also possible.
//!
//! The construction algorithm used is an incremental insertion based on flips, and the data structure is a cheap implementation of the star-based structure defined in [Blandford et al. (2003)](https://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.9.6823), cheap because the link of each vertex is stored a simple array (`Vec`) and not in an optimised blob like they did.
//! It results in a pretty fast library (comparison will come at some point), but it uses more space than the optimised one.
//!
//! The deletion of a vertex is also possible. The algorithm implemented is a modification of the one of [Mostafavi, Gold, and Dakowicz (2003)](https://doi.org/10.1016/S0098-3004(03)00017-7). The ears are filled by flipping, so it's in theory more robust.
//! I have also extended the algorithm to allow the deletion of vertices on the boundary of the convex hull.
//! The algorithm is sub-optimal, but in practice the number of neighbours of a given vertex in a DT is only 6, so it doesn't really matter.
//!
//! Robust arithmetic for the geometric predicates are used ([Shewchuk's predicates](https://www.cs.cmu.edu/~quake/robust.html), well the [Rust port of the code (robust crate)](https://crates.io/crates/robust)), so startin is robust and shouldn't crash (touch wood).
//!
//! There are a few interpolation functions implemented: (1) nearest-neighbour, (2) linear in TIN, (3) Laplace, (4) natural neighbour (aka Sibson's interpolation), (5) IDW.
//!
//!
//! # Web-demo with WebAssembly
//!
//! Rust can be compiled to [WebAssembly](https://www.rust-lang.org/what/wasm), and you can see a demo of some of the possibilities of startin (all computations are done locally and it's fast!).
//!
//! [--> web-demo](https://hugoledoux.github.io/startin_wasm/)
//!
//!
//! # Python bindings
//!
//! If you prefer Python, I made bindings: [https://github.com/hugoledoux/startinpy/](https://github.com/hugoledoux/startinpy/)
//!
//! There are a few more functions (eg reading GeoTIFF/LAZ, exporting GeoJSON/CityJSON) and it works with Numpy.
//!
//!
//! # C interface
//!
//! A basic C interface is available in `src/c_interface.rs`, to compile it:
//!
//! ```bash
//! cargo build --features c_api
//! ```
//!
//! # Usage
//!
//! ```rust
//! extern crate startin;
//!
//! fn main() {
//! let mut pts: Vec<[f64; 3]> = Vec::new();
//! pts.push([20.0, 30.0, 2.0]);
//! pts.push([120.0, 33.0, 12.5]);
//! pts.push([124.0, 222.0, 7.65]);
//! pts.push([20.0, 133.0, 21.0]);
//! pts.push([60.0, 60.0, 33.0]);
//! let mut dt = startin::Triangulation::new();
//! dt.insert(&pts, startin::InsertionStrategy::AsIs);
//! println!("{}", dt);
//! //-- print all the vertices
//! for (i, each) in dt.all_vertices().iter().enumerate() {
//! // skip the first one, the infinite vertex
//! if i > 0 {
//! println!("#{}: ({:.3}, {:.3}, {:.3})", i, each[0], each[1], each[2]);
//! }
//! }
//! //-- insert a new vertex
//! let re = dt.insert_one_pt(22.2, 33.3, 4.4);
//! match re {
//! Ok(_v) => println!(
//! "Inserted new point, now the DT has {} vertices",
//! dt.number_of_vertices()
//! ),
//! Err(v) => println!("Duplicate of vertex #{}, not inserted", v),
//! }
//! //-- remove it
//! match dt.remove(6) {
//! Ok(num) => println!("Vertex deleted, now the DT has {} vertices", num),
//! Err(why) => println!("!!! Deletion error: {:?}", why),
//! }
//! //-- get the convex hull
//! let ch = dt.convex_hull();
//! println!("Convex hull: {:?}", ch);
//! //-- fetch triangle containing (x, y)
//! match dt.locate(50.0, 50.0) {
//! Ok(tr) => println!("The triangle is {}", tr),
//! Err(why) => println!("Error: {:?}", why),
//! }
//! //-- interpolate with Laplace interpolation at 2 locations
//! let locs = vec![[51.0, 22.0], [50.3, 19.9]];
//! let interpolant = startin::interpolation::Laplace {};
//! let zs = startin::interpolation::interpolate(&interpolant, &mut dt, &locs);
//! for z in &zs {
//! match z {
//! Ok(value) => println!("z = {}", value),
//! Err(why) => println!("Interplation impossible: {:?}", why),
//! }
//! }
//!
//! //-- save the triangulation in geojson for debug purposes
//! let _re = dt.write_obj("/home/elvis/tr.obj".to_string());
//! }
//! ```
pub mod geom;
pub mod interpolation;
#[cfg(feature = "c_api")]
mod c_interface;
use rand::prelude::thread_rng;
use rand::Rng;
use std::fmt;
use std::fs::File;
use std::io::Write;
/// Errors that arise while using startin
#[derive(Debug, PartialEq)]
pub enum StartinError {
EmptyTriangulation,
OutsideConvexHull,
SearchCircleEmpty,
TriangleNotPresent,
VertexInfinite,
VertexRemoved,
VertexUnknown,
}
/// Possibilities for the insertion (with `insert()`)
pub enum InsertionStrategy {
AsIs,
BBox,
// Sprinkle,
// ConBRIO,
}
/// A Triangle is a triplet of indices
#[derive(Debug, PartialEq, Clone)]
pub struct Triangle {
pub v: [usize; 3],
}
impl Triangle {
/// Checks whether a Triangle is "infinite",
/// ie if one its vertices is the infinite vertex
fn is_infinite(&self) -> bool {
if self.v[0] == 0 || self.v[1] == 0 || self.v[2] == 0 {
return true;
}
return false;
}
}
impl fmt::Display for Triangle {
fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
write!(f, "({}, {}, {})", self.v[0], self.v[1], self.v[2])
}
}
//----------------------
#[repr(C)]
#[derive(Debug, Clone)]
struct Link(Vec<usize>);
impl Link {
fn new() -> Link {
// Link(Vec::new())
Link(Vec::with_capacity(8))
}
fn len(&self) -> usize {
self.0.len()
}
fn is_empty(&self) -> bool {
if self.0.len() == 0 {
true
} else {
false
}
}
fn add(&mut self, v: usize) {
self.0.push(v);
}
fn insert_after_v(&mut self, v: usize, after: usize) {
let pos = self.0.iter().position(|&x| x == after).unwrap();
self.0.insert(pos + 1, v);
}
fn delete(&mut self, v: usize) {
let re = self.0.iter().position(|&x| x == v);
if re != None {
self.0.remove(re.unwrap());
}
}
fn replace(&mut self, v: usize, newv: usize) {
let re = self.0.iter().position(|&x| x == v);
if re != None {
self.0[re.unwrap()] = newv;
// self.0.remove(re.unwrap());
}
}
fn infinite_first(&mut self) {
let re = self.0.iter().position(|&x| x == 0);
if re != None {
let posinf = re.unwrap();
if posinf == 0 {
return;
}
let mut newstar: Vec<usize> = Vec::new();
for j in posinf..self.0.len() {
newstar.push(self.0[j]);
}
for j in 0..posinf {
newstar.push(self.0[j]);
}
// println!("newstar: {:?}", newstar);
self.0 = newstar;
}
}
fn clear(&mut self) {
self.0.clear();
}
fn contains_infinite_vertex(&self) -> bool {
let pos = self.0.iter().position(|&x| x == 0);
if pos == None {
return false;
} else {
return true;
}
}
fn next_index(&self, i: usize) -> usize {
if i == (self.0.len() - 1) {
0
} else {
i + 1
}
}
fn prev_index(&self, i: usize) -> usize {
if i == 0 {
self.0.len() - 1
} else {
i - 1
}
}
fn get_index(&self, v: usize) -> Option<usize> {
return self.0.iter().position(|&x| x == v);
}
fn get_next_vertex(&self, v: usize) -> Option<usize> {
let re = self.get_index(v);
if re.is_none() {
return None;
}
let pos = re.unwrap();
if pos == (self.0.len() - 1) {
return Some(self.0[0]);
} else {
return Some(self.0[pos + 1]);
}
}
fn get_prev_vertex(&self, v: usize) -> Option<usize> {
let re = self.get_index(v);
if re.is_none() {
return None;
}
let pos = re.unwrap();
if pos == 0 {
return Some(self.0[self.0.len() - 1]);
} else {
return Some(self.0[pos - 1]);
}
}
fn iter(&self) -> Iter {
Iter(Box::new(self.0.iter()))
}
}
//-- taken from https://stackoverflow.com/questions/40668074/am-i-incorrectly-implementing-intoiterator-for-a-reference-or-is-this-a-rust-bug
struct Iter<'a>(Box<dyn Iterator<Item = &'a usize> + 'a>);
impl<'a> Iterator for Iter<'a> {
type Item = &'a usize;
fn next(&mut self) -> Option<Self::Item> {
self.0.next()
}
}
impl std::ops::Index<usize> for Link {
type Output = usize;
fn index(&self, idx: usize) -> &usize {
&self.0[idx as usize]
}
}
impl fmt::Display for Link {
fn fmt(&self, fmt: &mut fmt::Formatter) -> fmt::Result {
// fmt.write_str("pt: {}\n", self.pt)?;
fmt.write_str(&format!("link: {:?}\n", self.0))?;
Ok(())
}
}
/// A triangulation is a collection of Stars, each Star has its (x,y,z)
/// and a Link (an array of adjacent vertices, ordered CCW)
#[repr(C)]
struct Star {
pt: [f64; 3],
link: Link,
}
impl Star {
fn new(x: f64, y: f64, z: f64) -> Star {
let l = Link::new();
Star {
pt: [x, y, z],
link: l,
}
}
fn is_deleted(&self) -> bool {
self.link.is_empty()
}
}
/// Represents a triangulation
#[repr(C)]
pub struct Triangulation {
stars: Vec<Star>,
snaptol: f64,
cur: usize,
is_init: bool,
jump_and_walk: bool,
robust_predicates: bool,
removed_indices: Vec<usize>,
}
impl Triangulation {
pub fn new() -> Triangulation {
// TODO: allocate a certain number?
// let mut l: Vec<Star> = Vec::with_capacity(100000);
let mut l: Vec<Star> = Vec::new();
l.push(Star::new(f64::INFINITY, f64::INFINITY, f64::INFINITY));
let es: Vec<usize> = Vec::new();
Triangulation {
stars: l,
snaptol: 0.001,
cur: 0,
is_init: false,
jump_and_walk: false,
robust_predicates: true,
removed_indices: es,
}
}
fn insert_one_pt_init_phase(&mut self, x: f64, y: f64, z: f64) -> Result<usize, usize> {
let p: [f64; 3] = [x, y, z];
for i in 1..self.stars.len() {
if geom::distance2d_squared(&self.stars[i].pt, &p) <= (self.snaptol * self.snaptol) {
return Err(i);
}
}
self.collect_garbage();
//-- add point to Triangulation and create its empty star
self.stars.push(Star::new(x, y, z));
//-- form the first triangles (finite + infinite)
let l = self.stars.len();
if l >= 4 {
let a = l - 3;
let b = l - 2;
let c = l - 1;
let re = geom::orient2d(
&self.stars[a].pt,
&self.stars[b].pt,
&self.stars[c].pt,
self.robust_predicates,
);
if re == 1 {
// println!("init: ({},{},{})", a, b, c);
self.stars[0].link.add(a);
self.stars[0].link.add(c);
self.stars[0].link.add(b);
self.stars[a].link.add(0);
self.stars[a].link.add(b);
self.stars[a].link.add(c);
self.stars[b].link.add(0);
self.stars[b].link.add(c);
self.stars[b].link.add(a);
self.stars[c].link.add(0);
self.stars[c].link.add(a);
self.stars[c].link.add(b);
self.is_init = true;
} else if re == -1 {
// println!("init: ({},{},{})", a, c, b);
self.stars[0].link.add(a);
self.stars[0].link.add(b);
self.stars[0].link.add(c);
self.stars[a].link.add(0);
self.stars[a].link.add(c);
self.stars[a].link.add(b);
self.stars[b].link.add(0);
self.stars[b].link.add(a);
self.stars[b].link.add(c);
self.stars[c].link.add(0);
self.stars[c].link.add(b);
self.stars[c].link.add(a);
self.is_init = true;
}
}
self.cur = l - 1;
if self.is_init == true {
//-- insert the previous vertices in the dt
for j in 1..(l - 3) {
let tr = self.walk(&self.stars[j].pt);
// println!("found tr: {}", tr);
self.flip13(j, &tr);
self.update_dt(j);
}
}
Ok(self.cur)
}
/// Set a snap tolerance when inserting new points: if the newly inserted
/// one is closer than `snap_tol` to another one, then it is not inserted.
/// Avoids having very close vertices (like at 0.00007mm)
/// Default is 0.001unit (thus 1mm for most datasets).
pub fn set_snap_tolerance(&mut self, snaptol: f64) -> f64 {
if snaptol > 0.0 {
self.snaptol = snaptol;
}
self.snaptol
}
pub fn get_snap_tolerance(&self) -> f64 {
self.snaptol
}
/// Activate/deactive the jump-and-walk strategy for [`Triangulation::locate()`].
/// (deactivated by default)
/// If deactivated, then the walk starts from the last inserted triangle.
pub fn set_jump_and_walk(&mut self, b: bool) {
self.jump_and_walk = b;
}
/// Is using robut predicates (with [crate robust](https://docs.rs/robust))?
/// (activated by default)
pub fn is_using_robust_predicates(&self) -> bool {
self.robust_predicates
}
/// Activate/deactivate [robust predictates](https://docs.rs/robust)
pub fn use_robust_predicates(&mut self, b: bool) {
self.robust_predicates = b;
}
/// Insert a [`Vec`] of [`array`] (`[f64; 3]`) values.
/// If [`InsertionStrategy::AsIs`] is used, then [`Triangulation::insert_one_pt()`] is called
/// for each point in the order given.
///
/// # Arguments
///
/// * `pts` - a [`Vec`] of `[f64; 3]`
/// * `strategy` - the [`InsertionStrategy`] to use for the insertion
pub fn insert(&mut self, pts: &Vec<[f64; 3]>, strategy: InsertionStrategy) {
match strategy {
InsertionStrategy::BBox => {
//-- find the bbox
let mut bbox = geom::bbox2d(&pts);
//-- "padding" of the bbox to avoid conflicts
bbox[0] = bbox[0] - 10.0;
bbox[1] = bbox[1] - 10.0;
bbox[2] = bbox[2] + 10.0;
bbox[3] = bbox[3] + 10.0 ;
self.insert_with_bbox(&pts, &bbox);
}
InsertionStrategy::AsIs => {
for each in pts {
let _re = self.insert_one_pt(each[0], each[1], each[2]);
}
}
// InsertionStrategy::Sprinkle => println!("Sprinkle not implemented yet"),
// InsertionStrategy::ConBRIO => println!("ConBRIO not implemented yet"),
}
}
fn insert_with_bbox(&mut self, pts: &Vec<[f64; 3]>, bbox: &[f64; 4]) {
let mut c4: Vec<usize> = Vec::new();
//-- insert the 4 corners
c4.push(self.insert_one_pt(bbox[0], bbox[1], 0.0).unwrap());
c4.push(self.insert_one_pt(bbox[2], bbox[1], 0.0).unwrap());
c4.push(self.insert_one_pt(bbox[2], bbox[3], 0.0).unwrap());
c4.push(self.insert_one_pt(bbox[0], bbox[3], 0.0).unwrap());
for each in pts {
let _re = self.insert_one_pt(each[0], each[1], each[2]);
}
//-- remove the 4 corners
for each in &c4 {
let _re = self.remove(*each);
}
//-- collect garbage: remove the 4 added vertices and "shift" all the vertex ids
self.collect_garbage();
}
/// Insert the point (`px`, `py`, `pz`) in the triangulation.
/// Returns the vertex ID of the point if the vertex didn't exist.
/// If there was a vertex at that location, an Error is thrown with the already
/// existing vertex ID.
pub fn insert_one_pt(&mut self, px: f64, py: f64, pz: f64) -> Result<usize, usize> {
if self.is_init == false {
return self.insert_one_pt_init_phase(px, py, pz);
}
//-- walk
let p: [f64; 3] = [px, py, pz];
let tr = self.walk(&p);
// println!("STARTING TR: {}", tr);
if geom::distance2d_squared(&self.stars[tr.v[0]].pt, &p) <= (self.snaptol * self.snaptol) {
return Err(tr.v[0]);
}
if geom::distance2d_squared(&self.stars[tr.v[1]].pt, &p) <= (self.snaptol * self.snaptol) {
return Err(tr.v[1]);
}
if geom::distance2d_squared(&self.stars[tr.v[2]].pt, &p) <= (self.snaptol * self.snaptol) {
return Err(tr.v[2]);
}
//-- ok we now insert the point in the data structure
let pi: usize;
if self.removed_indices.is_empty() == true {
self.stars.push(Star::new(px, py, pz));
pi = self.stars.len() - 1;
} else {
// self.stars.push(Star::new(px, py, pz));
pi = self.removed_indices.pop().unwrap();
self.stars[pi].pt[0] = px;
self.stars[pi].pt[1] = py;
self.stars[pi].pt[2] = pz;
}
//-- flip13()
self.flip13(pi, &tr);
//-- update_dt()
self.update_dt(pi);
self.cur = pi;
Ok(pi)
}
fn update_dt(&mut self, pi: usize) {
// println!("--> Update DT");
let mut mystack: Vec<Triangle> = Vec::new();
mystack.push(Triangle {
v: [pi, self.stars[pi].link[0], self.stars[pi].link[1]],
});
mystack.push(Triangle {
v: [pi, self.stars[pi].link[1], self.stars[pi].link[2]],
});
mystack.push(Triangle {
v: [pi, self.stars[pi].link[2], self.stars[pi].link[0]],
});
loop {
let tr = match mystack.pop() {
None => break,
Some(x) => x,
};
let opposite = self.get_opposite_vertex(&tr);
// println!("stacked: {} {}", tr, opposite);
if tr.is_infinite() == true {
let mut a: i8 = 0;
if tr.v[0] == 0 {
a = geom::orient2d(
&self.stars[opposite].pt,
&self.stars[tr.v[1]].pt,
&self.stars[tr.v[2]].pt,
self.robust_predicates,
);
} else if tr.v[1] == 0 {
a = geom::orient2d(
&self.stars[tr.v[0]].pt,
&self.stars[opposite].pt,
&self.stars[tr.v[2]].pt,
self.robust_predicates,
);
} else if tr.v[2] == 0 {
a = geom::orient2d(
&self.stars[tr.v[0]].pt,
&self.stars[tr.v[1]].pt,
&self.stars[opposite].pt,
self.robust_predicates,
);
}
// println!("TODO: INCIRCLE FOR INFINITY {}", a);
if a > 0 {
// println!("FLIPPED0 {} {}", tr, opposite);
let (ret0, ret1) = self.flip22(&tr, opposite);
mystack.push(ret0);
mystack.push(ret1);
}
} else {
if opposite == 0 {
//- if insertion on CH then break the edge, otherwise do nothing
//-- TODO sure the flips are okay here?
if geom::orient2d(
&self.stars[tr.v[0]].pt,
&self.stars[tr.v[1]].pt,
&self.stars[tr.v[2]].pt,
self.robust_predicates,
) == 0
{
// println!("FLIPPED1 {} {}", tr, 0);
let (ret0, ret1) = self.flip22(&tr, 0);
mystack.push(ret0);
mystack.push(ret1);
}
} else {
if geom::incircle(
&self.stars[tr.v[0]].pt,
&self.stars[tr.v[1]].pt,
&self.stars[tr.v[2]].pt,
&self.stars[opposite].pt,
self.robust_predicates,
) > 0
{
// println!("FLIPPED2 {} {}", tr, opposite);
let (ret0, ret1) = self.flip22(&tr, opposite);
mystack.push(ret0);
mystack.push(ret1);
}
}
}
}
}
fn flip13(&mut self, pi: usize, tr: &Triangle) {
self.stars[pi].link.add(tr.v[0]);
self.stars[pi].link.add(tr.v[1]);
self.stars[pi].link.add(tr.v[2]);
self.stars[tr.v[0]].link.insert_after_v(pi, tr.v[1]);
self.stars[tr.v[1]].link.insert_after_v(pi, tr.v[2]);
self.stars[tr.v[2]].link.insert_after_v(pi, tr.v[0]);
//-- put infinite vertex first in list
// self.stars[pi].link.infinite_first();
}
fn flip31(&mut self, v: usize) {
// println!("FLIP31");
let mut ns: Vec<usize> = Vec::new();
for each in self.stars[v].link.iter() {
ns.push(*each);
}
for n in ns.iter() {
self.stars[*n].link.delete(v);
}
self.stars[v].link.clear();
self.stars[v].pt[0] = f64::NAN;
self.stars[v].pt[1] = f64::NAN;
self.stars[v].pt[2] = f64::NAN;
self.removed_indices.push(v);
if ns[0] != 0 {
self.cur = ns[0];
} else if ns[1] != 0 {
self.cur = ns[1];
} else if ns[2] != 0 {
self.cur = ns[2];
}
}
/// Returns the coordinates of the vertex v.
/// A [`StartinError`] is returned if `vi` doesn't exist
/// or is a removed vertex.
pub fn get_point(&self, vi: usize) -> Result<Vec<f64>, StartinError> {
match self.is_vertex_removed(vi) {
Err(why) => return Err(why),
Ok(b) => match b {
true => return Err(StartinError::VertexRemoved),
false => Ok(self.stars[vi].pt.to_vec()),
},
}
}
/// Returns the 3 adjacents (finite + infinite) [`Triangle`] to a triangle.
pub fn adjacent_triangles_to_triangle(
&self,
tr: &Triangle,
) -> Result<Vec<Triangle>, StartinError> {
if self.is_triangle(&tr) == false || tr.is_infinite() == true {
return Err(StartinError::TriangleNotPresent);
}
let mut trs: Vec<Triangle> = Vec::new();
let mut opp = self.stars[tr.v[2]].link.get_next_vertex(tr.v[1]).unwrap();
trs.push(Triangle {
v: [tr.v[1], opp, tr.v[2]],
});
opp = self.stars[tr.v[0]].link.get_next_vertex(tr.v[2]).unwrap();
trs.push(Triangle {
v: [tr.v[2], opp, tr.v[0]],
});
opp = self.stars[tr.v[1]].link.get_next_vertex(tr.v[0]).unwrap();
trs.push(Triangle {
v: [tr.v[0], opp, tr.v[1]],
});
Ok(trs)
}
/// Returns a [`Vec`] of [`Triangle`]s (finite + infinite) to the vertex `vi`.
pub fn incident_triangles_to_vertex(&self, vi: usize) -> Result<Vec<Triangle>, StartinError> {
match self.is_vertex_removed(vi) {
Err(why) => return Err(why),
Ok(b) => match b {
true => return Err(StartinError::VertexRemoved),
false => {
let mut trs: Vec<Triangle> = Vec::new();
for (i, each) in self.stars[vi].link.iter().enumerate() {
let j = self.stars[vi].link.next_index(i);
trs.push(Triangle {
v: [vi, *each, self.stars[vi].link[j]],
});
}
Ok(trs)
}
},
}
}
/// Returns the degree of the vertex with ID `vi`.
pub fn degree(&self, vi: usize) -> Result<usize, StartinError> {
match self.is_vertex_removed(vi) {
Err(why) => return Err(why),
Ok(b) => match b {
true => return Err(StartinError::VertexRemoved),
false => return Ok(self.stars[vi].link.len()),
},
}
}
/// Returns a list (`Vec<usize>`) (ordered CCW) of the adjacent vertices to `vi`.
pub fn adjacent_vertices_to_vertex(&self, vi: usize) -> Result<Vec<usize>, StartinError> {
match self.is_vertex_removed(vi) {
Err(why) => return Err(why),
Ok(b) => match b {
true => return Err(StartinError::VertexRemoved),
false => {
let mut adjs: Vec<usize> = Vec::new();
for each in self.stars[vi].link.iter() {
adjs.push(*each);
}
return Ok(adjs);
}
},
}
}
/// Returns whether a triplet of indices is a [`Triangle`] in the triangulation.
pub fn is_triangle(&self, tr: &Triangle) -> bool {
if tr.v[0] >= self.stars.len() || tr.v[1] >= self.stars.len() || tr.v[2] >= self.stars.len()
{
return false;
}
let re = self.stars[tr.v[0]].link.get_next_vertex(tr.v[1]);
if re.is_none() {
return false;
} else {
if re.unwrap() == tr.v[2] {
return true;
} else {
return false;
}
}
}
/// Returns whether a [`Triangle`] is finite, or not
pub fn is_finite(&self, tr: &Triangle) -> bool {
if tr.is_infinite() {
return false;
} else {
return true;
}
}
/// Returns some statistics about the triangulation.
pub fn statistics_degree(&self) -> (f64, usize, usize) {
let mut total: f64 = 0.0;
let mut min: usize = usize::max_value();
let mut max: usize = usize::min_value();
for i in 1..self.stars.len() {
total = total + self.stars[i].link.len() as f64;
if self.stars[i].link.len() > max {
max = self.stars[i].link.len();
}
if self.stars[i].link.len() < min {
min = self.stars[i].link.len();
}
}
total = total / (self.stars.len() - 2) as f64;
return (total, min, max);
}
/// Returns number of finite vertices in the triangulation.
/// The removed vertices are not counted.
pub fn number_of_vertices(&self) -> usize {
//-- number of finite vertices
self.stars.len() - 1 - self.removed_indices.len()
}
/// Returns number of finite triangles in the triangulation.
pub fn number_of_triangles(&self) -> usize {
//-- number of finite triangles
let mut count: usize = 0;
for (i, star) in self.stars.iter().enumerate() {
for (j, value) in star.link.iter().enumerate() {
if i < *value {
let k = star.link[star.link.next_index(j)];
if i < k {
let tr = Triangle { v: [i, *value, k] };
if tr.is_infinite() == false {
count = count + 1;
}
}
}
}
}
count
}
/// Returns the number of vertices which are marked as "removed"
pub fn number_of_removed_vertices(&self) -> usize {
self.removed_indices.len()
}
/// Returns whether the vertex `vi` is removed or not.
pub fn is_vertex_removed(&self, vi: usize) -> Result<bool, StartinError> {
if vi >= self.stars.len() {
return Err(StartinError::VertexUnknown);
}
Ok(self.stars[vi].is_deleted())
}
/// Returns the convex hull of the dataset, oriented CCW.
/// It is a list of vertex indices (first != last)
pub fn convex_hull(&self) -> Vec<usize> {
let mut re: Vec<usize> = Vec::new();
for x in self.stars[0].link.iter() {
re.push(*x);
}
re.reverse();
re
}
/// Returns the size (ie the number of vertices) of the convex hull of the dataset
pub fn number_of_vertices_on_convex_hull(&self) -> usize {
//-- number of finite vertices on the boundary of the convex hull
if self.is_init == false {
return 0;
}
return self.stars[0].link.len();
}
/// Returns `true` if the vertex `vi` is part of the boundary of the convex
/// hull of the dataset; `false` otherwise.
pub fn is_vertex_convex_hull(&self, vi: usize) -> bool {
if vi == 0 {
return false;
}
if self.is_vertex_valid(vi) == false {
return false;
}
self.stars[vi].link.contains_infinite_vertex()
}
/// Returns, if it exists, the [`Triangle`] containing `(px, py)`.
/// If it is direction on a vertex/edge, then one is randomly chosen.
pub fn locate(&self, px: f64, py: f64) -> Result<Triangle, StartinError> {
if self.is_init == false {
return Err(StartinError::EmptyTriangulation);
}
let p: [f64; 3] = [px, py, 0.0];
let re = self.walk(&p);
match re.is_infinite() {
true => return Err(StartinError::OutsideConvexHull),
false => return Ok(re),
}
}
/// Returns closest point (in 2D) to a query point `(x, y)`.
/// if `(px, py)` is outside the convex hull then [`StartinError::OutsideConvexHull`] is raised.
pub fn closest_point(&self, px: f64, py: f64) -> Result<usize, StartinError> {
let re = self.locate(px, py);
if re.is_err() {
return Err(re.err().unwrap());
}
let p: [f64; 3] = [px, py, 0.0];
let tr = re.unwrap();
let mut d = std::f64::MAX;
let mut closest: usize = 0;
//-- find closest vertex in the triangle containing p
for each in tr.v.iter() {
let dtmp = geom::distance2d_squared(&self.stars[*each].pt, &p);
if dtmp < d {
d = dtmp;
closest = *each;
}
}
loop {
let mut found_one_closer = false;
for each in self.stars[closest].link.iter() {
let dtmp = geom::distance2d_squared(&self.stars[*each].pt, &p);
if dtmp < d {
d = dtmp;
closest = *each;
found_one_closer = true;
break;
}
}
if found_one_closer == false {
break;
}
}
Ok(closest)
}
fn walk(&self, x: &[f64]) -> Triangle {
//-- find the starting tr
let mut cur = self.cur;
//-- jump-and-walk
if self.jump_and_walk == true {
let mut rng = thread_rng();
let mut d: f64 = geom::distance2d_squared(&self.stars[self.cur].pt, &x);
let n = (self.stars.len() as f64).powf(0.25);
// let n = (self.stars.len() as f64).powf(0.25) * 7.0;
for _i in 0..n as i32 {
let re: usize = rng.gen_range(1, self.stars.len());
// let dtemp = x.square_2d_distance(&self.stars[re].pt);
if self.stars[re].is_deleted() == true {
continue;
}
let dtemp = geom::distance2d_squared(&self.stars[re].pt, &x);
if dtemp < d {
cur = re;
d = dtemp;
}
}
}
let mut tr = Triangle { v: [0, 0, 0] };
// println!("cur: {}", cur);
//-- 1. find a finite triangle
tr.v[0] = cur;
let l = &self.stars[cur].link;
for i in 0..(l.len() - 1) {
if (l[i] != 0) && (l[i + 1] != 0) {
tr.v[1] = l[i];
tr.v[2] = l[i + 1];
break;
}
}
//-- 2. order it such that tr0-tr1-x is CCW
if geom::orient2d(
&self.stars[tr.v[0]].pt,
&self.stars[tr.v[1]].pt,
&x,
self.robust_predicates,
) == -1
{
if geom::orient2d(
&self.stars[tr.v[1]].pt,
&self.stars[tr.v[2]].pt,
&x,
self.robust_predicates,
) != -1
{
let tmp: usize = tr.v[0];
tr.v[0] = tr.v[1];
tr.v[1] = tr.v[2];
tr.v[2] = tmp;
} else {
let tmp: usize = tr.v[1];
tr.v[1] = tr.v[0];
tr.v[0] = tr.v[2];
tr.v[2] = tmp;
}
}
//-- 3. start the walk
//-- we know that tr0-tr1-x is CCW
loop {
if tr.is_infinite() == true {
break;
}
if geom::orient2d(
&self.stars[tr.v[1]].pt,
&self.stars[tr.v[2]].pt,
&x,
self.robust_predicates,
) != -1
{
if geom::orient2d(
&self.stars[tr.v[2]].pt,
&self.stars[tr.v[0]].pt,
&x,
self.robust_predicates,
) != -1
{
break;
} else {
//-- walk to incident to tr1,tr2
// println!("here");
let prev = self.stars[tr.v[2]].link.get_prev_vertex(tr.v[0]).unwrap();
tr.v[1] = tr.v[2];
tr.v[2] = prev;
}
} else {
//-- walk to incident to tr1,tr2
// a.iter().position(|&x| x == 2), Some(1)
let prev = self.stars[tr.v[1]].link.get_prev_vertex(tr.v[2]).unwrap();
tr.v[0] = tr.v[2];
tr.v[2] = prev;
}
}
return tr;
}
fn flip22(&mut self, tr: &Triangle, opposite: usize) -> (Triangle, Triangle) {
//-- step 1.
self.stars[tr.v[0]].link.insert_after_v(opposite, tr.v[1]);
//-- step 2.
self.stars[tr.v[1]].link.delete(tr.v[2]);
//-- step 3.
self.stars[opposite].link.insert_after_v(tr.v[0], tr.v[2]);
//-- step 4.
self.stars[tr.v[2]].link.delete(tr.v[1]);
//-- make 2 triangles to return (to stack)
let ret0 = Triangle {
v: [tr.v[0], tr.v[1], opposite],
};
let ret1 = Triangle {
v: [tr.v[0], opposite, tr.v[2]],
};
(ret0, ret1)
}
fn get_opposite_vertex(&self, tr: &Triangle) -> usize {
self.stars[tr.v[2]].link.get_next_vertex(tr.v[1]).unwrap()
}
/// Returns a [`Vec`]<[`Vec`]<[`f64`]>> of all the vertices
/// (including the infinite one and the removed ones)
pub fn all_vertices(&self) -> Vec<Vec<f64>> {
let mut pts: Vec<Vec<f64>> = Vec::with_capacity(self.stars.len() - 1);
for i in 0..self.stars.len() {
pts.push(self.stars[i].pt.to_vec());
}
pts
}
/// Returns a [`Vec`]<[`usize`]> of all the finite edges (implicitly grouped by 2)
pub fn all_finite_edges(&self) -> Vec<usize> {
let mut edges: Vec<usize> = Vec::new();
for i in 1..self.stars.len() {
for value in self.stars[i].link.iter() {
if (*value != 0) && (i < *value) {
edges.push(i);
edges.push(*value);
}
}
}
edges
}
/// Returns a [`Vec`]<[`Triangle`]> of all the (finite + infinite) triangles
pub fn all_triangles(&self) -> Vec<Triangle> {
let mut trs: Vec<Triangle> = Vec::new();
for (i, star) in self.stars.iter().enumerate() {
//-- reconstruct triangles
for (j, value) in star.link.iter().enumerate() {
if i < *value {
// let k = star.l[self.nexti(star.link.len(), j)];
let k = star.link[star.link.next_index(j)];
if i < k {
trs.push(Triangle { v: [i, *value, k] });
}
}
}
}
trs
}
/// Returns a [`Vec`]<[`Triangle`]> of all the finite triangles
pub fn all_finite_triangles(&self) -> Vec<Triangle> {
let alltrs = self.all_triangles();
let mut re: Vec<Triangle> = Vec::new();
for t in &alltrs {
if t.is_infinite() == false {
re.push(t.clone());
}
}
re
}
/// Validates the Delaunay triangulation:
/// (1) checks each triangle against each vertex (circumcircle tests); very slow
/// (2) checks whether the convex hull is really convex
pub fn is_valid(&self) -> bool {
self.is_valid_ch_convex() && self.is_valid_circumcircle()
}
fn is_valid_circumcircle(&self) -> bool {
let mut re = true;
let trs = self.all_finite_triangles();
for tr in trs.iter() {
for i in 1..self.stars.len() {
if self.stars[i].is_deleted() == false
&& geom::incircle(
&self.stars[tr.v[0]].pt,
&self.stars[tr.v[1]].pt,
&self.stars[tr.v[2]].pt,
&self.stars[i].pt,
self.robust_predicates,
) > 0
{
println!("NOT DELAUNAY FFS!");
println!("{} with {}", tr, i);
re = false
}
}
}
re
}
fn is_valid_ch_convex(&self) -> bool {
let mut re = true;
let ch = self.convex_hull();
for i in 0..ch.len() {
if geom::orient2d(
&self.stars[ch[i % ch.len()]].pt,
&self.stars[ch[(i + 1) % ch.len()]].pt,
&self.stars[ch[(i + 2) % ch.len()]].pt,
self.robust_predicates,
) == -1
{
re = false;
break;
}
}
if re == false {
println!("CONVEX NOT CONVEX");
}
return re;
}
fn remove_on_convex_hull(&mut self, v: usize) -> Result<usize, StartinError> {
// println!("!!! REMOVE ON CONVEX HULL");
let mut adjs: Vec<usize> = Vec::new();
//-- necessary because assumptions below for start-end line on CH
self.stars[v].link.infinite_first();
for each in self.stars[v].link.iter() {
adjs.push(*each);
}
// println!("adjs: {:?}", adjs);
let mut cur: usize = 0;
//-- 1. find and create finite triangles only
let mut nadjs = adjs.len();
let mut steps = 0;
while adjs.len() > 3 {
//-- control the loops to avoid infinite loop, when all options in a temp
//-- star have been tried it's because we're stuck (and done actually)
if steps == nadjs {
break;
}
if adjs.len() == nadjs {
steps += 1;
} else {
nadjs = adjs.len();
steps = 0;
}
//-- define the ear
let a = cur % adjs.len();
let b = (cur + 1) % adjs.len();
let c = (cur + 2) % adjs.len();
// println!("cur ear--> {:?} {}/{}/{}", adjs, a, b, c);
if adjs[a] == 0 || adjs[b] == 0 || adjs[c] == 0 {
//-- do not process infinite ear
cur += 1;
continue;
}
if (geom::orient2d(
&self.stars[adjs[a]].pt,
&self.stars[adjs[b]].pt,
&self.stars[adjs[c]].pt,
self.robust_predicates,
) == 1)
&& (geom::orient2d(
&self.stars[adjs[a]].pt,
&self.stars[adjs[c]].pt,
&self.stars[v].pt,
self.robust_predicates,
) >= 0)
{
// println!("ear {}-{}-{}", adjs[a], adjs[b], adjs[c]);
//-- test incircle with all other vertices in the "hole"
let cur2 = cur + 3;
let mut isdel = true;
for i in 0..adjs.len() - 3 {
// println!("test ear with {}", adjs[(cur2 + i) % adjs.len()]);
if adjs[(cur2 + i) % adjs.len()] != 0
&& geom::incircle(
&self.stars[adjs[a]].pt,
&self.stars[adjs[b]].pt,
&self.stars[adjs[c]].pt,
&self.stars[adjs[(cur2 + i) % adjs.len()]].pt,
self.robust_predicates,
) > 0
{
isdel = false;
break;
}
}
if isdel == true {
// println!("flip22");
let t = Triangle {
v: [adjs[a], adjs[b], v],
};
self.flip22(&t, adjs[c]);
adjs.remove((cur + 1) % adjs.len());
}
}
cur += 1;
}
//-- flip31 to remove the vertex
if adjs.len() == 3 {
self.flip31(v);
return Ok(self.stars.len() - 1);
} else {
//-- convex part is filled, and we need to apply a special "flip"
//-- to delete the vertex v and its incident edges
// println!("FLIP-FOR-CH");
self.stars[adjs[1]].link.delete(v);
self.stars[*(adjs.last().unwrap())].link.delete(v);
for i in 2..(adjs.len() - 1) {
if self.is_vertex_convex_hull(adjs[i]) == true {
//-- going back to a line, no triangles
//-- wipe it all and start the insert_init_phase again
for i in 0..self.stars.len() {
self.stars[i].link.clear();
}
self.stars[v].pt[0] = f64::NAN;
self.stars[v].pt[1] = f64::NAN;
self.stars[v].pt[2] = f64::NAN;
self.removed_indices.push(v);
self.is_init = false;
return Ok(self.stars.len() - 1);
}
self.stars[adjs[i]].link.replace(v, 0);
self.stars[adjs[i]].link.infinite_first();
}
let mut prev = v;
for i in 2..(adjs.len() - 1) {
self.stars[0].link.insert_after_v(adjs[i], prev);
prev = adjs[i];
}
self.stars[adjs[0]].link.delete(v);
self.stars[v].link.clear();
self.stars[v].pt[0] = f64::NAN;
self.stars[v].pt[1] = f64::NAN;
self.stars[v].pt[2] = f64::NAN;
self.removed_indices.push(v);
for i in 0..1000 {
if adjs[i] != 0 {
self.cur = adjs[i];
break;
}
}
Ok(self.stars.len() - 1)
}
}
/// Removes the vertex `vi` from the [`Triangulation`] and updates for the "Delaunay-ness".
///
/// The vertex is not removed from memory but flagged as removed, thus all the other vertices
/// keep their IDs.
/// The following insertion of a point will reuse this ID.
/// It is therefore possible to have an array that contains unused/removed vertices.
pub fn remove(&mut self, vi: usize) -> Result<usize, StartinError> {
// println!("REMOVE vertex {}", v);
if vi == 0 {
return Err(StartinError::VertexInfinite);
}
if self.is_init == false {
self.stars[vi].pt[0] = f64::NAN;
self.stars[vi].pt[1] = f64::NAN;
self.stars[vi].pt[2] = f64::NAN;
self.removed_indices.push(vi);
}
match self.is_vertex_removed(vi) {
Err(why) => return Err(why),
Ok(b) => {
if b == true {
return Err(StartinError::VertexRemoved);
}
}
}
if self.is_vertex_convex_hull(vi) {
return self.remove_on_convex_hull(vi);
}
let mut adjs: Vec<usize> = Vec::new();
for each in self.stars[vi].link.iter() {
adjs.push(*each);
}
// println!("adjs: {:?}", adjs);
let mut cur: usize = 0;
while adjs.len() > 3 {
let a = cur % adjs.len();
let b = (cur + 1) % adjs.len();
let c = (cur + 2) % adjs.len();
// println!("cur ear--> {:?} {}/{}/{}", adjs, a, b, c);
if (geom::orient2d(
&self.stars[adjs[a]].pt,
&self.stars[adjs[b]].pt,
&self.stars[adjs[c]].pt,
self.robust_predicates,
) == 1)
&& (geom::orient2d(
&self.stars[adjs[a]].pt,
&self.stars[adjs[c]].pt,
&self.stars[vi].pt,
self.robust_predicates,
) >= 0)
{
// println!("ear {}-{}-{}", adjs[a], adjs[b], adjs[c]);
//-- test incircle with all other vertices in the "hole"
let cur2 = cur + 3;
let mut isdel = true;
for i in 0..adjs.len() - 3 {
// println!("test ear with {}", adjs[(cur2 + i) % adjs.len()]);
if geom::incircle(
&self.stars[adjs[a]].pt,
&self.stars[adjs[b]].pt,
&self.stars[adjs[c]].pt,
&self.stars[adjs[(cur2 + i) % adjs.len()]].pt,
self.robust_predicates,
) > 0
{
isdel = false;
break;
}
}
if isdel == true {
// println!("flip22");
let t = Triangle {
v: [adjs[a], adjs[b], vi],
};
self.flip22(&t, adjs[c]);
adjs.remove((cur + 1) % adjs.len());
}
}
cur = cur + 1;
}
//-- flip31 to remove the vertex
self.flip31(vi);
Ok(self.stars.len() - 1)
}
/// Write an [OBJ file](https://en.wikipedia.org/wiki/Wavefront_.obj_file) to disk.
pub fn write_obj(&self, path: String) -> std::io::Result<()> {
let trs = self.all_finite_triangles();
let mut f = File::create(path)?;
let mut s = String::new();
//-- find one good vertice to replace the deleted one
let mut onegoodpt = vec![1.0, 1.0, 1.0];
for i in 1..self.stars.len() {
if self.stars[i].is_deleted() == false {
onegoodpt[0] = self.stars[i].pt[0];
onegoodpt[1] = self.stars[i].pt[1];
onegoodpt[2] = self.stars[i].pt[2];
break;
}
}
for i in 1..self.stars.len() {
if self.stars[i].is_deleted() == true {
s.push_str(&format!(
"v {} {} {}\n",
onegoodpt[0], onegoodpt[1], onegoodpt[2]
));
continue;
}
s.push_str(&format!(
"v {} {} {}\n",
self.stars[i].pt[0], self.stars[i].pt[1], self.stars[i].pt[2]
));
}
write!(f, "{}", s).unwrap();
let mut s = String::new();
for tr in trs.iter() {
s.push_str(&format!("f {} {} {}\n", tr.v[0], tr.v[1], tr.v[2]));
}
write!(f, "{}", s).unwrap();
// println!("write fobj: {:.2?}", starttime.elapsed());
Ok(())
}
/// Write a [PLY file](https://en.wikipedia.org/wiki/PLY_(file_format)) to disk.
pub fn write_ply(&self, path: String) -> std::io::Result<()> {
let trs = self.all_finite_triangles();
let mut f = File::create(path)?;
//-- header
write!(f, "ply\n").unwrap();
write!(f, "format ascii 1.0\n").unwrap();
write!(f, "comment made by startin\n").unwrap();
write!(f, "element vertex {}\n", self.stars.len() - 1).unwrap();
write!(f, "property float x\n").unwrap();
write!(f, "property float y\n").unwrap();
write!(f, "property float z\n").unwrap();
write!(f, "element face {}\n", trs.len()).unwrap();
write!(f, "property list uchar int vertex_indices\n").unwrap();
write!(f, "end_header\n").unwrap();
//-- find one good vertice to replace the deleted one
let mut onegoodpt = vec![1.0, 1.0, 1.0];
for i in 1..self.stars.len() {
if self.stars[i].is_deleted() == false {
onegoodpt[0] = self.stars[i].pt[0];
onegoodpt[1] = self.stars[i].pt[1];
onegoodpt[2] = self.stars[i].pt[2];
break;
}
}
let mut s = String::new();
for i in 1..self.stars.len() {
if self.stars[i].is_deleted() == true {
s.push_str(&format!(
"{} {} {}\n",
onegoodpt[0], onegoodpt[1], onegoodpt[2]
));
continue;
}
s.push_str(&format!(
"{} {} {}\n",
self.stars[i].pt[0], self.stars[i].pt[1], self.stars[i].pt[2]
));
}
write!(f, "{}", s).unwrap();
let mut s = String::new();
for tr in trs.iter() {
s.push_str(&format!(
"3 {} {} {}\n",
tr.v[0] - 1,
tr.v[1] - 1,
tr.v[2] - 1
));
}
write!(f, "{}", s).unwrap();
Ok(())
}
/// Returns a [`String`] containing different statistics about the triangulation.
pub fn printme(&self, withxyz: bool) -> String {
let mut s = String::from("**********\n");
// s.push_str(&format!("#pts: {}\n", self.number_pts()));
for (i, p) in self.stars.iter().enumerate() {
// s.push_str(&format!("{}: {}\n", i, self.stars[i].link));
s.push_str(&format!("{}: [", i));
for each in p.link.iter() {
s.push_str(&format!("{} - ", each));
}
s.push_str(&format!("]\n"));
if withxyz == true {
s.push_str(&format!("\t{:?}\n", self.stars[i].pt));
}
}
s.push_str("**********\n");
s
}
/// Returns the area of the Voronoi cell of `vi`.
///
/// # Arguments
///
/// * `vi` - the index of the vertex
/// * `ignore_infinity` - calculate the area even is `vi` is on the convex hull.
/// This is used by [`interpolation::NNI`] when neighbours have no area, this bounds
/// arbitrarily the area and because we take the different the interpolated value
/// is the same at the end.
pub fn voronoi_cell_area(&self, vi: usize, ignore_infinity: bool) -> Option<f64> {
if self.is_vertex_valid(vi) == false {
return None;
}
if (ignore_infinity == false) && (self.is_vertex_convex_hull(vi) == true) {
return Some(f64::INFINITY);
}
//-- process non-CH points that exists
let mut centres: Vec<Vec<f64>> = Vec::new();
let mut l = self.stars[vi].link.clone();
if l.contains_infinite_vertex() {
l.delete(0);
}
for (i, n) in l.iter().enumerate() {
let j = l.next_index(i);
centres.push(geom::circle_centre(
&self.stars[vi].pt,
&self.stars[*n].pt,
&self.stars[l[j]].pt,
));
}
//-- copy first to make circular
centres.push(vec![centres[0][0], centres[0][1]]);
let mut totalarea = 0.0_f64;
for c in centres.windows(2) {
totalarea += geom::area_triangle(&self.stars[vi].pt, &c[0], &c[1]);
}
Some(totalarea)
}
fn is_vertex_valid(&self, v: usize) -> bool {
let mut re = true;
if v >= self.stars.len() || self.stars[v].is_deleted() == true {
re = false;
}
re
}
/// Returns the (axis-aligned) bounding box of the triangulation.
pub fn get_bbox(&self) -> Vec<f64> {
let mut minx: f64 = std::f64::MAX;
let mut miny: f64 = std::f64::MAX;
let mut maxx: f64 = std::f64::MIN;
let mut maxy: f64 = std::f64::MIN;
for i in 1..self.stars.len() {
if self.is_vertex_removed(i).unwrap() == true {
continue;
}
if self.stars[i].pt[0] < minx {
minx = self.stars[i].pt[0];
}
if self.stars[i].pt[1] < miny {
miny = self.stars[i].pt[1];
}
if self.stars[i].pt[0] > maxx {
maxx = self.stars[i].pt[0];
}
if self.stars[i].pt[1] > maxy {
maxy = self.stars[i].pt[1];
}
}
vec![minx, miny, maxx, maxy]
}
/// Exaggerates vertically the z-values, used for visualisation mostly.
///
/// The value can be <1.0 to have negative exaggeration.
pub fn vertical_exaggeration(&mut self, factor: f64) {
let mut minz: f64 = std::f64::MAX;
for i in 1..self.stars.len() {
if self.stars[i].is_deleted() == true {
continue;
}
if self.stars[i].pt[2] < minz {
minz = self.stars[i].pt[2];
}
}
for i in 1..self.stars.len() {
if self.stars[i].is_deleted() == true {
continue;
}
let z2 = ((self.stars[i].pt[2] - minz) * factor) + minz;
self.stars[i].pt[2] = z2;
}
}
pub fn has_garbage(&self) -> bool {
if self.number_of_removed_vertices() > 0 {
true
} else {
false
}
}
/// Collect garbage, that is remove from memory the vertices
/// marked as removed.
///
/// Watch out: the vertices get new IDs (and thus the triangles) too. And this can
/// be a slow operation.
pub fn collect_garbage(&mut self) {
self.removed_indices.sort_unstable();
for star in self.stars.iter_mut() {
for value in star.link.0.iter_mut() {
let pos = self.removed_indices.binary_search(value).unwrap_err();
let newv = *value - pos;
*value = newv;
}
}
let mut offset = 0;
for each in &self.removed_indices {
self.stars.remove(each - offset);
offset += 1;
}
self.removed_indices.clear();
}
}
impl fmt::Display for Triangulation {
fn fmt(&self, fmt: &mut fmt::Formatter) -> fmt::Result {
fmt.write_str("======== TRIANGULATION ========\n")?;
fmt.write_str(&format!("# vertices: {:19}\n", self.number_of_vertices()))?;
fmt.write_str(&format!("# triangles: {:18}\n", self.number_of_triangles()))?;
fmt.write_str(&format!(
"# convex hull: {:16}\n",
self.number_of_vertices_on_convex_hull()
))?;
fmt.write_str(&format!("---\n"))?;
fmt.write_str(&format!("robust: {}\n", self.robust_predicates))?;
fmt.write_str(&format!("tolerance: {}\n", self.snaptol))?;
fmt.write_str("===============================\n")?;
Ok(())
}
}