STEP export: unique ids, plain curves, right revolution senses; seams on split faces

Vendored truck-stepio (writer only): an intersection curve's second surface
was written at the first one's id (duplicate ids for every boolean edge),
and revolved surfaces lost their orientation (the modeling-surface sense was
always true; inverted or mirrored revolutions came out inside out). Export
now writes intersection edges as C1 cubic B-splines through exact points,
so files carry plain 3D curves and no redundant surfaces.

Vendored truck-meshalgo: hinted parameter searches that run off a B-spline's
domain are rejected or clamped (meshing their grid never finished).

STEP import: straight revolution profiles are bounded in the meridian plane
and stop at the axis; loops on revolutions start away from the axis; face
extents include edge end points. Edges between pieces of one split surface
draw as seams.

kernel::step_validate checks a STEP file (duplicate ids, dangling
references, dense numbering); the oracle replay now exports every rebuilt
part, validates it and re-imports it with no warnings and the same volume
(all 32 rebuildable cases pass).

Co-Authored-By: Claude Opus 5.5 <[email protected]>
Claude-Session: https://claude.ai/code/session_01AS5Tn3xieewUSg3uxrMVdL
This commit is contained in:
bherbruckandClaude Opus 5.5 committed 2026-10-07 23:40:07 +00:00
1 parent 32e6b2a578
commit 08cf223ced
47 files changed
+10677 -39

No files matched your search

+113
View File
@@ -0,0 +1,113 @@
# THIS FILE IS AUTOMATICALLY GENERATED BY CARGO
#
# When uploading crates to the registry Cargo will automatically
# "normalize" Cargo.toml files for maximal compatibility
# with all versions of Cargo and also rewrite `path` dependencies
# to registry (e.g., crates.io) dependencies.
#
# If you are reading this file be aware that the original Cargo.toml
# will likely look very different (and much more reasonable).
# See Cargo.toml.orig for the original contents.
[package]
edition = "2021"
name = "truck-meshalgo"
version = "0.4.0"
build = false
autobins = false
autoexamples = false
autotests = false
autobenches = false
description = "Mesh algorigthms, include tessellations of the shape."
homepage = "https://github.com/ricosjp/truck"
readme = "README.md"
keywords = [
"truck",
"graphics",
"mesh",
]
categories = ["graphics"]
license = "Apache-2.0"
repository = "https://github.com/ricosjp/truck"
[lib]
name = "truck_meshalgo"
path = "src/lib.rs"
[dependencies.array-macro]
version = "2.1.8"
[dependencies.derive_more]
version = "0.99.18"
[dependencies.itertools]
version = "0.13.0"
[dependencies.rustc-hash]
version = "2.0.0"
[dependencies.spade]
version = "2.12.1"
optional = true
[dependencies.truck-base]
version = "0.5.0"
[dependencies.truck-geometry]
version = "0.5.0"
[dependencies.truck-polymesh]
version = "0.6.0"
[dependencies.truck-topology]
version = "0.6.0"
[dev-dependencies.proptest]
version = "1.5.0"
[dev-dependencies.rand]
version = "0.8.5"
[dev-dependencies.serde_json]
version = "1.0.128"
[dev-dependencies.truck-modeling]
version = "0.6.0"
[features]
analyzers = []
default = [
"analyzers",
"filters",
"tessellation",
"vtk",
]
filters = []
tessellation = [
"spade",
"filters",
]
vtk = ["vtkio"]
[target.'cfg(not(target_arch = "wasm32"))'.dependencies.rayon]
version = "1.10.0"
[target.'cfg(not(target_arch = "wasm32"))'.dependencies.vtkio]
version = "0.6.3"
optional = true
[target.'cfg(target_arch = "wasm32")'.dev-dependencies.getrandom]
version = "0.2.15"
features = ["js"]
+177
View File
@@ -0,0 +1,177 @@
Apache License
Version 2.0, January 2004
http://www.apache.org/licenses/
TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION
1. Definitions.
"License" shall mean the terms and conditions for use, reproduction,
and distribution as defined by Sections 1 through 9 of this document.
"Licensor" shall mean the copyright owner or entity authorized by
the copyright owner that is granting the License.
"Legal Entity" shall mean the union of the acting entity and all
other entities that control, are controlled by, or are under common
control with that entity. For the purposes of this definition,
"control" means (i) the power, direct or indirect, to cause the
direction or management of such entity, whether by contract or
otherwise, or (ii) ownership of fifty percent (50%) or more of the
outstanding shares, or (iii) beneficial ownership of such entity.
"You" (or "Your") shall mean an individual or Legal Entity
exercising permissions granted by this License.
"Source" form shall mean the preferred form for making modifications,
including but not limited to software source code, documentation
source, and configuration files.
"Object" form shall mean any form resulting from mechanical
transformation or translation of a Source form, including but
not limited to compiled object code, generated documentation,
and conversions to other media types.
"Work" shall mean the work of authorship, whether in Source or
Object form, made available under the License, as indicated by a
copyright notice that is included in or attached to the work
(an example is provided in the Appendix below).
"Derivative Works" shall mean any work, whether in Source or Object
form, that is based on (or derived from) the Work and for which the
editorial revisions, annotations, elaborations, or other modifications
represent, as a whole, an original work of authorship. For the purposes
of this License, Derivative Works shall not include works that remain
separable from, or merely link (or bind by name) to the interfaces of,
the Work and Derivative Works thereof.
"Contribution" shall mean any work of authorship, including
the original version of the Work and any modifications or additions
to that Work or Derivative Works thereof, that is intentionally
submitted to Licensor for inclusion in the Work by the copyright owner
or by an individual or Legal Entity authorized to submit on behalf of
the copyright owner. For the purposes of this definition, "submitted"
means any form of electronic, verbal, or written communication sent
to the Licensor or its representatives, including but not limited to
communication on electronic mailing lists, source code control systems,
and issue tracking systems that are managed by, or on behalf of, the
Licensor for the purpose of discussing and improving the Work, but
excluding communication that is conspicuously marked or otherwise
designated in writing by the copyright owner as "Not a Contribution."
"Contributor" shall mean Licensor and any individual or Legal Entity
on behalf of whom a Contribution has been received by Licensor and
subsequently incorporated within the Work.
2. Grant of Copyright License. Subject to the terms and conditions of
this License, each Contributor hereby grants to You a perpetual,
worldwide, non-exclusive, no-charge, royalty-free, irrevocable
copyright license to reproduce, prepare Derivative Works of,
publicly display, publicly perform, sublicense, and distribute the
Work and such Derivative Works in Source or Object form.
3. Grant of Patent License. Subject to the terms and conditions of
this License, each Contributor hereby grants to You a perpetual,
worldwide, non-exclusive, no-charge, royalty-free, irrevocable
(except as stated in this section) patent license to make, have made,
use, offer to sell, sell, import, and otherwise transfer the Work,
where such license applies only to those patent claims licensable
by such Contributor that are necessarily infringed by their
Contribution(s) alone or by combination of their Contribution(s)
with the Work to which such Contribution(s) was submitted. If You
institute patent litigation against any entity (including a
cross-claim or counterclaim in a lawsuit) alleging that the Work
or a Contribution incorporated within the Work constitutes direct
or contributory patent infringement, then any patent licenses
granted to You under this License for that Work shall terminate
as of the date such litigation is filed.
4. Redistribution. You may reproduce and distribute copies of the
Work or Derivative Works thereof in any medium, with or without
modifications, and in Source or Object form, provided that You
meet the following conditions:
(a) You must give any other recipients of the Work or
Derivative Works a copy of this License; and
(b) You must cause any modified files to carry prominent notices
stating that You changed the files; and
(c) You must retain, in the Source form of any Derivative Works
that You distribute, all copyright, patent, trademark, and
attribution notices from the Source form of the Work,
excluding those notices that do not pertain to any part of
the Derivative Works; and
(d) If the Work includes a "NOTICE" text file as part of its
distribution, then any Derivative Works that You distribute must
include a readable copy of the attribution notices contained
within such NOTICE file, excluding those notices that do not
pertain to any part of the Derivative Works, in at least one
of the following places: within a NOTICE text file distributed
as part of the Derivative Works; within the Source form or
documentation, if provided along with the Derivative Works; or,
within a display generated by the Derivative Works, if and
wherever such third-party notices normally appear. The contents
of the NOTICE file are for informational purposes only and
do not modify the License. You may add Your own attribution
notices within Derivative Works that You distribute, alongside
or as an addendum to the NOTICE text from the Work, provided
that such additional attribution notices cannot be construed
as modifying the License.
You may add Your own copyright statement to Your modifications and
may provide additional or different license terms and conditions
for use, reproduction, or distribution of Your modifications, or
for any such Derivative Works as a whole, provided Your use,
reproduction, and distribution of the Work otherwise complies with
the conditions stated in this License.
5. Submission of Contributions. Unless You explicitly state otherwise,
any Contribution intentionally submitted for inclusion in the Work
by You to the Licensor shall be under the terms and conditions of
this License, without any additional terms or conditions.
Notwithstanding the above, nothing herein shall supersede or modify
the terms of any separate license agreement you may have executed
with Licensor regarding such Contributions.
6. Trademarks. This License does not grant permission to use the trade
names, trademarks, service marks, or product names of the Licensor,
except as required for reasonable and customary use in describing the
origin of the Work and reproducing the content of the NOTICE file.
7. Disclaimer of Warranty. Unless required by applicable law or
agreed to in writing, Licensor provides the Work (and each
Contributor provides its Contributions) on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or
implied, including, without limitation, any warranties or conditions
of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A
PARTICULAR PURPOSE. You are solely responsible for determining the
appropriateness of using or redistributing the Work and assume any
risks associated with Your exercise of permissions under this License.
8. Limitation of Liability. In no event and under no legal theory,
whether in tort (including negligence), contract, or otherwise,
unless required by applicable law (such as deliberate and grossly
negligent acts) or agreed to in writing, shall any Contributor be
liable to You for damages, including any direct, indirect, special,
incidental, or consequential damages of any character arising as a
result of this License or out of the use or inability to use the
Work (including but not limited to damages for loss of goodwill,
work stoppage, computer failure or malfunction, or any and all
other commercial damages or losses), even if such Contributor
has been advised of the possibility of such damages.
9. Accepting Warranty or Additional Liability. While redistributing
the Work or Derivative Works thereof, You may choose to offer,
and charge a fee for, acceptance of support, warranty, indemnity,
or other liability obligations and/or rights consistent with this
License. However, in accepting such obligations, You may act only
on Your own behalf and on Your sole responsibility, not on behalf
of any other Contributor, and only if You agree to indemnify,
defend, and hold each Contributor harmless for any liability
incurred by, or claims asserted against, such Contributor by reason
of your accepting any such warranty or additional liability.
END OF TERMS AND CONDITIONS
+68
View File
@@ -0,0 +1,68 @@
# truck-meshalgo
[![Crates.io](https://img.shields.io/crates/v/truck-meshalgo.svg)](https://crates.io/crates/truck-meshalgo) [![Docs.rs](https://docs.rs/truck-meshalgo/badge.svg)](https://docs.rs/truck-meshalgo)
Mesh algorigthms, include tessellations of the shape.
## Sample Codes
### filleted-cube
An experiment to decompose a mesh into elements for future NURBS-shape approximation of the mesh.
Contains `doc(hidden)` methods.
- Input: filleted_cube.obj
- Output: planes.obj, lower.obj, upper.obj
### irregular-sphere
Add the normal to the sphere containing
the irregular normal generated from the NURBS containing the critical point.
- Input: irregular_sphere.obj
- Output: regular_sphere.obj
### octahedron-subdivision
Apply loop subdivision to regular octahedron.
- Input: hardcoded octahedron
- Output: octahedron.obj, subdivision-octahedron.obj
### requadrangulate-buddha
A benchmark that reads in heavy mesh data, applies triangulation and quadrangulation, and writes it out.
- Input: happy-buddha.obj
- Output: requadrangulated-buddha.obj
### smoothing-bunny
Adds smooth normals to the stanford bunny.
- Input: bunny.obj
- Output: smooth_bunny.obj
### splitting-sample
An experiment to decompose a mesh into elements for future NURBS-shape approximation of the mesh.
- Input: sample.obj
- Output: planes_parts_#.obj, others_parts_#.obj
### teapot
Adds smooth normals to and quadrangulate the famous teapot.
- Input: teapot.obj
- Output: quaded_pot.obj
### tessellate-shape
Tessellate a shape and output an obj file.
```bash
usage: tessellate-shape <input json file> <output json file>
```
The default `<output file>` is output.obj.
+275
View File
@@ -0,0 +1,275 @@
use super::*;
use array_macro::array;
/// Find collisions between two polygon meshes and extract interference lines.
///
/// # Details
///
/// The algorithm for mesh extraction is as follows.
/// All the triangles are projected toward the axis to form intervals.
/// By looking at the overlap of these intervals, we can narrow down the interfering triangles.
/// The overlap of the intervals is obtained by sorting the endpoints.
///
/// # Remarks
///
/// We see that the surfaces that are in contact are not interfering with each other.
/// Therefore, polygon meshes included in one plane are not interfering.
pub trait Collision {
/// If `self` and `other` collide, then returns only one interference line.
/// Otherwise, returns `None`.
fn collide_with(&self, other: &PolygonMesh) -> Option<(Point3, Point3)>;
/// Extract all interference lines between `self` and `other`.
/// # Remarks
/// The results is not arranged so that included lines make continuous maximal polyline curve.
fn extract_interference(&self, other: &PolygonMesh) -> Vec<(Point3, Point3)>;
}
impl Collision for PolygonMesh {
#[inline(always)]
fn collide_with(&self, other: &PolygonMesh) -> Option<(Point3, Point3)> {
are_colliding(self, other)
}
#[inline(always)]
fn extract_interference(&self, other: &PolygonMesh) -> Vec<(Point3, Point3)> {
collision(self, other)
}
}
#[repr(u8)]
#[derive(Clone, Copy, Debug, PartialEq, Eq, PartialOrd, Ord)]
enum EndPointType {
Front,
Back,
}
#[derive(Clone, Copy, Debug, PartialEq, PartialOrd)]
struct EndPoint {
entity: f64,
r#type: EndPointType,
segnum: usize,
index: usize,
}
impl EndPoint {
#[inline(always)]
fn new(entity: f64, r#type: EndPointType, segnum: usize, index: usize) -> EndPoint {
EndPoint {
entity,
r#type,
segnum,
index,
}
}
#[inline(always)]
fn from_seg(seg: (f64, f64), segnum: usize, index: usize) -> Vec<EndPoint> {
vec![
EndPoint::new(seg.0, EndPointType::Front, segnum, index),
EndPoint::new(seg.1, EndPointType::Back, segnum, index),
]
}
}
fn tri_to_seg(tri: [Point3; 3], unit: Vector3) -> (f64, f64) {
let a = tri[0].to_vec().dot(unit);
let b = tri[1].to_vec().dot(unit);
let c = tri[2].to_vec().dot(unit);
(f64::min(f64::min(a, b), c), f64::max(f64::max(a, b), c))
}
fn sorted_endpoints<I, J>(iter0: I, iter1: J, unit: Vector3) -> Vec<EndPoint>
where
I: IntoIterator<Item = [Point3; 3]>,
J: IntoIterator<Item = [Point3; 3]>, {
let mut res: Vec<EndPoint> = iter0
.into_iter()
.enumerate()
.filter(|(_, tri)| !(tri[1] - tri[0]).cross(tri[2] - tri[0]).so_small())
.flat_map(|(i, tri)| EndPoint::from_seg(tri_to_seg(tri, unit), 0, i))
.chain(
iter1
.into_iter()
.enumerate()
.filter(|(_, tri)| !(tri[1] - tri[0]).cross(tri[2] - tri[0]).so_small())
.flat_map(|(i, tri)| EndPoint::from_seg(tri_to_seg(tri, unit), 1, i)),
)
.collect();
res.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Greater));
res
}
fn colliding_segment_pairs(sort_endpoints: Vec<EndPoint>) -> impl Iterator<Item = (usize, usize)> {
let mut current = [Vec::<usize>::new(), Vec::<usize>::new()];
sort_endpoints
.into_iter()
.filter_map(
move |EndPoint {
r#type,
segnum,
index,
..
}| match r#type {
EndPointType::Front => {
current[segnum].push(index);
Some(current[1 - segnum].clone().into_iter().map(move |i| {
if segnum == 0 {
(index, i)
} else {
(i, index)
}
}))
}
EndPointType::Back => {
let i = current[segnum]
.iter()
.enumerate()
.find(|(_, idx)| **idx == index)
.unwrap()
.0;
current[segnum].swap_remove(i);
None
}
},
)
.flatten()
}
fn disjoint_bdbs(tri0: [Point3; 3], tri1: [Point3; 3]) -> bool {
let bdb0: BoundingBox<Point3> = tri0.iter().collect();
let bdb1: BoundingBox<Point3> = tri1.iter().collect();
bdb0.max()[0] < bdb1.min()[0]
|| bdb1.max()[0] < bdb0.min()[0]
|| bdb0.max()[1] < bdb1.min()[1]
|| bdb1.max()[1] < bdb0.min()[1]
|| bdb0.max()[2] < bdb1.min()[2]
|| bdb1.max()[2] < bdb0.min()[2]
}
fn collide_seg_triangle(seg: [Point3; 2], tri: [Point3; 3]) -> Option<Point3> {
let ab = tri[1] - tri[0];
let bc = tri[2] - tri[1];
let ca = tri[0] - tri[2];
let nor = ab.cross(ca);
if nor.so_small() {
return None;
}
let ap = seg[0] - tri[0];
let aq = seg[1] - tri[0];
let dotapnor = ap.dot(nor);
let dotaqnor = aq.dot(nor);
if dotapnor * dotaqnor > 0.0 {
return None;
}
let h = seg[0] + dotapnor / (dotapnor - dotaqnor) * (seg[1] - seg[0]);
if f64::signum(ab.cross(nor).dot(h - tri[0]) + TOLERANCE2)
+ f64::signum(bc.cross(nor).dot(h - tri[1]) + TOLERANCE2)
+ f64::signum(ca.cross(nor).dot(h - tri[2]) + TOLERANCE2)
>= 2.0
{
Some(h)
} else {
None
}
}
fn collide_triangles(tri0: [Point3; 3], tri1: [Point3; 3]) -> Option<(Point3, Point3)> {
let mut tuple = (None, None);
[
collide_seg_triangle([tri0[0], tri0[1]], tri1),
collide_seg_triangle([tri0[1], tri0[2]], tri1),
collide_seg_triangle([tri0[2], tri0[0]], tri1),
collide_seg_triangle([tri1[0], tri1[1]], tri0),
collide_seg_triangle([tri1[1], tri1[2]], tri0),
collide_seg_triangle([tri1[2], tri1[0]], tri0),
]
.into_iter()
.for_each(|pt| match tuple {
(None, _) => tuple.0 = pt,
(Some(_), None) => tuple.1 = pt,
(Some(ref mut p), Some(ref mut q)) => {
if let Some(pt) = pt {
let dist0 = pt.distance2(*p);
let dist1 = pt.distance2(*q);
let dist2 = p.distance2(*q);
if dist2 < dist0 {
*q = pt;
} else if dist2 < dist1 {
*p = pt;
}
}
}
});
match tuple {
(Some(a), Some(b)) => Some((a, b)),
_ => None,
}
}
fn make_pos_tri(poly: &PolygonMesh, face: [StandardVertex; 3]) -> [Point3; 3] {
array![i => poly.positions()[face[i].pos]; 3]
}
fn collision(poly0: &PolygonMesh, poly1: &PolygonMesh) -> Vec<(Point3, Point3)> {
let unit = match poly0.positions().is_empty() {
true => return Vec::new(),
false => hash::take_one_unit(poly0.positions()[0]),
};
let tris0 = poly0.faces().triangle_iter().collect::<Vec<_>>();
let tris1 = poly1.faces().triangle_iter().collect::<Vec<_>>();
let iter0 = tris0.iter().map(|face| make_pos_tri(poly0, *face));
let iter1 = tris1.iter().map(|face| make_pos_tri(poly1, *face));
colliding_segment_pairs(sorted_endpoints(iter0, iter1, unit))
.filter_map(|(idx0, idx1)| {
let tri0 = make_pos_tri(poly0, tris0[idx0]);
let tri1 = make_pos_tri(poly1, tris1[idx1]);
match disjoint_bdbs(tri0, tri1) {
true => None,
false => collide_triangles(tri0, tri1),
}
})
.collect()
}
fn are_colliding(poly0: &PolygonMesh, poly1: &PolygonMesh) -> Option<(Point3, Point3)> {
let unit = match poly0.positions().is_empty() {
true => return None,
false => hash::take_one_unit(poly0.positions()[0]),
};
let tris0 = poly0.faces().triangle_iter().collect::<Vec<_>>();
let tris1 = poly1.faces().triangle_iter().collect::<Vec<_>>();
let iter0 = tris0.iter().map(|face| make_pos_tri(poly0, *face));
let iter1 = tris1.iter().map(|face| make_pos_tri(poly1, *face));
colliding_segment_pairs(sorted_endpoints(iter0, iter1, unit)).find_map(|(idx0, idx1)| {
let tri0 = make_pos_tri(poly0, tris0[idx0]);
let tri1 = make_pos_tri(poly1, tris1[idx1]);
match disjoint_bdbs(tri0, tri1) {
true => None,
false => collide_triangles(tri0, tri1),
}
})
}
#[test]
fn collide_triangles_test() {
let tri0 = [
Point3::origin(),
Point3::new(1.0, 0.0, 0.0),
Point3::new(0.0, 1.0, 0.0),
];
let tri1 = [
Point3::new(0.0, 0.0, -1.0),
Point3::new(-1.0, -1.0, 1.0),
Point3::new(1.0, 1.0, 1.0),
];
assert!(collide_triangles(tri0, tri1).is_some());
let tri0 = [
Point3::new(0.0, 0.0, 0.0),
Point3::new(1.0, 0.0, 0.0),
Point3::new(0.0, 1.0, 0.0),
];
let tri1 = [
Point3::new(0.0, 0.0, 0.5),
Point3::new(1.0, 0.0, 1.0),
Point3::new(1.0, 1.0, 1.0),
];
assert!(collide_triangles(tri0, tri1).is_none());
}
+140
View File
@@ -0,0 +1,140 @@
use super::*;
use derive_more::{Deref, DerefMut};
#[derive(Clone, Copy, Debug)]
struct Ray {
origin: Point3,
direction: Vector3,
}
#[derive(Clone, Copy, Debug, Deref, DerefMut)]
struct Triangle([Point3; 3]);
impl Triangle {
#[inline(always)]
fn normal(self) -> Vector3 { (self[1] - self[0]).cross(self[2] - self[0]).normalize() }
fn is_crossing(self, ray: Ray) -> bool {
let a = self[0] - self[1];
let b = self[0] - self[2];
let mat = Matrix3::from_cols(a, b, ray.direction);
if mat.determinant().so_small() {
false
} else {
let inv = mat.invert().unwrap();
let uvt = inv * (self[0] - ray.origin);
uvt[0] > 0.0 && uvt[1] > 0.0 && uvt[0] + uvt[1] < 1.0 && uvt[2] > 0.0
}
}
}
/// whether a point is in a domain rounded by a closed polygon.
pub trait IncludingPointInDomain {
/// Count signed number of faces crossing ray with origin `point` and direction `ray_direction`.
/// Counter increase if the dot product of the ray and the normal of a face is positive,
/// and decrease if it is negative.
///
/// # Examples
/// ```
/// use truck_meshalgo::prelude::*;
/// let simplex = PolygonMesh::new(
/// StandardAttributes {
/// positions: vec![
/// Point3::new(0.0, 0.0, 0.0),
/// Point3::new(1.0, 0.0, 0.0),
/// Point3::new(0.0, 1.0, 0.0),
/// Point3::new(0.0, 0.0, 1.0),
/// ],
/// ..Default::default()
/// },
/// Faces::from_iter(vec![
/// [0, 2, 1],
/// [0, 1, 3],
/// [0, 3, 2],
/// [1, 2, 3],
/// ]),
/// );
///
/// assert_eq!(simplex.signed_crossing_faces(Point3::new(0.1, 0.1, 0.1), Vector3::unit_x()), 1);
/// assert_eq!(simplex.signed_crossing_faces(Point3::new(-0.1, 0.1, 0.1), Vector3::unit_x()), 0);
/// ```
fn signed_crossing_faces(&self, point: Point3, ray_direction: Vector3) -> isize;
/// whether `point` is in a domain rounded by a closed polygon.
///
/// # Examples
/// ```
/// use truck_meshalgo::prelude::*;
/// let simplex = PolygonMesh::new(
/// StandardAttributes {
/// positions: vec![
/// Point3::new(0.0, 0.0, 0.0),
/// Point3::new(1.0, 0.0, 0.0),
/// Point3::new(0.0, 1.0, 0.0),
/// Point3::new(0.0, 0.0, 1.0),
/// ],
/// ..Default::default()
/// },
/// Faces::from_iter(vec![
/// [0, 2, 1],
/// [0, 1, 3],
/// [0, 3, 2],
/// [1, 2, 3],
/// ]),
/// );
///
/// assert!(simplex.inside(Point3::new(0.1, 0.1, 0.1)));
/// assert!(!simplex.inside(Point3::new(-0.1, 0.1, 0.1)));
/// ```
fn inside(&self, point: Point3) -> bool;
}
impl IncludingPointInDomain for PolygonMesh {
fn signed_crossing_faces(&self, point: Point3, ray_direction: Vector3) -> isize {
let ray = Ray {
origin: point,
direction: ray_direction,
};
self.face_iter().fold(0, |mut counter, face| {
for i in 2..face.len() {
let tri = Triangle([
self.positions()[face[0].pos],
self.positions()[face[i - 1].pos],
self.positions()[face[i].pos],
]);
if tri.is_crossing(ray) {
counter += f64::signum(tri.normal().dot(ray.direction)) as isize;
}
}
counter
})
}
fn inside(&self, point: Point3) -> bool {
let dir = hash::take_one_unit(point);
self.signed_crossing_faces(point, dir) >= 1
}
}
#[test]
fn inside100() {
let positions = vec![
Point3::new(0.0, 0.0, 0.0),
Point3::new(1.0, 0.0, 0.0),
Point3::new(0.0, 1.0, 0.0),
Point3::new(0.0, 0.0, 1.0),
];
let faces: Faces = vec![[0, 2, 1], [0, 1, 3], [0, 3, 2], [1, 2, 3]]
.into_iter()
.collect();
let simplex = PolygonMesh::new(
StandardAttributes {
positions,
..Default::default()
},
faces,
);
for _ in 0..100 {
assert!(simplex.inside(Point3::new(0.1, 0.1, 0.1)));
assert!(!simplex.inside(Point3::new(-0.1, 0.1, 0.1)));
}
}
+17
View File
@@ -0,0 +1,17 @@
use crate::*;
mod collision;
mod in_out_judge;
mod point_cloud;
mod splitting;
mod topology;
mod volume;
pub use collision::Collision;
pub use in_out_judge::IncludingPointInDomain;
pub use point_cloud::WithPointCloud;
pub use splitting::ExperimentalSplitters;
pub use splitting::Splitting;
pub use topology::Topology;
pub use truck_topology::shell::ShellCondition;
pub use volume::CalcVolume;
@@ -0,0 +1,283 @@
use super::*;
use array_macro::array;
#[derive(Clone, Debug)]
pub struct HashedPointCloud {
space: Vec<Vec<Point3>>,
size: [usize; 3],
range: [[f64; 2]; 3],
num_points: usize,
}
#[allow(dead_code)]
impl HashedPointCloud {
#[inline(always)]
pub fn new(size: [usize; 3], range: [[f64; 2]; 3]) -> Self {
HashedPointCloud {
space: vec![Vec::new(); size[0] * size[1] * size[2]],
size,
range,
num_points: 0,
}
}
#[inline(always)]
pub fn from_points<'a, I>(points: I, inf: f64) -> Self
where
I: IntoIterator + Clone,
I::IntoIter: Iterator<Item = &'a Point3>, {
let mut bdb = BoundingBox::<Point3>::new();
let len = points.clone().into_iter().copied().fold(0, |counter, pt| {
bdb.push(pt);
counter + 1
});
let size = 1.0 + f64::powf(len as f64, 1.0 / 3.0);
let size = bdb.diagonal().map(|a| {
if a / size < inf {
1 + (a / inf) as usize
} else {
size as usize
}
});
let mut res =
HashedPointCloud::new(size.into(), array![i => [bdb.min()[i], bdb.max()[i]]; 3]);
points.into_iter().for_each(|pt| res.push(*pt));
res
}
#[inline(always)]
pub fn size(&self) -> [usize; 3] { self.size }
#[inline(always)]
pub fn push(&mut self, point: Point3) {
self.num_points += 1;
let idx = point.hash(self);
self[idx].push(point);
}
#[inline(always)]
pub fn distance(&self, t: impl DistanceWithPointCloud) -> f64 { t.distance(self) }
#[inline(always)]
pub fn distance2(&self, t: impl DistanceWithPointCloud) -> f64 { t.distance2(self) }
#[inline(always)]
pub fn is_colliding(&self, t: impl DistanceWithPointCloud, tol: f64) -> bool {
t.is_colliding(self, tol)
}
}
impl std::ops::Index<[usize; 3]> for HashedPointCloud {
type Output = Vec<Point3>;
#[inline(always)]
fn index(&self, idcs: [usize; 3]) -> &Vec<Point3> {
&self.space[(idcs[0] * self.size[1] + idcs[1]) * self.size[2] + idcs[2]]
}
}
impl std::ops::IndexMut<[usize; 3]> for HashedPointCloud {
#[inline(always)]
fn index_mut(&mut self, idcs: [usize; 3]) -> &mut Vec<Point3> {
&mut self.space[(idcs[0] * self.size[1] + idcs[1]) * self.size[2] + idcs[2]]
}
}
pub trait SpaceHash {
fn hash(self, space: &HashedPointCloud) -> [usize; 3];
}
impl SpaceHash for Point3 {
#[inline(always)]
fn hash(self, space: &HashedPointCloud) -> [usize; 3] {
let x = (self[0] - space.range[0][0]) / (space.range[0][1] - space.range[0][0]);
let ix = f64::clamp(x * space.size[0] as f64, 0.0, space.size[0] as f64 - 1.0) as usize;
let y = (self[1] - space.range[1][0]) / (space.range[1][1] - space.range[1][0]);
let iy = f64::clamp(y * space.size[1] as f64, 0.0, space.size[1] as f64 - 1.0) as usize;
let z = (self[2] - space.range[2][0]) / (space.range[2][1] - space.range[2][0]);
let iz = f64::clamp(z * space.size[2] as f64, 0.0, space.size[2] as f64 - 1.0) as usize;
[ix, iy, iz]
}
}
impl SpaceHash for usize {
#[inline(always)]
fn hash(self, space: &HashedPointCloud) -> [usize; 3] {
[
self / (space.size[1] * space.size[2]),
self % (space.size[1] * space.size[2]) / space.size[2],
self % space.size[2],
]
}
}
pub trait DistanceWithPointCloud: Sized {
fn distance2(&self, space: &HashedPointCloud) -> f64;
fn distance(&self, space: &HashedPointCloud) -> f64 { f64::sqrt(self.distance2(space)) }
fn is_colliding(&self, space: &HashedPointCloud, tol: f64) -> bool {
nonpositive_tolerance!(tol, 0.0);
self.distance2(space) < tol * tol
}
}
impl DistanceWithPointCloud for Point3 {
fn distance2(&self, space: &HashedPointCloud) -> f64 {
let idcs = self.hash(space);
let closure = |dist2: f64, pt: &Point3| f64::min(dist2, MetricSpace::distance2(*self, *pt));
let mut dist2 = space[idcs].iter().fold(f64::INFINITY, closure);
if idcs[0] > 0 {
dist2 = space[[idcs[0] - 1, idcs[1], idcs[2]]]
.iter()
.fold(dist2, closure);
}
if idcs[0] + 1 < space.size[0] {
dist2 = space[[idcs[0] + 1, idcs[1], idcs[2]]]
.iter()
.fold(dist2, closure);
}
if idcs[1] > 0 {
dist2 = space[[idcs[0], idcs[1] - 1, idcs[2]]]
.iter()
.fold(dist2, closure);
}
if idcs[1] + 1 < space.size[1] {
dist2 = space[[idcs[0], idcs[1] + 1, idcs[2]]]
.iter()
.fold(dist2, closure);
}
if idcs[2] > 0 {
dist2 = space[[idcs[0], idcs[1], idcs[2] - 1]]
.iter()
.fold(dist2, closure);
}
if idcs[2] + 1 < space.size[2] {
dist2 = space[[idcs[0], idcs[1], idcs[2] + 1]]
.iter()
.fold(dist2, closure);
}
dist2
}
}
impl DistanceWithPointCloud for [Point3; 3] {
fn distance2(&self, space: &HashedPointCloud) -> f64 {
let bdd: BoundingBox<Point3> = self.iter().collect();
let mut range: [[usize; 3]; 2] = [bdd.min().hash(space), bdd.max().hash(space)];
range[0][0] = usize::max(range[0][0], 1) - 1;
range[0][1] = usize::max(range[0][1], 1) - 1;
range[0][2] = usize::max(range[0][2], 1) - 1;
range[1][0] = usize::min(range[1][0] + 1, space.size[0] - 1);
range[1][1] = usize::min(range[1][1] + 1, space.size[1] - 1);
range[1][2] = usize::min(range[1][2] + 1, space.size[2] - 1);
let mut dist2 = f64::INFINITY;
(range[0][0]..=range[1][0]).for_each(|ix| {
(range[0][1]..=range[1][1]).for_each(|iy| {
(range[0][2]..=range[1][2]).for_each(|iz| {
dist2 = space[[ix, iy, iz]].iter().fold(dist2, |dist2, pt| {
f64::min(dist2, distance2_point_triangle(*pt, *self))
});
})
})
});
dist2
}
}
impl<'a> DistanceWithPointCloud for &'a PolygonMesh {
fn distance2(&self, space: &HashedPointCloud) -> f64 {
let dist2 = self.faces().triangle_iter().fold(-1.0, |dist2, tri| {
let tri = array![i => self.positions()[tri[i].pos]; 3];
f64::max(dist2, tri.distance2(space))
});
match dist2 < 0.0 {
true => f64::INFINITY,
false => dist2,
}
}
fn is_colliding(&self, space: &HashedPointCloud, tol: f64) -> bool {
nonpositive_tolerance!(tol, 0.0);
self.faces().triangle_iter().any(|tri| {
let tri = array![i => self.positions()[tri[i].pos]; 3];
tri.distance2(space) < tol * tol
})
}
}
#[test]
fn distance_point_triangle_test() {
let triangle = [
Point3::new(0.0, 0.0, 0.0),
Point3::new(1.0, 0.0, 0.0),
Point3::new(0.0, 1.0, 0.0),
];
assert_eq!(
12.45 * 12.45,
distance2_point_triangle(Point3::new(0.25, 0.25, 12.45), triangle)
);
assert_eq!(
12.45 * 12.45,
distance2_point_triangle(Point3::new(0.25, 0.25, -12.45), triangle)
);
assert_eq!(
2.0,
distance2_point_triangle(Point3::new(-1.0, 0.5, -1.0), triangle)
);
assert_eq!(
2.0,
distance2_point_triangle(Point3::new(0.5, -1.0, 1.0), triangle)
);
assert_eq!(
1.5,
distance2_point_triangle(Point3::new(1.0, 1.0, 1.0), triangle)
);
assert_eq!(
3.0,
distance2_point_triangle(Point3::new(-1.0, -1.0, -1.0), triangle)
);
}
#[cfg(test)]
fn exec_space_division_distance() {
const SPACE_SIZE: f64 = 100.0;
const NUM_POINTS: usize = 1000;
const NUM_TRIANGLES: usize = 10;
const TRIANGLE_DISPLACEMENT: f64 = 1.0;
let points = (0..NUM_POINTS)
.map(|_| {
Point3::new(
SPACE_SIZE * rand::random::<f64>() - SPACE_SIZE / 2.0,
SPACE_SIZE * rand::random::<f64>() - SPACE_SIZE / 2.0,
SPACE_SIZE * rand::random::<f64>() - SPACE_SIZE / 2.0,
)
})
.collect::<Vec<_>>();
let triangles = (0..NUM_TRIANGLES)
.map(|_| {
let pt = Point3::new(
SPACE_SIZE * rand::random::<f64>() - SPACE_SIZE / 2.0,
SPACE_SIZE * rand::random::<f64>() - SPACE_SIZE / 2.0,
SPACE_SIZE * rand::random::<f64>() - SPACE_SIZE / 2.0,
);
[
pt,
pt + Vector3::new(
TRIANGLE_DISPLACEMENT * 2.0 * rand::random::<f64>() - TRIANGLE_DISPLACEMENT,
TRIANGLE_DISPLACEMENT * 2.0 * rand::random::<f64>() - TRIANGLE_DISPLACEMENT,
TRIANGLE_DISPLACEMENT * 2.0 * rand::random::<f64>() - TRIANGLE_DISPLACEMENT,
),
pt + Vector3::new(
TRIANGLE_DISPLACEMENT * 2.0 * rand::random::<f64>() - TRIANGLE_DISPLACEMENT,
TRIANGLE_DISPLACEMENT * 2.0 * rand::random::<f64>() - TRIANGLE_DISPLACEMENT,
TRIANGLE_DISPLACEMENT * 2.0 * rand::random::<f64>() - TRIANGLE_DISPLACEMENT,
),
]
})
.collect::<Vec<_>>();
let hashed = HashedPointCloud::from_points(&points, 1.0);
let dist_0 = triangles.iter().fold(std::f64::INFINITY, |dist, triangle| {
f64::min(dist, hashed.distance2(*triangle))
});
let dist_1 = points.iter().fold(std::f64::INFINITY, |dist2, pt| {
triangles.iter().fold(dist2, |dist2, triangle| {
f64::min(dist2, distance2_point_triangle(*pt, *triangle))
})
});
assert_eq!(dist_0, dist_1);
}
#[test]
fn space_division_distance() { (0..10).for_each(|_| exec_space_division_distance()) }
+84
View File
@@ -0,0 +1,84 @@
use super::*;
mod hashed_point_cloud;
use hashed_point_cloud::HashedPointCloud;
mod sort_end_points;
/// Investigates positional relation between polygon mesh and point cloud.
pub trait WithPointCloud {
/// Whether all faces of the polygon mesh `self` has intersection with the neighborhood of `point_cloud`.
/// # Arguments
/// - `tol`: the radius of the neighborhoods of points in point cloud.
/// # Panics
/// `tol` must be more than `TOLERANCE`.
/// # Examples
/// ```
/// use truck_meshalgo::prelude::*;
/// let mesh = PolygonMesh::new(
/// StandardAttributes {
/// positions: vec![
/// Point3::new(0.0, 0.0, 0.0),
/// Point3::new(1.0, 0.0, 0.0),
/// Point3::new(0.0, 1.0, 0.0),
/// Point3::new(0.0, 0.0, 2.0),
/// Point3::new(1.0, 0.0, 2.0),
/// Point3::new(0.0, 1.0, 2.0),
/// ],
/// ..Default::default()
/// },
/// Faces::from_iter(vec![[0, 1, 2], [3, 4, 5]]),
/// );
///
/// let mut point_cloud = vec![Point3::new(0.25, 0.25, 0.0)];
/// assert!(!mesh.is_clung_to_by(&point_cloud, 0.001));
/// point_cloud.push(Point3::new(0.25, 0.25, 2.0));
/// assert!(mesh.is_clung_to_by(&point_cloud, 0.001));
/// ```
fn is_clung_to_by(&self, point_cloud: &[Point3], tol: f64) -> bool;
/// Whether the neighborhood of the polygon mesh `self` includes `point_cloud`.
/// # Panics
/// `tol` must be more than `TOLERANCE`.
fn neighborhood_include(&self, point_cloud: &[Point3], tol: f64) -> bool;
/// Whether the polygon mesh `self` and `point_cloud` collides.
/// # Panics
/// `tol` must be more than `TOLERANCE`.
fn collide_with_neiborhood_of(&self, point_cloud: &[Point3], tol: f64) -> bool;
}
impl WithPointCloud for PolygonMesh {
#[inline(always)]
fn is_clung_to_by(&self, point_cloud: &[Point3], tol: f64) -> bool {
nonpositive_tolerance!(tol);
HashedPointCloud::from_points(point_cloud, tol * 2.0).distance2(self) < tol * tol
}
#[inline(always)]
fn collide_with_neiborhood_of(&self, point_cloud: &[Point3], tol: f64) -> bool {
HashedPointCloud::from_points(point_cloud, tol * 2.0).is_colliding(self, tol)
}
#[inline(always)]
fn neighborhood_include(&self, point_cloud: &[Point3], tol: f64) -> bool {
sort_end_points::pointcloud_in_polygon_neighborhood(self, point_cloud, tol)
}
}
// https://iquilezles.org/www/articles/distfunctions/distfunctions.htm
fn distance2_point_triangle(point: Point3, triangle: [Point3; 3]) -> f64 {
let ab = triangle[1] - triangle[0];
let ap = point - triangle[0];
let bc = triangle[2] - triangle[1];
let bp = point - triangle[1];
let ca = triangle[0] - triangle[2];
let cp = point - triangle[2];
let nor = ab.cross(ca);
let coef = f64::signum(ab.cross(nor).dot(ap))
+ f64::signum(bc.cross(nor).dot(bp))
+ f64::signum(ca.cross(nor).dot(cp));
if coef < 2.0 || nor.magnitude().so_small() {
let a = (ap - ab * f64::clamp(ab.dot(ap) / ab.dot(ab), 0.0, 1.0)).magnitude2();
let b = (bp - bc * f64::clamp(bc.dot(bp) / bc.dot(bc), 0.0, 1.0)).magnitude2();
let c = (cp - ca * f64::clamp(ca.dot(cp) / ca.dot(ca), 0.0, 1.0)).magnitude2();
f64::min(f64::min(a, b), c)
} else {
nor.dot(ap) * nor.dot(ap) / nor.magnitude2()
}
}
@@ -0,0 +1,121 @@
use super::*;
use array_macro::array;
#[repr(u8)]
#[derive(Clone, Copy, Debug, PartialEq, Eq, PartialOrd, Ord)]
enum EndPointType {
Front,
Middle,
Back,
}
#[derive(Clone, Copy, Debug, PartialEq, PartialOrd)]
struct EndPoint {
entity: f64,
r#type: EndPointType,
index: usize,
}
impl EndPoint {
#[inline(always)]
fn new(entity: f64, r#type: EndPointType, index: usize) -> EndPoint {
EndPoint {
entity,
r#type,
index,
}
}
#[inline(always)]
fn from_seg(seg: (f64, f64), index: usize) -> Vec<EndPoint> {
vec![
EndPoint::new(seg.0, EndPointType::Front, index),
EndPoint::new(seg.1, EndPointType::Back, index),
]
}
}
fn tri_to_seg(tri: [Point3; 3], unit: Vector3, tol: f64) -> (f64, f64) {
let a = tri[0].to_vec().dot(unit);
let b = tri[1].to_vec().dot(unit);
let c = tri[2].to_vec().dot(unit);
(
f64::min(f64::min(a, b), c) - tol,
f64::max(f64::max(a, b), c) + tol,
)
}
fn sorted_endpoints<'a, I, J>(iter0: I, iter1: J, unit: Vector3, tol: f64) -> Vec<EndPoint>
where
I: IntoIterator<Item = [Point3; 3]>,
J: IntoIterator<Item = &'a Point3>, {
let mut res: Vec<EndPoint> = iter0
.into_iter()
.enumerate()
.filter(|(_, tri)| !(tri[1] - tri[0]).cross(tri[2] - tri[0]).so_small())
.flat_map(|(i, tri)| EndPoint::from_seg(tri_to_seg(tri, unit, tol), i))
.chain(
iter1
.into_iter()
.enumerate()
.map(|(i, point)| EndPoint::new(point.to_vec().dot(unit), EndPointType::Middle, i)),
)
.collect();
res.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Greater));
res
}
fn sorted_endpoints_by_polymesh_points(
polygon: &PolygonMesh,
points: &[Point3],
tol: f64,
) -> Vec<EndPoint> {
let unit = if !polygon.positions().is_empty() {
hash::take_one_unit(polygon.positions()[0])
} else if !points.is_empty() {
hash::take_one_unit(points[0])
} else {
return Vec::new();
};
sorted_endpoints(
polygon
.faces()
.triangle_iter()
.map(|tri| array![i => polygon.positions()[tri[i].pos]; 3]),
points.iter(),
unit,
tol,
)
}
pub fn pointcloud_in_polygon_neighborhood(
polygon: &PolygonMesh,
points: &[Point3],
tol: f64,
) -> bool {
nonpositive_tolerance!(tol, 0.0);
let mut current = Vec::new();
let triangles = polygon.faces().triangle_iter().collect::<Vec<_>>();
sorted_endpoints_by_polymesh_points(polygon, points, tol)
.into_iter()
.all(move |EndPoint { r#type, index, .. }| match r#type {
EndPointType::Front => {
current.push(index);
true
}
EndPointType::Back => {
let i = current
.iter()
.enumerate()
.find(|(_, idx)| **idx == index)
.unwrap()
.0;
current.swap_remove(i);
true
}
EndPointType::Middle => current.iter().any(|i| {
let tri = triangles[*i];
let tri = array![i => polygon.positions()[tri[i].pos]; 3];
distance2_point_triangle(points[index], tri) < tol * tol
}),
})
}
+343
View File
@@ -0,0 +1,343 @@
use super::*;
use std::f64::consts::PI;
/// Splitting the faces into several clusters.
pub trait Splitting {
/// Creates a sub mesh by the face indices.
/// # Examples
/// ```
/// use truck_polymesh::*;
/// use truck_meshalgo::analyzers::*;
///
/// // cube
/// let mesh = PolygonMesh::new(
/// StandardAttributes {
/// positions: vec![
/// Point3::new(0.0, 0.0, 0.0),
/// Point3::new(1.0, 0.0, 0.0),
/// Point3::new(1.0, 1.0, 0.0),
/// Point3::new(0.0, 1.0, 0.0),
/// Point3::new(0.0, 0.0, 1.0),
/// Point3::new(1.0, 0.0, 1.0),
/// Point3::new(1.0, 1.0, 1.0),
/// Point3::new(0.0, 1.0, 1.0),
/// ],
/// ..Default::default()
/// },
/// Faces::from_iter(&[
/// &[3, 2, 1, 0], &[0, 1, 5, 4], &[1, 2, 6, 5],
/// &[2, 3, 7, 6], &[3, 0, 4, 7], &[4, 5, 6, 7],
/// ]),
/// );
///
/// let submesh = mesh.create_mesh_by_face_indices(&[0, 3, 5]);
///
/// // the same attributes vector
/// assert_eq!(mesh.positions(), submesh.positions());
///
/// let faces0 = Faces::from_iter(&[
/// &[3, 2, 1, 0], &[2, 3, 7, 6], &[4, 5, 6, 7],
/// ]);
/// assert_eq!(submesh.faces().len(), faces0.len());
/// assert_eq!(submesh.faces()[0], faces0[0]);
/// assert_eq!(submesh.faces()[1], faces0[1]);
/// assert_eq!(submesh.faces()[2], faces0[2]);
/// ```
fn create_mesh_by_face_indices(&self, indices: &[usize]) -> PolygonMesh;
/// Extracts polygons such that there exists its normal is the same as its face normal.
/// # Arguments
/// - `tol`: tolerance to be regarded as the same normal as the face normal
/// # Returns
/// - The first polygon consists the faces included in planes.
/// - The second polygon is the extracted remainder.
/// # Panics
/// `tol` must be more than `TOLERANCE`.
/// # Examples
/// ```
/// use truck_polymesh::*;
/// use truck_meshalgo::analyzers::*;
/// let mesh = PolygonMesh::new(
/// StandardAttributes {
/// positions: vec![
/// Point3::new(0.0, 0.5, 0.0),
/// Point3::new(0.0, 0.5, 1.0),
/// Point3::new(1.0, 0.5, 1.0),
/// Point3::new(1.0, 0.5, 0.0),
/// Point3::new(0.0, 0.0, 2.0),
/// Point3::new(1.0, 0.0, 2.0),
/// ],
/// normals: vec![
/// Vector3::new(0.0, 1.0, 0.0),
/// // displaced normals for smooth rendering
/// Vector3::new(0.0, 1.0, 1.0).normalize(),
/// ],
/// ..Default::default()
/// },
/// Faces::from_iter(&[
/// &[(0, None, Some(0)), (1, None, Some(0)), (2, None, Some(0)), (3, None, Some(0))],
/// &[(2, None, Some(0)), (1, None, Some(0)), (4, None, Some(1)), (5, None, Some(1))],
/// ]),
/// );
/// let (plane, remained) = mesh.extract_planes(TOLERANCE); // TOLERANCE == 1.0e-7
/// assert_eq!(plane.len(), 1); assert_eq!(plane[0], 0);
/// assert_eq!(remained.len(), 1); assert_eq!(remained[0], 1);
/// ```
fn extract_planes(&self, tol: f64) -> (Vec<usize>, Vec<usize>);
/// Splits into the components.
/// # Details
/// Two polygons are considered to be in the same component if they share an edge
/// whose vertices has the same positions (and normals if `use_normal == true`).
/// # Examples
/// ```
/// use truck_polymesh::*;
/// use truck_meshalgo::{analyzers::*, filters::*};
///
/// // cube consisting tri_faces
/// let mut mesh = PolygonMesh::new(
/// StandardAttributes {
/// positions: vec![
/// Point3::new(0.0, 0.0, 0.0),
/// Point3::new(1.0, 0.0, 0.0),
/// Point3::new(1.0, 1.0, 0.0),
/// Point3::new(0.0, 1.0, 0.0),
/// Point3::new(0.0, 0.0, 1.0),
/// Point3::new(1.0, 0.0, 1.0),
/// Point3::new(1.0, 1.0, 1.0),
/// Point3::new(0.0, 1.0, 1.0),
/// ],
/// ..Default::default()
/// },
/// Faces::from_iter(&[
/// &[3, 2, 0], &[1, 0, 2], &[0, 1, 4], &[5, 4, 1],
/// &[1, 2, 5], &[6, 5, 2], &[2, 3, 6], &[7, 6, 3],
/// &[3, 0, 7], &[4, 7, 0], &[4, 5, 7], &[6, 7, 5],
/// ]),
/// );
///
/// // sign up normals
/// mesh.add_naive_normals(true).put_together_same_attrs(TOLERANCE);
///
/// // into components with normals
/// let components = mesh.components(true);
/// // The number of components is six because the mesh is a cube.
/// assert_eq!(components.len(), 6);
///
/// // into components without normals
/// let components = mesh.components(false);
/// // The number of components is one because the mesh is a connected cube.
/// assert_eq!(components.len(), 1);
/// ```
fn components(&self, use_normal: bool) -> Vec<Vec<usize>>;
}
impl Splitting for PolygonMesh {
fn create_mesh_by_face_indices(&self, indices: &[usize]) -> PolygonMesh {
let positions = self.positions().clone();
let uv_coords = self.uv_coords().clone();
let normals = self.normals().clone();
let faces: Faces = indices.iter().map(|i| &self.faces()[*i]).collect();
PolygonMesh::new(
StandardAttributes {
positions,
uv_coords,
normals,
},
faces,
)
}
fn extract_planes(&self, tol: f64) -> (Vec<usize>, Vec<usize>) {
nonpositive_tolerance!(tol, 0.0);
self.faces_into_two_clusters(|face: &[Vertex]| {
is_in_the_plane(self.positions(), self.normals(), face, tol * tol)
})
}
fn components(&self, use_normal: bool) -> Vec<Vec<usize>> {
let face_adjacency = self.faces().face_adjacency(use_normal);
get_components(&face_adjacency)
}
}
#[doc(hidden)]
pub trait ExperimentalSplitters {
fn faces_into_two_clusters<F: Fn(&[Vertex]) -> bool>(
&self,
func: F,
) -> (Vec<usize>, Vec<usize>);
fn clustering_faces_by_gcurvature(
&self,
threshold: f64,
preferred_upper: bool,
) -> (Vec<usize>, Vec<usize>);
fn get_gcurve(&self) -> Vec<f64>;
}
/// splitting the global faces
impl ExperimentalSplitters for PolygonMesh {
/// separate all faces into two clusters, one with `func` returning true and the other with false.
/// # Returns
/// (the vector of all faces `f` which `func(f) == true`, the one of the other faces)
fn faces_into_two_clusters<F: Fn(&[Vertex]) -> bool>(
&self,
func: F,
) -> (Vec<usize>, Vec<usize>) {
let mut true_faces = Vec::new();
let mut false_faces = Vec::new();
for (i, face) in self.face_iter().enumerate() {
match func(face) {
true => true_faces.push(i),
false => false_faces.push(i),
}
}
(true_faces, false_faces)
}
fn clustering_faces_by_gcurvature(
&self,
threshold: f64,
preferred_upper: bool,
) -> (Vec<usize>, Vec<usize>) {
let gcurve = self.get_gcurve();
self.faces_into_two_clusters(|face: &[Vertex]| {
is_signed_up_upper(face, &gcurve, preferred_upper, threshold)
})
}
fn get_gcurve(&self) -> Vec<f64> {
let positions = self.positions();
let mut angles = vec![0.0; positions.len()];
let mut weights = vec![0.0; positions.len()];
for face in self.tri_faces() {
angles[face[0].pos] += get_angle(positions, face, 0, 1, 2);
angles[face[1].pos] += get_angle(positions, face, 1, 2, 0);
angles[face[2].pos] += get_angle(positions, face, 2, 0, 1);
add_weights(&mut weights, positions, face);
}
for face in self.quad_faces() {
angles[face[0].pos] += get_angle(positions, face, 0, 1, 3);
angles[face[1].pos] += get_angle(positions, face, 1, 2, 0);
angles[face[2].pos] += get_angle(positions, face, 2, 3, 1);
angles[face[3].pos] += get_angle(positions, face, 3, 0, 1);
add_weights(&mut weights, positions, face);
}
for face in self.other_faces() {
let n = face.len() - 1;
angles[face[0].pos] += get_angle(positions, face, 0, 1, n);
for i in 1..n {
angles[face[i].pos] += get_angle(positions, face, i, i + 1, i - 1);
}
angles[face[n].pos] += get_angle(positions, face, n, 0, n - 1);
add_weights(&mut weights, positions, face);
}
angles
.into_iter()
.zip(weights)
.map(|(ang, weight)| (PI * 2.0 - ang) / weight)
.collect()
}
}
/// divide the graph to the connected components.
/// # Arguments
/// * adjacency - the adjacency matrix
/// # Return
/// * the list of the indices of faces contained in each components
fn get_components(adjacency: &[Vec<usize>]) -> Vec<Vec<usize>> {
let mut unchecked = vec![true; adjacency.len()];
let mut components = Vec::new();
loop {
let first = match unchecked.iter().position(|x| *x) {
Some(idx) => idx,
None => return components,
};
let mut stack = vec![first];
let mut component = vec![first];
unchecked[first] = false;
while let Some(cursor) = stack.pop() {
for i in &adjacency[cursor] {
if unchecked[*i] {
unchecked[*i] = false;
component.push(*i);
stack.push(*i);
}
}
}
components.push(component);
}
}
fn is_in_the_plane(positions: &[Point3], normals: &[Vector3], face: &[Vertex], tol2: f64) -> bool {
let n = FaceNormal::new(positions, face, 0).normal;
for v in face {
if let Some(nor) = v.nor {
if n.distance2(normals[nor]) < tol2 {
return true;
}
}
}
false
}
fn is_signed_up_upper(
face: &[Vertex],
gcurve: &[f64],
preferred_upper: bool,
threshold: f64,
) -> bool {
if preferred_upper {
face.as_ref().iter().any(|v| gcurve[v.pos] > threshold)
} else {
face.as_ref().iter().all(|v| gcurve[v.pos] > threshold)
}
}
fn get_angle(positions: &[Point3], face: &[Vertex], idx0: usize, idx1: usize, idx2: usize) -> f64 {
let vec0 = positions[face[idx1].pos] - positions[face[idx0].pos];
let vec1 = positions[face[idx2].pos] - positions[face[idx0].pos];
vec0.angle(vec1).0
}
fn add_weights(weights: &mut [f64], positions: &[Point3], face: &[Vertex]) {
let area = (2..face.len()).fold(0.0, |sum, i| {
let vec0 = positions[face[i - 1].pos] - positions[face[0].pos];
let vec1 = positions[face[i].pos] - positions[face[0].pos];
sum + (vec0.cross(vec1)).magnitude() / 2.0
}) / (face.len() as f64);
for v in face {
weights[v.pos] += area;
}
}
#[test]
fn into_components_test() {
let faces = vec![
[(0, None, Some(0)), (1, None, Some(1)), (2, None, Some(2))].as_ref(),
&[(0, None, Some(0)), (5, None, Some(5)), (1, None, Some(1))],
&[(0, None, Some(7)), (5, None, Some(5)), (6, None, Some(6))],
&[(5, None, Some(5)), (6, None, Some(6)), (7, None, Some(7))],
&[
(1, None, Some(1)),
(4, None, Some(4)),
(3, None, Some(3)),
(2, None, Some(2)),
],
]
.into_iter()
.collect();
let positions = vec![Point3::origin(); 8];
let normals = vec![Vector3::unit_x(); 8];
let mesh = PolygonMesh::new(
StandardAttributes {
positions,
normals,
..Default::default()
},
faces,
);
let comp = mesh.components(true);
assert_eq!(comp.len(), 2);
assert_eq!(comp[0], vec![0, 1, 4]);
assert_eq!(comp[1], vec![2, 3]);
}
+117
View File
@@ -0,0 +1,117 @@
use super::*;
use itertools::Itertools;
use rustc_hash::{FxHashMap as HashMap, FxHashSet as HashSet};
use truck_topology::shell::ShellCondition;
/// Extracts boundaries, and check shell condition: closed or orientation.
pub trait Topology {
/// Returns a vector of all boundaries as line strip.
fn extract_boundaries(&self) -> Vec<Vec<usize>>;
/// Determines the shell conditions: non-regular, regular, oriented, or closed.
/// The complexity increases in proportion to the number of edges.
///
/// Examples for each condition can be found on the page of
/// [`ShellCondition`](https://docs.rs/truck-topology/0.2.0/truck_topology/shell/enum.ShellCondition.html).
fn shell_condition(&self) -> ShellCondition;
}
#[derive(Clone, Debug)]
struct Boundaries {
checked: HashSet<[usize; 2]>,
boundary: HashMap<[usize; 2], bool>,
condition: ShellCondition,
}
impl Boundaries {
#[inline(always)]
fn new() -> Self {
Boundaries {
checked: Default::default(),
boundary: Default::default(),
condition: ShellCondition::Oriented,
}
}
#[inline(always)]
fn insert(&mut self, edge: [Vertex; 2]) {
let ori = edge[0].pos < edge[1].pos;
let edge = match ori {
true => [edge[0].pos, edge[1].pos],
false => [edge[1].pos, edge[0].pos],
};
self.condition = self.condition
& match (self.checked.insert(edge), self.boundary.insert(edge, ori)) {
(true, None) => ShellCondition::Oriented,
(false, None) => ShellCondition::Irregular,
(true, Some(_)) => panic!("unexpected case!"),
(false, Some(ori0)) => {
self.boundary.remove(&edge);
match ori == ori0 {
true => ShellCondition::Regular,
false => ShellCondition::Oriented,
}
}
};
}
#[inline(always)]
fn condition(&self) -> ShellCondition {
if self.condition == ShellCondition::Oriented && self.boundary.is_empty() {
ShellCondition::Closed
} else {
self.condition
}
}
}
impl FromIterator<[Vertex; 2]> for Boundaries {
fn from_iter<I: IntoIterator<Item = [Vertex; 2]>>(iter: I) -> Boundaries {
let mut boundaries = Boundaries::new();
iter.into_iter().for_each(|edge| boundaries.insert(edge));
boundaries
}
}
fn face_edge_iter<T: Copy>(face: &[T]) -> impl Iterator<Item = [T; 2]> + '_ {
face.iter().circular_tuple_windows().map(|(v, w)| [*v, *w])
}
impl Topology for Faces {
fn extract_boundaries(&self) -> Vec<Vec<usize>> {
let boundaries: Boundaries = self.face_iter().flat_map(face_edge_iter).collect();
let mut vemap: HashMap<usize, usize> = boundaries
.boundary
.into_iter()
.map(|(edge, ori)| match ori {
true => (edge[0], edge[1]),
false => (edge[1], edge[0]),
})
.collect();
let mut res = Vec::new();
while !vemap.is_empty() {
let mut wire = Vec::new();
let front = vemap.iter().next().unwrap();
let front = (*front.0, *front.1);
vemap.remove(&front.0);
wire.push(front.0);
let mut cursor = front.1;
while cursor != front.0 {
wire.push(cursor);
cursor = vemap.remove(&cursor).unwrap_or(front.0);
}
res.push(wire);
}
res
}
fn shell_condition(&self) -> ShellCondition {
self.face_iter()
.flat_map(face_edge_iter)
.collect::<Boundaries>()
.condition()
}
}
impl Topology for PolygonMesh {
fn extract_boundaries(&self) -> Vec<Vec<usize>> { self.faces().extract_boundaries() }
fn shell_condition(&self) -> ShellCondition { self.faces().shell_condition() }
}
+175
View File
@@ -0,0 +1,175 @@
use super::*;
use array_macro::array;
/// Calculate the volume and the center of gravity of mesh.
pub trait CalcVolume {
/// Returns the volume of the mesh if the mesh is closed.
///
/// # Details
/// More strictly, it returns the integral value of `xdydz`. Because of the linearity of this operation,
/// if multiple meshes represent a single closed geometry, the overall volume can be calculated by computing
/// this value for each and adding them together.
///
/// # Examples
/// ```
/// // The regular trihedron in the unit sphere.
/// use truck_meshalgo::prelude::*;
/// let a = 1.0 / f64::sqrt(3.0);
/// let positions = vec![
/// Point3::new(-a, -a, -a),
/// Point3::new(a, a, -a),
/// Point3::new(a, -a, a),
/// Point3::new(-a, a, a),
/// ];
/// let attrs = StandardAttributes {
/// positions,
/// ..Default::default()
/// };
/// let faces = Faces::from_iter(&[
/// [0, 1, 2],
/// [1, 3, 2],
/// [1, 0, 3],
/// [0, 2, 3],
/// ]);
/// let polygon = PolygonMesh::new(attrs, faces);
/// assert_near!(polygon.volume(), 8.0 * f64::sqrt(3.0) / 27.0);
/// ```
/// ```
/// // Cube with a hollow inside, defined by two boundaries.
/// use truck_meshalgo::prelude::*;
/// let positions0 = vec![
/// Point3::new(0.0, 0.0, 0.0),
/// Point3::new(3.0, 0.0, 0.0),
/// Point3::new(0.0, 3.0, 0.0),
/// Point3::new(0.0, 0.0, 3.0),
/// Point3::new(0.0, 3.0, 3.0),
/// Point3::new(3.0, 0.0, 3.0),
/// Point3::new(3.0, 3.0, 0.0),
/// Point3::new(3.0, 3.0, 3.0),
/// ];
/// let attrs0 = StandardAttributes {
/// positions: positions0,
/// ..Default::default()
/// };
/// let faces0 = Faces::from_iter(&[
/// [0, 2, 6, 1],
/// [0, 1, 5, 3],
/// [1, 6, 7, 5],
/// [6, 2, 4, 7],
/// [2, 0, 3, 4],
/// [3, 5, 7, 4],
/// ]);
/// let polygon0 = PolygonMesh::new(attrs0, faces0);
/// let positions1 = vec![
/// Point3::new(1.0, 1.0, 1.0),
/// Point3::new(2.0, 1.0, 1.0),
/// Point3::new(1.0, 2.0, 1.0),
/// Point3::new(1.0, 1.0, 2.0),
/// Point3::new(1.0, 2.0, 2.0),
/// Point3::new(2.0, 1.0, 2.0),
/// Point3::new(2.0, 2.0, 1.0),
/// Point3::new(2.0, 2.0, 2.0),
/// ];
/// let attrs1 = StandardAttributes {
/// positions: positions1,
/// ..Default::default()
/// };
/// let faces1 = Faces::from_iter(&[
/// [1, 6, 2, 0],
/// [3, 5, 1, 0],
/// [5, 7, 6, 1],
/// [7, 4, 2, 6],
/// [4, 3, 0, 2],
/// [4, 7, 5, 3],
/// ]);
/// let polygon1 = PolygonMesh::new(attrs1, faces1);
/// // We consider the solid defined by two boundaries(polygon0 and polygon1).
/// // The volume will be calculated by `position0.volume() + position1.volume()`.
/// assert_near!(polygon0.volume() + polygon1.volume(), 26.0);
/// ```
fn volume(&self) -> f64;
/// Returns the center of gravity of the mesh by homogeneous coordinate if the mesh is closed.
///
/// # Details
/// More strictly, it returns the integral value of `(1/2 x^2 dydz, 1/2 y^2dzdx, 1/2 z^2 dxdy, xdydz)`.
/// Because of the linearity of this operation, if multiple meshes represent a single closed geometry,
/// the overall volume can be calculated by computing this value for each and adding them together.
///
/// # Examples
/// ```
/// // The triangular prism
/// use truck_meshalgo::prelude::*;
/// let positions = vec![
/// Point3::new(0.0, 0.0, 0.0),
/// Point3::new(2.0, 5.0, 0.0),
/// Point3::new(-5.0, 1.0, 0.0),
/// Point3::new(0.0, 0.0, 2.0),
/// Point3::new(2.0, 5.0, 2.0),
/// Point3::new(-5.0, 1.0, 2.0),
/// ];
/// let attrs = StandardAttributes {
/// positions,
/// ..Default::default()
/// };
/// let faces = Faces::from_iter(&[
/// [2, 1, 0].as_slice(),
/// &[0, 1, 4, 3],
/// &[1, 2, 5, 4],
/// &[2, 0, 3, 5],
/// &[3, 4, 5],
/// ]);
/// let polygon = PolygonMesh::new(attrs, faces);
/// let homog = polygon.center_of_gravity();
/// assert_near!(homog.to_point(), Point3::new(-1.0, 2.0, 1.0));
/// ```
fn center_of_gravity(&self) -> Vector4;
}
impl CalcVolume for PolygonMesh {
fn volume(&self) -> f64 {
point_triangles(self).fold(0.0, |sum, [p, q, r]| {
sum + (p.x + q.x + r.x) * ((q.y - p.y) * (r.z - p.z) - (r.y - p.y) * (q.z - p.z))
}) / 6.0
}
fn center_of_gravity(&self) -> Vector4 {
let arr = point_triangles(self).fold([0.0; 4], |sum, [p, q, r]| {
let det = array![i => {
let (j, k) = ((i + 1) % 3, (i + 2) % 3);
(q[j] - p[j]) * (r[k] - p[k]) - (r[j] - p[j]) * (q[k] - p[k])
}; 3];
let vals = array![i => {
let s = p[i] + q[i] + r[i];
s * s - p[i] * q[i] - q[i] * r[i] - r[i] * p[i]
}; 3];
let res = array![i => sum[i] + vals[i] * det[i]; 3];
let res3 = sum[3] + (p.x + q.x + r.x) * det[0];
[res[0], res[1], res[2], res3]
});
Vector4::new(arr[0] / 24.0, arr[1] / 24.0, arr[2] / 24.0, arr[3] / 6.0)
}
}
fn point_triangles(poly: &PolygonMesh) -> impl Iterator<Item = [Point3; 3]> + '_ {
poly.faces()
.triangle_iter()
.map(|faces| array![i => poly.positions()[faces[i].pos]; 3])
}
impl CalcVolume for truck_topology::Solid<Point3, PolylineCurve<Point3>, PolygonMesh> {
fn volume(&self) -> f64 {
self.face_iter()
.map(|face| match face.orientation() {
true => face.surface().volume(),
false => -face.surface().volume(),
})
.sum::<f64>()
}
fn center_of_gravity(&self) -> Vector4 {
self.face_iter()
.map(|face| match face.orientation() {
true => face.surface().center_of_gravity(),
false => -face.surface().center_of_gravity(),
})
.sum::<Vector4>()
}
}
+77
View File
@@ -0,0 +1,77 @@
use super::*;
use rustc_hash::{FxHashMap as HashMap, FxHashSet as HashSet};
pub trait Adjacency {
#[allow(dead_code)]
fn vertex_adjacency(&self, num_of_vertices: usize) -> Vec<Vec<usize>>;
/// create the adjacency list of the faces
fn face_adjacency(&self, use_normal: bool) -> Vec<Vec<usize>>;
}
impl Adjacency for Faces {
fn vertex_adjacency(&self, num_of_vertices: usize) -> Vec<Vec<usize>> {
let mut already = HashSet::default();
let mut res = vec![Vec::new(); num_of_vertices];
for face in self.face_iter() {
face.windows(2)
.chain(std::iter::once([face[face.len() - 1], face[0]].as_ref()))
.map(|edge| [edge[0].pos, edge[1].pos])
.for_each(|edge| {
let first = match edge[0] < edge[1] {
true => already.insert((edge[0], edge[1])),
false => already.insert((edge[1], edge[0])),
};
if first {
res[edge[0]].push(edge[1]);
res[edge[1]].push(edge[0]);
}
})
}
res
}
fn face_adjacency(&self, use_normal: bool) -> Vec<Vec<usize>> {
let len = self.len();
let mut face_adjacency = vec![Vec::<usize>::new(); len];
let mut edge_face_map = HashMap::default();
for (i, face) in self.face_iter().enumerate() {
face.windows(2)
.chain(std::iter::once([face[face.len() - 1], face[0]].as_ref()))
.for_each(|v| {
signup_adjacency(
i,
v[0],
v[1],
&mut face_adjacency,
&mut edge_face_map,
use_normal,
)
})
}
face_adjacency
}
}
fn signup_adjacency(
i: usize,
v0: Vertex,
v1: Vertex,
face_adjacency: &mut [Vec<usize>],
edge_face_map: &mut HashMap<[(usize, Option<usize>); 2], usize>,
use_normal: bool,
) {
let edge = match (v0.pos < v1.pos, use_normal) {
(true, true) => [(v0.pos, v0.nor), (v1.pos, v1.nor)],
(false, true) => [(v1.pos, v1.nor), (v0.pos, v0.nor)],
(true, false) => [(v0.pos, None), (v1.pos, None)],
(false, false) => [(v1.pos, None), (v0.pos, None)],
};
match edge_face_map.get(&edge) {
Some(j) => {
face_adjacency[i].push(*j);
face_adjacency[*j].push(i);
}
None => {
edge_face_map.insert(edge, i);
}
}
}
+26
View File
@@ -0,0 +1,26 @@
use super::*;
#[derive(Clone, Copy, Debug)]
pub struct FaceNormal {
pub face_id: usize,
pub normal: Vector3,
}
impl FaceNormal {
pub fn new(positions: &[Point3], face: &[Vertex], face_id: usize) -> FaceNormal {
let center = face
.iter()
.fold(Vector3::zero(), |sum, v| sum + positions[v.pos].to_vec())
/ face.len() as f64;
let normal = face
.windows(2)
.chain(std::iter::once([face[face.len() - 1], face[0]].as_ref()))
.fold(Vector3::zero(), |sum, v| {
let vec0 = positions[v[0].pos].to_vec() - center;
let vec1 = positions[v[1].pos].to_vec() - center;
sum + vec0.cross(vec1)
})
.normalize();
FaceNormal { face_id, normal }
}
}
+6
View File
@@ -0,0 +1,6 @@
use crate::*;
mod adjacency;
mod face_normal;
pub(super) use adjacency::Adjacency;
pub(super) use face_normal::FaceNormal;
+105
View File
@@ -0,0 +1,105 @@
use super::*;
/// STL face generate from `PolygonMesh`
#[derive(Debug)]
pub struct TriangleIter<'a> {
tri_faces: std::slice::Iter<'a, [Vertex; 3]>,
quad_faces: std::slice::Iter<'a, [Vertex; 4]>,
other_faces: std::slice::Iter<'a, [Vertex; 3]>,
current_face: Option<&'a [Vertex]>,
current_vertex: usize,
len: usize,
}
impl<'a> Iterator for TriangleIter<'a> {
type Item = [Vertex; 3];
fn next(&mut self) -> Option<[Vertex; 3]> {
if self.current_face.is_none() {
self.current_face = if let Some(face) = self.tri_faces.next() {
Some(face)
} else if let Some(face) = self.quad_faces.next() {
Some(face)
} else if let Some(face) = self.other_faces.next() {
Some(face)
} else {
return None;
}
}
let face = self.current_face.unwrap();
let res = [
face[0],
face[self.current_vertex + 1],
face[self.current_vertex + 2],
];
if self.current_vertex + 3 == face.len() {
self.current_face = None;
self.current_vertex = 0;
} else {
self.current_vertex += 1;
}
self.len -= 1;
Some(res)
}
#[inline(always)]
fn size_hint(&self) -> (usize, Option<usize>) { (self.len, Some(self.len)) }
}
impl<'a> ExactSizeIterator for TriangleIter<'a> {}
#[derive(Clone, Debug)]
pub struct Triangulate<'a> {
entity: &'a PolygonMesh,
other_faces: Vec<[Vertex; 3]>,
}
impl<'a> Triangulate<'a> {
#[inline(always)]
pub fn new(entity: &'a PolygonMesh) -> Self {
fn divide_mesh<T: Copy>(face: &Vec<T>) -> impl Iterator<Item = [T; 3]> + '_ {
face.windows(2).skip(1).map(move |a| [face[0], a[0], a[1]])
}
Triangulate {
other_faces: entity.other_faces().iter().flat_map(divide_mesh).collect(),
entity,
}
}
#[inline(always)]
pub fn get(&self, idx: usize) -> [Vertex; 3] {
let tri_and_quad = self.entity.tri_faces().len() + 2 * self.entity.quad_faces().len();
if idx < self.entity.tri_faces().len() {
self.entity.tri_faces()[idx]
} else if idx < tri_and_quad {
let idx = idx - self.entity.tri_faces().len();
let face = self.entity.quad_faces()[idx / 2];
if idx % 2 == 0 {
[face[0], face[1], face[2]]
} else {
[face[0], face[2], face[3]]
}
} else {
self.other_faces[idx - tri_and_quad]
}
}
#[inline(always)]
pub fn entity(&self) -> &PolygonMesh { self.entity }
}
impl<'a> IntoIterator for &'a Triangulate<'a> {
type Item = [Vertex; 3];
type IntoIter = TriangleIter<'a>;
#[inline(always)]
fn into_iter(self) -> TriangleIter<'a> {
let len = self
.entity
.face_iter()
.fold(0, |len, face| len + face.len());
TriangleIter {
tri_faces: self.entity.tri_faces().iter(),
quad_faces: self.entity.quad_faces().iter(),
other_faces: self.other_faces.iter(),
current_face: None,
current_vertex: 0,
len,
}
}
}
+12
View File
@@ -0,0 +1,12 @@
use crate::*;
use polygon_mesh::PolygonMeshEditor;
mod normal_filters;
mod optimizing;
mod structuring;
mod subdivision;
pub use normal_filters::NormalFilters;
pub use optimizing::OptimizingFilter;
pub use structuring::StructuringFilter;
pub use subdivision::Subdivision;
+396
View File
@@ -0,0 +1,396 @@
use super::*;
use rustc_hash::FxHashMap as HashMap;
/// Filters for adding normals
pub trait NormalFilters {
/// Normalize all normals and assign `None` to the `nor` index of the vertices
/// that has irregular normals.
/// # Examples
/// ```
/// use truck_polymesh::*;
/// use truck_meshalgo::filters::*;
///
/// // Morbid data only for testing
/// let mut mesh = PolygonMesh::new(
/// StandardAttributes {
/// positions: vec![Point3::new(0.0, 0.0, 0.0)],
/// normals: vec![
/// Vector3::new(100.0, 20.0, 56.0),
/// Vector3::new(1.0e-12, 3.536e10, std::f64::NAN),
/// Vector3::new(0.0, 1.0, 0.0),
/// ],
/// ..Default::default()
/// },
/// Faces::from_iter(&[&[
/// (0, None, Some(0)),
/// (0, None, Some(1)),
/// (0, None, Some(2))
/// ]]),
/// );
///
/// mesh.normalize_normals();
/// assert!(mesh.normals()[0].magnitude().near(&1.0));
/// assert_eq!(mesh.faces()[0][1].nor, None);
/// ```
fn normalize_normals(&mut self) -> &mut Self;
/// Adds face normals to each vertices.
/// # Arguments
/// - If `overwrite == true`, clear all normals and update all normals in vertices.
/// - If `overwrite == false`, add normals only for `nor` is `None`.
/// # Examples
/// Compare with the examples of [`add_smooth_normals`](./trait.NormalFilters.html#tymethod.add_smooth_normals).
/// ```
/// use truck_meshalgo::prelude::*;
/// let mut mesh = PolygonMesh::new(
/// StandardAttributes {
/// positions: vec![
/// Point3::new(-5.0, 0.0, 0.0),
/// Point3::new(0.0, 2.0, -2.0),
/// Point3::new(0.0, 2.0, 0.0),
/// Point3::new(0.0, 2.0, 2.0),
/// Point3::new(5.0, 0.0, 0.0),
/// ],
/// ..Default::default()
/// },
/// Faces::from_iter(&[
/// &[0, 2, 1],
/// &[0, 3, 2],
/// &[1, 2, 4],
/// &[2, 3, 4],
/// ]),
/// );
///
/// mesh.add_naive_normals(true);
/// let v0: StandardVertex = mesh.faces()[0][1];
/// let v1: StandardVertex = mesh.faces()[3][0];
///
/// // those vertices are at position with the index 2.
/// assert_eq!(v0.pos, 2); assert_eq!(v1.pos, 2);
///
/// // Each normal is each face normal.
/// assert!(mesh.normals()[v0.nor.unwrap()].near(&Vector3::new(-2.0, 5.0, 0.0).normalize()));
/// assert!(mesh.normals()[v1.nor.unwrap()].near(&Vector3::new(2.0, 5.0, 0.0).normalize()));
/// ```
fn add_naive_normals(&mut self, overwrite: bool) -> &mut Self;
/// add the smooth normal vectors to the mesh.
/// # Details
/// For each vertex, apply the following algorithm:
/// 1. prepare vectors that enumerate the normals of the faces containing the target vertices in order.
/// 1. cluster each normal `n` in turn in the following manner.
/// * If there is an existing cluster `A` in which the angle between the weighted
/// average of `A` and `n` is less than or equal to `tol_ang`, add `n` to `A`.
/// * If cluster `A` as described above does not exist, create a new cluster that contains only `n`.
/// # Arguments
/// - If `overwrite == true`, clear all normals and update all normals in vertices.
/// - If `overwrite == false`, add normals only for `nor` is `None`.
/// # Examples
/// Compare with the examples of [`add_smooth_normals`](./trait.NormalFilters.html#tymethod.add_smooth_normals).
/// ```
/// use truck_meshalgo::prelude::*;
/// let mut mesh = PolygonMesh::new(
/// StandardAttributes {
/// positions: vec![
/// Point3::new(-5.0, 0.0, 0.0),
/// Point3::new(0.0, 2.0, -2.0),
/// Point3::new(0.0, 2.0, 0.0),
/// Point3::new(0.0, 2.0, 2.0),
/// Point3::new(5.0, 0.0, 0.0),
/// ],
/// ..Default::default()
/// },
/// Faces::from_iter(&[
/// &[0, 2, 1], &[0, 3, 2], &[1, 2, 4], &[2, 3, 4],
/// ]),
/// );
///
/// mesh.add_smooth_normals(0.8, true);
/// let v0: StandardVertex = mesh.faces()[0][1];
/// let v1: StandardVertex = mesh.faces()[3][0];
///
/// // those vertices are at position with the index 2.
/// assert_eq!(v0.pos, 2); assert_eq!(v1.pos, 2);
///
/// // Normals are avaraged!
/// assert!(mesh.normals()[v0.nor.unwrap()].near(&Vector3::new(0.0, 1.0, 0.0)));
/// assert!(mesh.normals()[v1.nor.unwrap()].near(&Vector3::new(0.0, 1.0, 0.0)));
/// assert_eq!(v0.nor, v1.nor);
///
/// // If the tolerance is enough little, the faces are recognized as edges.
/// mesh.add_smooth_normals(0.6, true); // Normals are overwritten!
/// let v0: StandardVertex = mesh.faces()[0][1];
/// let v1: StandardVertex = mesh.faces()[3][0];
/// assert!(mesh.normals()[v0.nor.unwrap()].near(&Vector3::new(-2.0, 5.0, 0.0).normalize()));
/// assert!(mesh.normals()[v1.nor.unwrap()].near(&Vector3::new(2.0, 5.0, 0.0).normalize()));
/// ```
fn add_smooth_normals(&mut self, tol_ang: f64, overwrite: bool) -> &mut Self;
/// Makes the orientation of faces compatible to the normal vectors.
/// # Examples
/// ```
/// use truck_polymesh::*;
/// use truck_meshalgo::filters::*;
/// let mut polymesh = PolygonMesh::new(
/// StandardAttributes {
/// positions: vec![
/// Point3::new(0.0, 0.0, 0.0),
/// Point3::new(1.0, 0.0, 0.0),
/// Point3::new(0.0, 1.0, 0.0),
/// ],
/// normals: vec![Vector3::new(0.0, 0.0, -1.0)],
/// ..Default::default()
/// },
/// Faces::from_tri_and_quad_faces(
/// vec![[
/// (0, None, Some(0)).into(),
/// (1, None, Some(0)).into(),
/// (2, None, Some(0)).into(),
/// ]],
/// Vec::new(),
/// ),
/// );
/// polymesh.make_face_compatible_to_normal();
/// assert_eq!(polymesh.faces()[0][0].pos, 2);
/// assert_eq!(polymesh.faces()[0][1].pos, 1);
/// assert_eq!(polymesh.faces()[0][2].pos, 0);
/// ```
fn make_face_compatible_to_normal(&mut self) -> &mut Self;
/// Makes the orientation of faces compatible to the normal vectors.
/// # Examples
/// ```
/// use truck_polymesh::*;
/// use truck_meshalgo::filters::*;
/// let mut polymesh = PolygonMesh::new(
/// StandardAttributes {
/// positions: vec![
/// Point3::new(0.0, 0.0, 0.0),
/// Point3::new(1.0, 0.0, 0.0),
/// Point3::new(0.0, 1.0, 0.0),
/// ],
/// normals: vec![Vector3::new(0.0, 0.0, -1.0)],
/// ..Default::default()
/// },
/// Faces::from_tri_and_quad_faces(
/// vec![[
/// (0, None, Some(0)).into(),
/// (1, None, Some(0)).into(),
/// (2, None, Some(0)).into(),
/// ]],
/// Vec::new(),
/// ),
/// );
/// polymesh.make_normal_compatible_to_face();
/// let nor = polymesh.faces()[0][0].nor.unwrap();
/// assert_eq!(polymesh.normals()[nor], Vector3::new(0.0, 0.0, 1.0));
/// ```
fn make_normal_compatible_to_face(&mut self) -> &mut Self;
}
impl NormalFilters for PolygonMesh {
fn normalize_normals(&mut self) -> &mut Self {
let mut mesh = self.debug_editor();
let PolygonMeshEditor {
attributes: StandardAttributes { normals, .. },
faces,
..
} = &mut mesh;
normals
.iter_mut()
.for_each(move |normal| *normal = normal.normalize());
faces.face_iter_mut().flatten().for_each(|v| {
if let Some(idx) = v.nor {
if !normals[idx].magnitude2().near(&1.0) {
v.nor = None;
}
}
});
drop(mesh);
self
}
fn make_face_compatible_to_normal(&mut self) -> &mut Self {
let mut mesh = self.debug_editor();
let PolygonMeshEditor {
attributes: StandardAttributes {
positions, normals, ..
},
faces,
..
} = &mut mesh;
for face in faces.face_iter_mut() {
let normal = face.iter().fold(Vector3::zero(), |normal, v| {
normal + v.nor.map(|i| normals[i]).unwrap_or_else(Vector3::zero)
});
let face_normal = FaceNormal::new(positions, face, 0).normal;
if normal.dot(face_normal) < 0.0 {
face.reverse();
}
}
drop(mesh);
self
}
fn make_normal_compatible_to_face(&mut self) -> &mut Self {
let mut mesh = self.debug_editor();
let PolygonMeshEditor {
attributes: StandardAttributes {
positions, normals, ..
},
faces,
..
} = &mut mesh;
for face in faces.face_iter_mut() {
let face_normal = FaceNormal::new(positions, face, 0).normal;
face.iter_mut().for_each(|v| {
v.nor = v.nor.map(|idx| {
if normals[idx].dot(face_normal) < 0.0 {
normals.push(-normals[idx]);
normals.len() - 1
} else {
idx
}
});
})
}
drop(mesh);
self
}
fn add_naive_normals(&mut self, overwrite: bool) -> &mut Self {
let mut mesh = self.debug_editor();
let PolygonMeshEditor {
attributes: StandardAttributes {
positions, normals, ..
},
faces,
..
} = &mut mesh;
if overwrite {
normals.clear()
}
faces.face_iter_mut().for_each(move |face| {
let normal = FaceNormal::new(positions, face, 0).normal;
let mut added = false;
face.iter_mut().for_each(|v| {
if v.nor.is_none() || overwrite {
if !added {
normals.push(normal);
added = true;
}
v.nor = Some(normals.len() - 1);
}
});
});
drop(mesh);
self
}
fn add_smooth_normals(&mut self, tol_ang: f64, overwrite: bool) -> &mut Self {
let vnmap = self.clustering_noraml_faces(tol_ang.cos());
self.reflect_normal_clusters(vnmap, overwrite);
self
}
}
trait SubNormalFilter {
fn clustering_noraml_faces(&self, inf: f64) -> HashMap<usize, Vec<Vec<FaceNormal>>>;
fn reflect_normal_clusters(
&mut self,
vnmap: HashMap<usize, Vec<Vec<FaceNormal>>>,
overwrite: bool,
);
}
impl SubNormalFilter for PolygonMesh {
fn clustering_noraml_faces(&self, inf: f64) -> HashMap<usize, Vec<Vec<FaceNormal>>> {
let positions = self.positions();
let mut vnmap = HashMap::default();
self.face_iter()
.enumerate()
.for_each(|(i, face)| add_face_normal(positions, i, face, &mut vnmap, inf));
vnmap
}
fn reflect_normal_clusters(
&mut self,
vnmap: HashMap<usize, Vec<Vec<FaceNormal>>>,
overwrite: bool,
) {
let mut mesh = self.debug_editor();
let PolygonMeshEditor {
attributes: StandardAttributes { normals, .. },
faces,
..
} = &mut mesh;
if overwrite {
normals.clear();
}
for (pos_id, vecs) in vnmap.into_iter() {
for vec in vecs {
let normal = vec
.iter()
.fold(Vector3::zero(), |sum, x| sum + x.normal)
.normalize();
for FaceNormal { face_id, .. } in vec {
signup_vertex_normal(pos_id, face_id, normals, normal, faces, overwrite);
}
}
}
}
}
fn add_face_normal(
positions: &[Point3],
face_id: usize,
face: &[Vertex],
vnmap: &mut HashMap<usize, Vec<Vec<FaceNormal>>>,
inf: f64,
) {
let face_normal = FaceNormal::new(positions, face, face_id);
face.iter().for_each(|v| {
add_to_vnmap(v.pos, face_normal, vnmap, inf);
})
}
fn add_to_vnmap(
pos_id: usize,
face_normal: FaceNormal,
vnmap: &mut HashMap<usize, Vec<Vec<FaceNormal>>>,
inf: f64,
) {
match vnmap.get_mut(&pos_id) {
Some(vecs) => {
for vec in vecs.iter_mut() {
let normal = vec
.iter()
.fold(Vector3::zero(), |sum, x| sum + x.normal)
.normalize();
if face_normal.normal.dot(normal) > inf {
vec.push(face_normal);
return;
}
}
vecs.push(vec![face_normal]);
}
None => {
let vecs = vec![vec![face_normal]];
vnmap.insert(pos_id, vecs);
}
}
}
fn signup_vertex_normal(
pos_id: usize,
face_id: usize,
normals: &mut Vec<Vector3>,
normal: Vector3,
faces: &mut Faces,
overwrite: bool,
) {
let face = faces[face_id].as_mut();
let j = (0..face.len()).find(|j| face[*j].pos == pos_id).unwrap();
if face[j].nor.is_none() || overwrite {
if let Some(n) = normals.last() {
if n != &normal {
normals.push(normal);
}
} else {
normals.push(normal);
}
face[j].nor = Some(normals.len() - 1);
}
}
+441
View File
@@ -0,0 +1,441 @@
use super::*;
use rustc_hash::FxHashMap as HashMap;
use std::iter::Iterator;
use std::ops::{Div, Mul};
/// Filters for optimizing data
pub trait OptimizingFilter {
/// remove all unused position, texture coordinates, and normal vectors.
/// # Examples
/// ```
/// use truck_polymesh::*;
/// use truck_meshalgo::filters::*;
/// let mut mesh = PolygonMesh::new(
/// StandardAttributes {
/// positions: vec![
/// Point3::new(0.0, 0.0, 0.0),
/// Point3::new(1.0, 0.0, 0.0),
/// Point3::new(0.0, 1.0, 0.0),
/// Point3::new(0.0, 0.0, 1.0),
/// ],
/// ..Default::default()
/// },
/// // 0 is not used!
/// Faces::from_iter(&[&[1, 2, 3]]),
/// );
///
/// assert_eq!(mesh.positions().len(), 4);
/// mesh.remove_unused_attrs();
/// assert_eq!(mesh.positions().len(), 3);
/// ```
fn remove_unused_attrs(&mut self) -> &mut Self;
/// Removes degenerate polygons.
/// # Examples
/// ```
/// use truck_polymesh::*;
/// use truck_meshalgo::filters::*;
/// let mut mesh = PolygonMesh::new(
/// StandardAttributes {
/// positions: vec![
/// Point3::new(0.0, 0.0, 0.0),
/// Point3::new(1.0, 0.0, 0.0),
/// Point3::new(0.0, 1.0, 0.0),
/// Point3::new(0.0, 0.0, 1.0),
/// ],
/// ..Default::default()
/// },
/// Faces::from_iter(&[
/// &[0, 1, 2],
/// &[2, 1, 2], // degenerate face!
/// &[2, 1, 3],
/// ]),
/// );
///
/// assert_eq!(mesh.faces().len(), 3);
/// mesh.remove_degenerate_faces();
/// assert_eq!(mesh.faces().len(), 2);
/// ```
fn remove_degenerate_faces(&mut self) -> &mut Self;
/// Gives the same indices to the same positions, texture coordinate, and normal vectors, respectively.
/// # Remarks
/// No longer needed attributes are NOT autoremoved.
/// One can remove such attributes by running [`remove_unused_attrs`] manually.
///
/// [`remove_unused_attrs`]: ./trait.WasteEliminatingFilter.html#tymethod.remove_unused_attrs
///
/// # Examples
/// ```
/// use truck_meshalgo::prelude::*;
/// let mut mesh = PolygonMesh::new(
/// StandardAttributes {
/// positions: vec![
/// Point3::new(0.0, 0.0, 0.0),
/// Point3::new(1.0, 0.0, 0.0),
/// Point3::new(0.0, 1.0, 0.0),
/// Point3::new(1.0, 1.0, 0.0),
/// Point3::new(0.0, 1.0, 0.0),
/// Point3::new(1.0, 0.0, 0.0),
/// ],
/// ..Default::default()
/// },
/// Faces::from_iter(&[
/// &[0, 1, 2],
/// &[3, 4, 5],
/// ]),
/// );
///
/// assert_eq!(mesh.faces()[1][1], StandardVertex { pos: 4, uv: None, nor: None });
/// mesh.put_together_same_attrs(TOLERANCE);
/// assert_eq!(mesh.faces()[1][1], StandardVertex { pos: 2, uv: None, nor: None });
///
/// // Remarks: No longer needed attributes are NOT autoremoved!
/// assert_eq!(mesh.positions().len(), 6);
/// mesh.remove_unused_attrs();
/// assert_eq!(mesh.positions().len(), 4);
/// ```
fn put_together_same_attrs(&mut self, tol: f64) -> &mut Self;
}
fn all_pos_mut(faces: &mut Faces) -> impl Iterator<Item = &mut usize> {
faces.face_iter_mut().flatten().map(move |v| &mut v.pos)
}
fn all_uv_mut(faces: &mut Faces) -> impl Iterator<Item = &mut usize> {
faces
.face_iter_mut()
.flatten()
.filter_map(move |v| v.uv.as_mut())
}
fn all_nor_mut(faces: &mut Faces) -> impl Iterator<Item = &mut usize> {
faces
.face_iter_mut()
.flatten()
.filter_map(move |v| v.nor.as_mut())
}
impl OptimizingFilter for PolygonMesh {
fn remove_unused_attrs(&mut self) -> &mut Self {
let mut mesh = self.debug_editor();
let PolygonMeshEditor {
attributes:
StandardAttributes {
positions,
uv_coords,
normals,
},
faces,
..
} = &mut mesh;
let pos_iter = all_pos_mut(faces);
let idcs = sub_remove_unused_attrs(pos_iter, positions.len());
*positions = idcs.iter().map(|i| positions[*i]).collect();
let uv_iter = all_uv_mut(faces);
let idcs = sub_remove_unused_attrs(uv_iter, uv_coords.len());
*uv_coords = idcs.iter().map(|i| uv_coords[*i]).collect();
let nor_iter = all_nor_mut(faces);
let idcs = sub_remove_unused_attrs(nor_iter, normals.len());
*normals = idcs.iter().map(|i| normals[*i]).collect();
drop(mesh);
self
}
fn remove_degenerate_faces(&mut self) -> &mut Self {
let mesh = self.debug_editor();
let mut faces = Faces::default();
for tri in mesh.faces.tri_faces() {
if !degenerate_triangle(*tri) {
faces.push(tri);
}
}
for quad in mesh.faces.quad_faces() {
match degenerate_quadrangle(*quad) {
QuadrangleType::TotallyDegenerate => {}
QuadrangleType::Triangle(tri) => faces.push(tri),
QuadrangleType::NonDegenerate => faces.push(quad),
}
}
for face in mesh.faces.other_faces() {
faces.extend(split_into_nondegenerate(face.clone()));
}
*mesh.faces = faces;
drop(mesh);
self
}
fn put_together_same_attrs(&mut self, tol: f64) -> &mut Self {
let mut mesh = self.debug_editor();
let PolygonMeshEditor {
attributes:
StandardAttributes {
positions,
uv_coords,
normals,
},
faces,
..
} = &mut mesh;
let bnd_box: BoundingBox<_> = positions.iter().collect();
let center = bnd_box.center();
let diag = bnd_box.diagonal().map(|a| f64::max(a.abs(), 1.0));
let normalized_positions = positions
.iter()
.map(move |position| 2.0 * (position - center).zip(diag, |a, b| a / b))
.collect::<Vec<_>>();
let pos_map = sub_put_together_same_attrs(&normalized_positions, tol);
all_pos_mut(faces).for_each(|idx| *idx = pos_map[*idx]);
let uv_map = sub_put_together_same_attrs(uv_coords, tol);
all_uv_mut(faces).for_each(|idx| *idx = uv_map[*idx]);
let nor_map = sub_put_together_same_attrs(normals, tol);
all_nor_mut(faces).for_each(|idx| *idx = nor_map[*idx]);
drop(mesh);
self
}
}
fn sub_remove_unused_attrs<'a, I: Iterator<Item = &'a mut usize>>(
iter: I,
old_len: usize,
) -> Vec<usize> {
let mut new2old = Vec::new();
let mut old2new = vec![None; old_len];
for idx in iter {
*idx = match old2new[*idx] {
Some(k) => k,
None => {
let k = new2old.len();
new2old.push(*idx);
old2new[*idx] = Some(k);
k
}
};
}
new2old
}
trait SameAttr: Copy + CastIntVector + MetricSpace<Metric = f64> {}
impl<T> SameAttr for T where T: Copy + CastIntVector + MetricSpace<Metric = f64> {}
fn sub_put_together_same_attrs<T: SameAttr>(attrs: &[T], tol: f64) -> Vec<usize> {
let map = create_blocks(attrs, tol);
let adj = create_adjacency(map, attrs, tol);
let components = create_components(adj, attrs);
centralize(components, attrs)
}
type BlockMap<'a, T> = HashMap<<T as CastIntVector>::IntVector, Vec<usize>>;
fn create_blocks<T: SameAttr>(attrs: &[T], tol: f64) -> BlockMap<'_, T> {
let mut map = HashMap::default();
attrs.iter().enumerate().for_each(|(i, attr)| {
attr.round(tol).neighborhood().into_iter().for_each(|idx| {
map.entry(idx).or_insert_with(Vec::new).push(i);
});
});
map
}
fn create_adjacency<T: SameAttr>(map: BlockMap<'_, T>, attrs: &[T], tol: f64) -> Vec<Vec<usize>> {
let mut adj = vec![Vec::new(); attrs.len()];
map.into_iter().for_each(|(_, vec)| {
if vec.len() > 1 {
vec.iter().copied().enumerate().for_each(|(k, i)| {
vec[(k + 1)..].iter().copied().for_each(|j| {
if attrs[i].distance2(attrs[j]) < tol * tol {
adj[i].push(j);
adj[j].push(i);
}
});
});
}
});
adj
}
fn create_components<T: SameAttr>(
mut adj: Vec<Vec<usize>>,
attrs: &[T],
) -> impl Iterator<Item = Vec<usize>> {
let mut already = vec![false; attrs.len()];
(0..adj.len()).filter_map(move |i| {
if already[i] {
return None;
}
let mut stack = vec![i];
let mut component = Vec::new();
while let Some(i) = stack.pop() {
if !already[i] {
component.push(i);
already[i] = true;
stack.append(&mut adj[i]);
}
}
Some(component)
})
}
fn centralize<T: SameAttr>(
components: impl Iterator<Item = Vec<usize>>,
attrs: &[T],
) -> Vec<usize> {
let mut res = vec![0; attrs.len()];
components.for_each(|component| {
let center = component
.iter()
.fold(T::zero(), |sum, j| sum.add_element_wise(attrs[*j]))
/ component.len() as f64;
let (mut mindist, mut mindx) = (f64::INFINITY, 0);
component.iter().copied().for_each(|i| {
let dist = attrs[i].distance2(center);
if dist < mindist {
mindist = dist;
mindx = i;
}
});
component.into_iter().for_each(|i| {
res[i] = mindx;
});
});
res
}
fn degenerate_triangle(tri: [Vertex; 3]) -> bool {
tri[0].pos == tri[1].pos || tri[1].pos == tri[2].pos || tri[2].pos == tri[0].pos
}
enum QuadrangleType {
NonDegenerate,
Triangle([Vertex; 3]),
TotallyDegenerate,
}
fn degenerate_quadrangle(quad: [Vertex; 4]) -> QuadrangleType {
if (quad[0].pos == quad[2].pos || quad[1].pos == quad[3].pos)
|| (quad[0].pos == quad[1].pos && quad[2].pos == quad[3].pos)
|| quad[1].pos == quad[2].pos && quad[3].pos == quad[0].pos
{
QuadrangleType::TotallyDegenerate
} else if quad[0].pos == quad[1].pos || quad[1].pos == quad[2].pos {
QuadrangleType::Triangle([quad[0], quad[2], quad[3]])
} else if quad[2].pos == quad[3].pos || quad[3].pos == quad[0].pos {
QuadrangleType::Triangle([quad[0], quad[1], quad[2]])
} else {
QuadrangleType::NonDegenerate
}
}
fn split_into_nondegenerate(poly: Vec<Vertex>) -> Vec<Vec<Vertex>> {
for i in 0..poly.len() {
for j in (i + 1)..poly.len() {
if poly[i].pos == poly[j].pos {
let polygon0: Vec<_> = (0..(j - i)).map(|k| poly[k + i]).collect();
let polygon1: Vec<_> = ((j - i)..poly.len())
.map(|k| poly[(k + i) % poly.len()])
.collect();
let mut result = split_into_nondegenerate(polygon0);
result.extend(split_into_nondegenerate(polygon1));
return result;
}
}
}
vec![poly]
}
trait CastIntVector:
Sized + ElementWise<Self> + Mul<f64, Output = Self> + Div<f64, Output = Self> {
type IntVector: IntVector;
fn cast_int(self) -> Self::IntVector;
fn round(self, tol: f64) -> Self::IntVector;
fn zero() -> Self;
}
macro_rules! impl_cast_int {
($typename: ident, $n: expr) => {
impl CastIntVector for $typename {
type IntVector = [i64; $n];
fn cast_int(self) -> [i64; $n] {
self.cast::<i64>()
.unwrap_or_else(|| panic!("failed to cast: {self:?}"))
.into()
}
fn round(self, tol: f64) -> Self::IntVector { (self / tol).cast_int() }
fn zero() -> Self { Self::from([0.0; $n]) }
}
};
}
impl_cast_int!(Vector2, 2);
impl_cast_int!(Vector3, 3);
impl_cast_int!(Point3, 3);
trait IntVector: Copy + std::hash::Hash + Eq {
type Neighborhood: IntoIterator<Item = Self>;
fn neighborhood(self) -> Self::Neighborhood;
}
impl IntVector for [i64; 2] {
type Neighborhood = [[i64; 2]; 9];
fn neighborhood(self) -> Self::Neighborhood {
let [a, b] = self;
[
[a - 1, b - 1],
[a, b - 1],
[a + 1, b - 1],
[a - 1, b],
[a, b],
[a + 1, b],
[a - 1, b + 1],
[a, b + 1],
[a + 1, b + 1],
]
}
}
impl IntVector for [i64; 3] {
type Neighborhood = [[i64; 3]; 27];
fn neighborhood(self) -> Self::Neighborhood {
let [a, b, c] = self;
[
[a - 1, b - 1, c - 1],
[a, b - 1, c - 1],
[a + 1, b - 1, c - 1],
[a - 1, b, c - 1],
[a, b, c - 1],
[a + 1, b, c - 1],
[a - 1, b + 1, c - 1],
[a, b + 1, c - 1],
[a + 1, b + 1, c - 1],
[a - 1, b - 1, c],
[a, b - 1, c],
[a + 1, b - 1, c],
[a - 1, b, c],
[a, b, c],
[a + 1, b, c],
[a - 1, b + 1, c],
[a, b + 1, c],
[a + 1, b + 1, c],
[a - 1, b - 1, c + 1],
[a, b - 1, c + 1],
[a + 1, b - 1, c + 1],
[a - 1, b, c + 1],
[a, b, c + 1],
[a + 1, b, c + 1],
[a - 1, b + 1, c + 1],
[a, b + 1, c + 1],
[a + 1, b + 1, c + 1],
]
}
}
#[cfg(test)]
mod tests {
use super::*;
fn into_vertices(iter: &[usize]) -> Vec<Vertex> { iter.iter().map(|i| i.into()).collect() }
#[test]
fn degenerate_polygon_test() {
let poly = into_vertices(&[0, 1, 2, 0, 3, 4, 5, 6, 3, 7, 8, 9]);
let polys = split_into_nondegenerate(poly);
assert_eq!(polys[0], into_vertices(&[0, 1, 2]));
assert_eq!(polys[1], into_vertices(&[3, 4, 5, 6]));
assert_eq!(polys[2], into_vertices(&[3, 7, 8, 9, 0]));
}
}
+242
View File
@@ -0,0 +1,242 @@
use super::*;
/// triangulation, quadrangulation, give a structure
pub trait StructuringFilter {
/// triangulate all n-gons
/// # Examples
/// ```
/// use truck_polymesh::*;
/// use truck_meshalgo::filters::*;
///
/// // cube consisting quad faces
/// let mut mesh = PolygonMesh::new(
/// StandardAttributes {
/// positions: vec![
/// Point3::new(0.0, 0.0, 0.0),
/// Point3::new(1.0, 0.0, 0.0),
/// Point3::new(1.0, 1.0, 0.0),
/// Point3::new(0.0, 1.0, 0.0),
/// Point3::new(0.0, 0.0, 1.0),
/// Point3::new(1.0, 0.0, 1.0),
/// Point3::new(1.0, 1.0, 1.0),
/// Point3::new(0.0, 1.0, 1.0),
/// ],
/// ..Default::default()
/// },
/// Faces::from_iter(&[
/// &[3, 2, 1, 0], &[0, 1, 5, 4], &[1, 2, 6, 5],
/// &[2, 3, 7, 6], &[3, 0, 4, 7], &[4, 5, 6, 7],
/// ]),
/// );
///
/// // the number of face becomes twice since each quadrangle decompose into two triangles.
/// assert_eq!(mesh.faces().len(), 6);
/// mesh.triangulate();
/// assert_eq!(mesh.faces().len(), 12);
/// ```
fn triangulate(&mut self) -> &mut Self;
/// join two triangles into one quadrangle.
/// # Arguments
/// * `plane_tol` - the tolerance for determining that four points are in the same plane
/// * `score_tol` - The upper limit of the score to determine if the four points form an uncrushed rectangle
/// # Details
/// The overview of the algorithm is the following:
/// 1. make the list of pairs of triangles satisfying the following conditions:
/// * two faces are adjacent by one edge,
/// * the pair of faces consists four vertices in the same plane
/// 1. for each joined quadrangle, calculate the score by the sum of the absolute value of
/// the cosine for each angles,
/// 1. sort the list of the pairs of triangles by the score
/// 1. take a pair of triangles in order from the top of the list and register a new one
/// if it doesn't conflict with the one has been already registered.
/// # Panics
/// `plane_tol` and `score_tol` must be more than `TOLERANCE`.
/// # Examples
/// ```
/// use truck_polymesh::*;
/// use truck_meshalgo::filters::*;
///
/// // cube consisting tri_faces
/// let mut mesh = PolygonMesh::new(
/// StandardAttributes {
/// positions: vec![
/// Point3::new(0.0, 0.0, 0.0),
/// Point3::new(1.0, 0.0, 0.0),
/// Point3::new(1.0, 1.0, 0.0),
/// Point3::new(0.0, 1.0, 0.0),
/// Point3::new(0.0, 0.0, 1.0),
/// Point3::new(1.0, 0.0, 1.0),
/// Point3::new(1.0, 1.0, 1.0),
/// Point3::new(0.0, 1.0, 1.0),
/// ],
/// ..Default::default()
/// },
/// Faces::from_iter(&[
/// &[3, 2, 0], &[1, 0, 2], &[0, 1, 4], &[5, 4, 1],
/// &[1, 2, 5], &[6, 5, 2], &[2, 3, 6], &[7, 6, 3],
/// &[3, 0, 7], &[4, 7, 0], &[4, 5, 7], &[6, 7, 5],
/// ]),
/// );
///
/// // The number of faces becomes a half since each pair of triangles is combined.
/// assert_eq!(mesh.faces().len(), 12);
/// mesh.quadrangulate(0.01, 0.1);
/// assert_eq!(mesh.faces().len(), 6);
/// ```
fn quadrangulate(&mut self, plane_tol: f64, score_tol: f64) -> &mut Self;
}
impl StructuringFilter for PolygonMesh {
fn triangulate(&mut self) -> &mut Self {
let tri_faces = self.faces().triangle_iter().collect::<Vec<_>>();
*self.debug_editor().faces = Faces::from_tri_and_quad_faces(tri_faces, Vec::new());
self
}
fn quadrangulate(&mut self, plane_tol: f64, score_tol: f64) -> &mut Self {
nonpositive_tolerance!(plane_tol, 0.0);
nonpositive_tolerance!(score_tol, 0.0);
let list = self.create_face_edge_list(plane_tol, score_tol);
self.reflect_face_edge_list(list);
self
}
}
trait SubStructureFilter {
fn create_face_edge_list(&self, plane_tol: f64, score_tol: f64) -> Vec<FaceEdge>;
fn reflect_face_edge_list(&mut self, list: Vec<FaceEdge>);
fn get_face_edge(
&self,
face0_id: usize,
face1_id: usize,
plane_tol: f64,
score_tol: f64,
) -> Option<FaceEdge>;
}
impl SubStructureFilter for PolygonMesh {
fn create_face_edge_list(&self, plane_tol: f64, score_tol: f64) -> Vec<FaceEdge> {
let face_adjacency = self.faces().face_adjacency(true);
let mut passed = Vec::new();
let tri_len = self.faces().tri_faces().len();
for (i, face) in face_adjacency.iter().enumerate().take(tri_len) {
for j in face {
if i > *j {
continue;
} else if let Some(face_edge) = self.get_face_edge(i, *j, plane_tol, score_tol) {
passed.push(face_edge);
}
}
}
passed.sort_by(|x, y| x.score.partial_cmp(&y.score).unwrap());
passed
}
fn reflect_face_edge_list(&mut self, list: Vec<FaceEdge>) {
let mut used = vec![false; self.faces().tri_faces().len()];
let mut quad_faces = self.faces().quad_faces().clone();
quad_faces.extend(list.into_iter().filter_map(|face_edge| {
let (i, j) = face_edge.faces;
if used[i] || used[j] {
None
} else {
used[i] = true;
used[j] = true;
Some(face_edge.positions)
}
}));
let tri_faces = self.faces().tri_faces();
let tri_faces = used
.into_iter()
.enumerate()
.filter_map(move |(i, flag)| match flag {
true => None,
false => Some(tri_faces[i]),
})
.collect::<Vec<_>>();
*self.debug_editor().faces = Faces::from_tri_and_quad_faces(tri_faces, quad_faces);
}
fn get_face_edge(
&self,
face0_id: usize,
face1_id: usize,
plane_tol: f64,
score_tol: f64,
) -> Option<FaceEdge> {
let face0 = self.faces().tri_faces()[face0_id];
let face1 = self.faces().tri_faces()[face1_id];
let k = (0..3)
.find(|k| face0.iter().all(|x| x.pos != face1[*k].pos))
.unwrap();
let vec0 = self.positions()[face0[1].pos] - self.positions()[face0[0].pos];
let vec1 = self.positions()[face0[2].pos] - self.positions()[face0[0].pos];
let mut n = vec0.cross(vec1);
n /= n.magnitude();
let vec2 = self.positions()[face1[k].pos] - self.positions()[face0[0].pos];
let mat = Matrix3::from_cols(vec0, vec1, n);
let coef = mat.invert().unwrap() * vec2;
if coef[2] > plane_tol {
None
} else if coef[0] > 0.0 && coef[1] > 0.0 {
let score = calc_score(vec0, vec2 - vec0, vec1 - vec2, vec1);
if score < score_tol {
Some(FaceEdge {
faces: (face0_id, face1_id),
positions: [face0[0], face0[1], face1[k], face0[2]],
score,
})
} else {
None
}
} else if coef[0] < 0.0 && coef[1] > 0.0 && coef[0] + coef[1] < 1.0 {
let score = calc_score(vec0, vec1 - vec0, vec2 - vec1, vec2);
if score < score_tol {
Some(FaceEdge {
faces: (face0_id, face1_id),
positions: [face0[0], face0[1], face0[2], face1[k]],
score,
})
} else {
None
}
} else if coef[0] > 0.0 && coef[1] < 0.0 && coef[0] + coef[1] < 1.0 {
let score = calc_score(vec2, vec0 - vec2, vec1 - vec0, vec1);
if score < score_tol {
Some(FaceEdge {
faces: (face0_id, face1_id),
positions: [face0[0], face1[k], face0[1], face0[2]],
score,
})
} else {
None
}
} else {
None
}
}
}
struct FaceEdge {
faces: (usize, usize),
positions: [Vertex; 4],
score: f64,
}
#[inline(always)]
fn calc_score(edge0: Vector3, edge1: Vector3, edge2: Vector3, edge3: Vector3) -> f64 {
edge0.cos_angle(edge1).abs()
+ edge1.cos_angle(edge2).abs()
+ edge2.cos_angle(edge3).abs()
+ edge3.cos_angle(edge0).abs()
}
trait CosAngle {
fn cos_angle(self, other: Self) -> f64;
}
impl CosAngle for Vector3 {
fn cos_angle(self, other: Self) -> f64 {
self.dot(other) / (self.magnitude() * other.magnitude())
}
}
+199
View File
@@ -0,0 +1,199 @@
use super::*;
use rustc_hash::FxHashMap as HashMap;
use std::f64::consts::PI;
/// subdivision surface algorithms
pub trait Subdivision {
/// Extended Loop method
///
/// # Remarks
/// Confirm:
/// - All faces are triangles.
/// - `self.shell_condition()` is `Oriented` or `Closed` before use.
///
/// This method does NOT check these conditions.
fn loop_subdivision(&mut self) -> &mut Self;
}
#[derive(Clone, Copy, Debug, PartialEq, Eq, PartialOrd, Ord, Hash)]
struct Edge(usize, usize);
#[derive(Clone, Copy, Debug)]
struct EdgeInfo {
idx: usize,
first_wing: usize,
second_wing: Option<usize>,
}
impl Edge {
#[inline]
fn new(v0: usize, v1: usize) -> Self { Edge(usize::min(v0, v1), usize::max(v0, v1)) }
}
impl EdgeInfo {
#[inline]
fn new(idx: usize, first_wing: usize) -> Self {
Self {
idx,
first_wing,
second_wing: None,
}
}
#[inline]
fn add_second_wing(&mut self, second_wing: usize) { self.second_wing = Some(second_wing); }
}
impl Subdivision for PolygonMesh {
fn loop_subdivision(&mut self) -> &mut Self {
let mut edges = HashMap::default();
let mut vertex_adjacency = vec![Vec::new(); self.positions().len()];
self.tri_faces().iter().for_each(|face| {
let face = [face[0].pos, face[1].pos, face[2].pos];
add_vertex_edge(&mut edges, &mut vertex_adjacency, face[1], face[2], face[0]);
add_vertex_edge(&mut edges, &mut vertex_adjacency, face[2], face[0], face[1]);
add_vertex_edge(&mut edges, &mut vertex_adjacency, face[0], face[1], face[2]);
});
let positions = vertex_adjacency
.iter()
.enumerate()
.map(|t| calc_new_position(t.0, t.1, &edges, self.positions()))
.chain(edge_positions(&edges, &vertex_adjacency, self.positions()))
.collect::<Vec<_>>();
let tri_faces = self
.tri_faces()
.iter()
.flat_map(|v| {
let len = self.positions().len();
let e: [StandardVertex; 3] = [
(edges.get(&Edge::new(v[1].pos, v[2].pos)).unwrap().idx + len).into(),
(edges.get(&Edge::new(v[2].pos, v[0].pos)).unwrap().idx + len).into(),
(edges.get(&Edge::new(v[0].pos, v[1].pos)).unwrap().idx + len).into(),
];
vec![
[v[0], e[2], e[1]],
[e[0], v[2], e[1]],
[e[0], e[2], v[1]],
[e[0], e[1], e[2]],
]
})
.collect::<Vec<_>>();
{
let editor = self.debug_editor();
editor.attributes.positions = positions;
*editor.faces = Faces::from_tri_and_quad_faces(tri_faces, Vec::new());
}
self
}
}
fn add_vertex_edge(
edges: &mut HashMap<Edge, EdgeInfo>,
vertex_adjacency: &mut [Vec<usize>],
v0: usize,
v1: usize,
v2: usize,
) {
let idx = edges.len();
edges
.entry(Edge::new(v0, v1))
.and_modify(|tuple| tuple.add_second_wing(v2))
.or_insert_with(|| {
vertex_adjacency[v0].push(v1);
vertex_adjacency[v1].push(v0);
EdgeInfo::new(idx, v2)
});
}
#[derive(Clone, Copy, Debug)]
enum VertexBoundaryCondition {
#[allow(dead_code)]
Corner(usize, usize),
Boundary(usize, usize),
Inner,
}
impl VertexBoundaryCondition {
fn new(v: usize, adjacency: &[usize], edges: &HashMap<Edge, EdgeInfo>) -> Self {
let binfo = adjacency.iter().copied().fold((None, None), |binfo, w| {
let edge = edges.get(&Edge::new(v, w)).unwrap();
match (edge.second_wing, binfo) {
(None, (None, _)) => (Some(w), None),
(None, (Some(x), None)) => (Some(x), Some(w)),
_ => binfo,
}
});
if let (Some(x), Some(y)) = binfo {
if adjacency.len() == 2 {
Self::Corner(x, y)
} else {
Self::Boundary(x, y)
}
} else {
Self::Inner
}
}
}
fn calc_new_position(
v: usize,
adjacency: &[usize],
edges: &HashMap<Edge, EdgeInfo>,
positions: &[Point3],
) -> Point3 {
use VertexBoundaryCondition::*;
let point = positions[v];
match VertexBoundaryCondition::new(v, adjacency, edges) {
Corner(_, _) => point,
Boundary(w0, w1) => point + ((positions[w0] - point) + (positions[w1] - point)) / 8.0,
Inner => {
let alpha = 3.0 / 8.0 + f64::cos(2.0 * PI / adjacency.len() as f64) / 4.0;
let alpha = (5.0 / 8.0 - alpha * alpha) / adjacency.len() as f64;
point
+ adjacency
.iter()
.copied()
.map(|w| positions[w] - point)
.sum::<Vector3>()
* alpha
}
}
}
fn edge_positions(
edges: &HashMap<Edge, EdgeInfo>,
vertex_adjacency: &[Vec<usize>],
positions: &[Point3],
) -> Vec<Point3> {
use VertexBoundaryCondition::*;
let mut res = vec![Point3::origin(); edges.len()];
edges.iter().for_each(|tuple| {
let Edge(v0, v1) = *tuple.0;
let EdgeInfo {
idx,
first_wing,
second_wing,
} = *tuple.1;
res[idx] = match second_wing {
Some(second_wing) => {
let gamma = match VertexBoundaryCondition::new(v0, &vertex_adjacency[v0], edges) {
Boundary(_, _) => {
5.0 / 8.0 - f64::cos(PI / vertex_adjacency[v0].len() as f64) / 4.0
}
Inner => match VertexBoundaryCondition::new(v1, &vertex_adjacency[v1], edges) {
Boundary(_, _) => {
3.0 / 8.0 + f64::cos(PI / vertex_adjacency[v1].len() as f64) / 4.0
}
Inner => 0.5,
Corner(_, _) => unreachable!(),
},
Corner(_, _) => unreachable!(),
};
positions[v0]
+ (positions[v1] - positions[v0]) * gamma
+ (positions[first_wing] - positions[v0]) / 8.0
+ (positions[second_wing] - positions[v1]) / 8.0
}
None => positions[v0].midpoint(positions[v1]),
};
});
res
}
+56
View File
@@ -0,0 +1,56 @@
//! Mesh algorigthms, include tessellations of the shape.
// SolveCraft: compiler warnings no longer fail the build.
#![deny(clippy::all, rust_2018_idioms)]
#![warn(
missing_docs,
missing_debug_implementations,
trivial_casts,
trivial_numeric_casts,
unsafe_code,
unstable_features,
unused_import_braces,
unused_qualifications
)]
use common::*;
/// re-export polymesh
pub mod rexport_polymesh {
pub use truck_polymesh::*;
}
use truck_polymesh::{StandardVertex as Vertex, *};
/// polygon mesh analizers, including
///
/// - determines topological properties: connectivity, boundary extraction, or shell conditions (closed or oriented)
/// - detects collisions between two meshes and extracts interference lines
/// - investigates positional relations between mesh and point clouds.
#[cfg(feature = "analyzers")]
pub mod analyzers;
mod common;
/// Edits meshes. Add normals, optimizing data, and so on.
#[cfg(feature = "filters")]
pub mod filters;
/// Tessellates shapes.
#[cfg(feature = "tessellation")]
pub mod tessellation;
/// VTK Output
#[cfg(feature = "vtk")]
#[cfg(not(target_arch = "wasm32"))]
pub mod vtk;
/// This module contains all traits and re-exports `truck_polymesh`.
pub mod prelude {
#[cfg(feature = "analyzers")]
pub use crate::analyzers::*;
#[cfg(feature = "filters")]
pub use crate::filters::*;
pub use crate::rexport_polymesh::*;
#[cfg(feature = "tessellation")]
pub use crate::tessellation::*;
#[cfg(feature = "vtk")]
#[cfg(not(target_arch = "wasm32"))]
pub use crate::vtk::*;
}
+360
View File
@@ -0,0 +1,360 @@
use crate::*;
use spade::{iterators::*, *};
use truck_topology::{compress::*, *};
#[cfg(not(target_arch = "wasm32"))]
mod parallelizable {
/// Parallelizable by `rayon`.
pub trait Parallelizable: Send + Sync {}
impl<T: Send + Sync> Parallelizable for T {}
}
#[cfg(target_arch = "wasm32")]
mod parallelizable {
/// No parallelization in the case of wasm.
pub trait Parallelizable {}
impl<T> Parallelizable for T {}
}
pub use parallelizable::*;
/// Gathered the traits used in tessellation.
pub trait PolylineableCurve:
ParametricCurve3D + BoundedCurve + ParameterDivision1D<Point = Point3> + Parallelizable {
}
impl<
C: ParametricCurve3D + BoundedCurve + ParameterDivision1D<Point = Point3> + Parallelizable,
> PolylineableCurve for C
{
}
/// It can be meshed, but not necessarily trimmed.
pub trait PreMeshableSurface: ParametricSurface3D + ParameterDivision2D + Parallelizable {}
impl<S: ParametricSurface3D + ParameterDivision2D + Parallelizable> PreMeshableSurface for S {}
/// The generated mesh can be trimmed only if the boundary curves ride strictly on a surface.
pub trait MeshableSurface: PreMeshableSurface + SearchParameter<D2, Point = Point3> {}
impl<S: PreMeshableSurface + SearchParameter<D2, Point = Point3>> MeshableSurface for S {}
/// The generated mesh can be trimmed if the boundary curves does not ride strictly on a surface.
pub trait RobustMeshableSurface:
MeshableSurface + SearchNearestParameter<D2, Point = Point3> {
}
impl<S: MeshableSurface + SearchNearestParameter<D2, Point = Point3>> RobustMeshableSurface for S {}
type PolylineCurve = truck_polymesh::PolylineCurve<Point3>;
/// Trait for converting tessellated shape into polygon.
pub trait MeshedShape {
/// Converts tessellated shape into polygon.
fn to_polygon(&self) -> PolygonMesh;
}
impl MeshedShape for Shell<Point3, PolylineCurve, PolygonMesh> {
fn to_polygon(&self) -> PolygonMesh {
let mut polygon = PolygonMesh::default();
self.face_iter().for_each(|face| {
polygon.merge(face.oriented_surface());
});
polygon
}
}
impl MeshedShape for Shell<Point3, PolylineCurve, Option<PolygonMesh>> {
fn to_polygon(&self) -> PolygonMesh {
let mut polygon = PolygonMesh::default();
self.face_iter().for_each(|face| {
if let Some(mut poly) = face.surface() {
if !face.orientation() {
poly.invert();
}
polygon.merge(poly);
}
});
polygon
}
}
impl<P, C, S> MeshedShape for Solid<P, C, S>
where Shell<P, C, S>: MeshedShape
{
fn to_polygon(&self) -> PolygonMesh {
let mut polygon = PolygonMesh::default();
self.boundaries().iter().for_each(|shell| {
polygon.merge(shell.to_polygon());
});
polygon
}
}
impl MeshedShape for CompressedShell<Point3, PolylineCurve, PolygonMesh> {
fn to_polygon(&self) -> PolygonMesh {
let mut polygon = PolygonMesh::default();
self.faces.iter().for_each(|face| match face.orientation {
true => polygon.merge(face.surface.clone()),
false => polygon.merge(face.surface.inverse()),
});
polygon
}
}
impl MeshedShape for CompressedShell<Point3, PolylineCurve, Option<PolygonMesh>> {
fn to_polygon(&self) -> PolygonMesh {
let mut polygon = PolygonMesh::default();
self.faces.iter().for_each(|face| {
if let Some(surface) = &face.surface {
match face.orientation {
true => polygon.merge(surface.clone()),
false => polygon.merge(surface.inverse()),
}
}
});
polygon
}
}
impl<P, C, S> MeshedShape for CompressedSolid<P, C, S>
where CompressedShell<P, C, S>: MeshedShape
{
fn to_polygon(&self) -> PolygonMesh {
let mut polygon = PolygonMesh::default();
self.boundaries.iter().for_each(|shell| {
polygon.merge(shell.to_polygon());
});
polygon
}
}
/// Trait for tessellating `Shell` and `Solid`.
pub trait MeshableShape {
/// Shape whose edges are made polylines and faces polygon surface.
type MeshedShape: MeshedShape;
/// Tessellates shapes. The division of curves and surfaces are by `ParameterDivision1D` and `ParameterDivision2D`,
/// and the constrained Delauney triangulation is based on the crate [`spade`](https://crates.io/crates/spade).
///
/// # Panics
///
/// `tol` must be more than `TOLERANCE`.
///
/// # Remarks
///
/// - The tessellated mesh is not necessarily closed even if `self` is `Solid`.
/// If you want to get closed mesh, use [`OptimizingFilter::put_together_same_attrs`].
/// - This method requires that the curve ride strictly on a surface. If not, try [`RobustMeshableShape`].
///
/// [`OptimizingFilter::put_together_same_attrs`]: crate::filters::OptimizingFilter::put_together_same_attrs
///
/// # Examples
/// ```
/// use truck_meshalgo::prelude::*;
/// use truck_modeling::builder;
/// use truck_topology::shell::ShellCondition;
///
/// // modeling a unit cube
/// let v = builder::vertex(Point3::origin());
/// let e = builder::tsweep(&v, Vector3::unit_x());
/// let f = builder::tsweep(&e, Vector3::unit_y());
/// let cube = builder::tsweep(&f, Vector3::unit_z());
///
/// // cube is Solid, however, the tessellated mesh is not closed.
/// let mut mesh = cube.triangulation(0.01).to_polygon();
/// assert_ne!(mesh.shell_condition(), ShellCondition::Closed);
///
/// // use optimization filters!
/// mesh.put_together_same_attrs(TOLERANCE);
/// assert_eq!(mesh.shell_condition(), ShellCondition::Closed);
/// ```
fn triangulation(&self, tol: f64) -> Self::MeshedShape;
}
/// Trait for tessellating `Shell` and `Solid` in `truck-modeling`.
pub trait RobustMeshableShape {
/// Shape whose edges are made polylines and faces polygon surface.
type MeshedShape: MeshedShape;
/// Tessellates shapes. The division of curves and surfaces are by `ParameterDivision1D` and `ParameterDivision2D`,
/// and the constrained Delauney triangulation is based on the crate [`spade`](https://crates.io/crates/spade).
///
/// # Panics
///
/// `tol` must be more than `TOLERANCE`.
///
/// # Remarks
///
/// Since polyline vertices are projected onto the surface, processing speed is often slower than with [`MeshableShape::triangulation`].
/// It also does not close the mesh of a solid even if one uses [`OptimizingFilter::put_together_same_attrs`].
///
/// [`OptimizingFilter::put_together_same_attrs`]: crate::filters::OptimizingFilter::put_together_same_attrs
///
/// # Examples
/// ```
/// use truck_meshalgo::prelude::*;
/// use truck_modeling::*;
/// use truck_topology::shell::ShellCondition;
///
/// // manual modeling an open half cylinder
///
/// // points
/// let p = [
/// Point3::new(1.0, 0.0, 0.0),
/// Point3::new(-1.0, 0.0, 0.0),
/// Point3::new(1.0, 0.0, 1.0),
/// Point3::new(-1.0, 0.0, 1.0),
/// ];
/// // vertices
/// let v = Vertex::news(&p);
/// // Curves that do not ride on a cylinder
/// let bsp0: Curve = BSplineCurve::new(
/// KnotVec::bezier_knot(3),
/// vec![
/// p[0],
/// Point3::new(1.0, 4.0 / 3.0, 0.0),
/// Point3::new(-1.0, 4.0 / 3.0, 0.0),
/// p[1],
/// ],
/// )
/// .into();
/// let bsp1: Curve = BSplineCurve::new(
/// KnotVec::bezier_knot(3),
/// vec![
/// p[3],
/// Point3::new(-1.0, 4.0 / 3.0, 1.0),
/// Point3::new(1.0, 4.0 / 3.0, 1.0),
/// p[2],
/// ],
/// )
/// .into();
/// // wire
/// let w: Wire = vec![
/// builder::line(&v[2], &v[0]),
/// Edge::new(&v[0], &v[1], bsp0),
/// builder::line(&v[1], &v[3]),
/// Edge::new(&v[3], &v[2], bsp1),
/// ]
/// .into();
/// // revoluted curve
/// let surface_raw = RevolutedCurve::by_revolution(
/// Curve::Line(Line(p[2], p[0])),
/// Point3::origin(),
/// Vector3::unit_z(),
/// );
/// let surface: Surface = Processor::new(surface_raw).into();
/// // shell
/// let shell: Shell = vec![Face::new(vec![w], surface)].into();
///
/// // Simple triangulation fails since some edges do not ride on a cylinder
/// let poly_shell = shell.triangulation(0.01);
/// assert!(poly_shell[0].surface().is_none());
///
/// // Robust triangulation!
/// let poly_shell = shell.robust_triangulation(0.01);
/// let poly = poly_shell[0].surface().unwrap();
/// assert!(!poly.positions().is_empty());
/// ```
fn robust_triangulation(&self, tol: f64) -> Self::MeshedShape;
}
impl<C: PolylineableCurve, S: MeshableSurface> MeshableShape for Shell<Point3, C, S> {
type MeshedShape = Shell<Point3, PolylineCurve, Option<PolygonMesh>>;
fn triangulation(&self, tol: f64) -> Self::MeshedShape {
nonpositive_tolerance!(tol);
#[cfg(not(target_arch = "wasm32"))]
let res = triangulation::shell_tessellation(self, tol, triangulation::by_search_parameter);
#[cfg(target_arch = "wasm32")]
let res = triangulation::shell_tessellation_single_thread(
self,
tol,
triangulation::by_search_parameter,
);
res
}
}
impl<C: PolylineableCurve, S: RobustMeshableSurface> RobustMeshableShape for Shell<Point3, C, S> {
type MeshedShape = Shell<Point3, PolylineCurve, Option<PolygonMesh>>;
fn robust_triangulation(&self, tol: f64) -> Self::MeshedShape {
nonpositive_tolerance!(tol);
#[cfg(not(target_arch = "wasm32"))]
let res = triangulation::shell_tessellation(
self,
tol,
triangulation::by_search_nearest_parameter,
);
#[cfg(target_arch = "wasm32")]
let res = triangulation::shell_tessellation_single_thread(
self,
tol,
triangulation::by_search_nearest_parameter,
);
res
}
}
impl<C: PolylineableCurve, S: MeshableSurface> MeshableShape for Solid<Point3, C, S> {
type MeshedShape = Solid<Point3, PolylineCurve, Option<PolygonMesh>>;
fn triangulation(&self, tol: f64) -> Self::MeshedShape {
let boundaries = self
.boundaries()
.iter()
.map(|shell| shell.triangulation(tol))
.collect::<Vec<_>>();
Solid::new(boundaries)
}
}
impl<C: PolylineableCurve, S: RobustMeshableSurface> RobustMeshableShape for Solid<Point3, C, S> {
type MeshedShape = Solid<Point3, PolylineCurve, Option<PolygonMesh>>;
fn robust_triangulation(&self, tol: f64) -> Self::MeshedShape {
let boundaries = self
.boundaries()
.iter()
.map(|shell| shell.robust_triangulation(tol))
.collect::<Vec<_>>();
Solid::new(boundaries)
}
}
impl<C: PolylineableCurve, S: MeshableSurface> MeshableShape for CompressedShell<Point3, C, S> {
type MeshedShape = CompressedShell<Point3, PolylineCurve, Option<PolygonMesh>>;
fn triangulation(&self, tol: f64) -> Self::MeshedShape {
nonpositive_tolerance!(tol);
triangulation::cshell_tessellation(self, tol, triangulation::by_search_parameter)
}
}
impl<C: PolylineableCurve, S: RobustMeshableSurface> RobustMeshableShape
for CompressedShell<Point3, C, S>
{
type MeshedShape = CompressedShell<Point3, PolylineCurve, Option<PolygonMesh>>;
fn robust_triangulation(&self, tol: f64) -> Self::MeshedShape {
nonpositive_tolerance!(tol);
triangulation::cshell_tessellation(self, tol, triangulation::by_search_nearest_parameter)
}
}
impl<C: PolylineableCurve, S: MeshableSurface> MeshableShape for CompressedSolid<Point3, C, S> {
type MeshedShape = CompressedSolid<Point3, PolylineCurve, Option<PolygonMesh>>;
fn triangulation(&self, tol: f64) -> Self::MeshedShape {
let boundaries = self
.boundaries
.iter()
.map(|shell| shell.triangulation(tol))
.collect::<Vec<_>>();
CompressedSolid { boundaries }
}
}
impl<C: PolylineableCurve, S: RobustMeshableSurface> RobustMeshableShape
for CompressedSolid<Point3, C, S>
{
type MeshedShape = CompressedSolid<Point3, PolylineCurve, Option<PolygonMesh>>;
fn robust_triangulation(&self, tol: f64) -> Self::MeshedShape {
let boundaries = self
.boundaries
.iter()
.map(|shell| shell.robust_triangulation(tol))
.collect::<Vec<_>>();
CompressedSolid { boundaries }
}
}
mod triangulation;
+705
View File
@@ -0,0 +1,705 @@
#![allow(clippy::many_single_char_names)]
use super::*;
use crate::filters::NormalFilters;
use crate::Point2;
use array_macro::array;
use itertools::Itertools;
use rustc_hash::FxHashMap as HashMap;
#[cfg(not(target_arch = "wasm32"))]
use rayon::prelude::*;
type SPoint2 = spade::Point2<f64>;
type Cdt = ConstrainedDelaunayTriangulation<SPoint2>;
type MeshedShell = Shell<Point3, PolylineCurve, Option<PolygonMesh>>;
type MeshedCShell = CompressedShell<Point3, PolylineCurve, Option<PolygonMesh>>;
pub(super) trait SP<S>:
Fn(&S, Point3, Option<(f64, f64)>) -> Option<(f64, f64)> + Parallelizable {
}
impl<S, F> SP<S> for F where F: Fn(&S, Point3, Option<(f64, f64)>) -> Option<(f64, f64)> + Parallelizable {}
/// SolveCraft: a parameter inside the surface's bounded, non-periodic ranges.
fn in_domain<S: ParametricSurface>(surface: &S, (u, v): (f64, f64)) -> bool {
let ok = |x: f64, range: Option<(f64, f64)>, period: Option<f64>| {
x.is_finite()
&& match (range, period) {
(_, Some(_)) | (None, None) => true,
(Some((a, b)), None) => {
let m = (b - a).abs() * 1e-6 + 1e-9;
x >= a.min(b) - m && x <= a.max(b) + m
}
}
};
let (ur, vr) = surface.try_range_tuple();
ok(u, ur, surface.u_period()) && ok(v, vr, surface.v_period())
}
/// SolveCraft: the parameter clamped into the surface's bounded, non-periodic ranges.
fn clamped<S: ParametricSurface>(surface: &S, (u, v): (f64, f64)) -> (f64, f64) {
let c = |x: f64, range: Option<(f64, f64)>, period: Option<f64>| match (range, period) {
(Some((a, b)), None) => x.clamp(a.min(b), a.max(b)),
_ => x,
};
let (ur, vr) = surface.try_range_tuple();
(c(u, ur, surface.u_period()), c(v, vr, surface.v_period()))
}
/// SolveCraft: a nearest-parameter search can run off a B-spline's domain and extrapolate to
/// thousands of parameter units (meshing that grid never finishes). Outside the nominal domain
/// keep the result only if it is really closer to the point than the clamped one (planes report
/// a nominal unit range but extend for ever).
fn sane<S: RobustMeshableSurface>(surface: &S, point: Point3, uv: (f64, f64)) -> (f64, f64) {
if in_domain(surface, uv) {
return uv;
}
let c = clamped(surface, uv);
let d = |(u, v): (f64, f64)| surface.subs(u, v).distance2(point);
if d(uv) <= d(c) { uv } else { c }
}
pub(super) fn by_search_parameter<S>(
surface: &S,
point: Point3,
hint: Option<(f64, f64)>,
) -> Option<(f64, f64)>
where
S: MeshableSurface,
{
surface
.search_parameter(point, hint, 100)
.or_else(|| surface.search_parameter(point, None, 100))
}
pub(super) fn by_search_nearest_parameter<S>(
surface: &S,
point: Point3,
hint: Option<(f64, f64)>,
) -> Option<(f64, f64)>
where
S: RobustMeshableSurface,
{
surface
.search_parameter(point, hint, 100)
.or_else(|| surface.search_parameter(point, None, 100))
.or_else(|| surface.search_nearest_parameter(point, hint, 100).map(|uv| sane(surface, point, uv)))
.or_else(|| surface.search_nearest_parameter(point, None, 100).map(|uv| sane(surface, point, uv)))
}
/// Tessellates faces
#[cfg(not(target_arch = "wasm32"))]
pub(super) fn shell_tessellation<'a, C, S>(
shell: &Shell<Point3, C, S>,
tol: f64,
sp: impl SP<S>,
) -> MeshedShell
where
C: PolylineableCurve + 'a,
S: PreMeshableSurface + 'a,
{
let vmap: HashMap<_, _> = shell
.vertex_par_iter()
.map(|v| (v.id(), v.mapped(Point3::clone)))
.collect();
let eset: HashMap<_, _> = shell.edge_par_iter().map(move |e| (e.id(), e)).collect();
let edge_map: HashMap<_, _> = eset
.into_par_iter()
.map(move |(id, edge)| {
let v0 = vmap.get(&edge.absolute_front().id()).unwrap();
let v1 = vmap.get(&edge.absolute_back().id()).unwrap();
let curve = edge.curve();
let poly = PolylineCurve::from_curve(&curve, curve.range_tuple(), tol);
(id, Edge::debug_new(v0, v1, poly))
})
.collect();
let create_edge = |edge: &Edge<Point3, C>| -> Edge<_, _> {
let new_edge = edge_map.get(&edge.id()).unwrap();
match edge.orientation() {
true => new_edge.clone(),
false => new_edge.inverse(),
}
};
let create_boundary =
|wire: &Wire<Point3, C>| -> Wire<_, _> { wire.edge_iter().map(create_edge).collect() };
let create_face = move |face: &Face<Point3, C, S>| -> Face<_, _, _> {
let wires: Vec<_> = face
.absolute_boundaries()
.iter()
.map(create_boundary)
.collect();
shell_create_polygon(&face.surface(), wires, face.orientation(), tol, &sp)
};
shell.face_par_iter().map(create_face).collect()
}
/// Tessellates faces
#[cfg(any(target_arch = "wasm32", test))]
pub(super) fn shell_tessellation_single_thread<'a, C, S>(
shell: &'a Shell<Point3, C, S>,
tol: f64,
sp: impl SP<S>,
) -> MeshedShell
where
C: PolylineableCurve + 'a,
S: PreMeshableSurface + 'a,
{
use truck_base::entry_map::FxEntryMap as EntryMap;
use truck_topology::Vertex as TVertex;
let mut vmap = EntryMap::new(
move |v: &TVertex<Point3>| v.id(),
move |v| v.mapped(Point3::clone),
);
let mut edge_map = EntryMap::new(
move |edge: &'a Edge<Point3, C>| edge.id(),
move |edge| {
let vf = edge.absolute_front();
let v0 = vmap.entry_or_insert(vf).clone();
let vb = edge.absolute_back();
let v1 = vmap.entry_or_insert(vb).clone();
let curve = edge.curve();
let poly = PolylineCurve::from_curve(&curve, curve.range_tuple(), tol);
Edge::debug_new(&v0, &v1, poly)
},
);
let mut create_edge = move |edge: &'a Edge<Point3, C>| -> Edge<_, _> {
let new_edge = edge_map.entry_or_insert(edge);
match edge.orientation() {
true => new_edge.clone(),
false => new_edge.inverse(),
}
};
let mut create_boundary = move |wire: &'a Wire<Point3, C>| -> Wire<_, _> {
wire.edge_iter().map(&mut create_edge).collect()
};
let create_face = move |face: &'a Face<Point3, C, S>| -> Face<_, _, _> {
let wires: Vec<_> = face
.absolute_boundaries()
.iter()
.map(&mut create_boundary)
.collect();
shell_create_polygon(&face.surface(), wires, face.orientation(), tol, &sp)
};
shell.face_iter().map(create_face).collect()
}
/// Tessellates faces
pub(super) fn cshell_tessellation<'a, C, S>(
shell: &CompressedShell<Point3, C, S>,
tol: f64,
sp: impl SP<S>,
) -> MeshedCShell
where
C: PolylineableCurve + 'a,
S: PreMeshableSurface + 'a,
{
let vertices = shell.vertices.clone();
let tessellate_edge = |edge: &CompressedEdge<C>| {
let curve = &edge.curve;
CompressedEdge {
vertices: edge.vertices,
curve: PolylineCurve::from_curve(curve, curve.range_tuple(), tol),
}
};
#[cfg(not(target_arch = "wasm32"))]
let edges: Vec<_> = shell.edges.par_iter().map(tessellate_edge).collect();
#[cfg(target_arch = "wasm32")]
let edges: Vec<_> = shell.edges.iter().map(tessellate_edge).collect();
let tessellate_face = |face: &CompressedFace<S>| {
let boundaries = face.boundaries.clone();
let surface = &face.surface;
let create_edge = |edge_idx: &CompressedEdgeIndex| match edge_idx.orientation {
true => Some(edges.get(edge_idx.index)?.curve.clone()),
false => Some(edges.get(edge_idx.index)?.curve.inverse()),
};
let create_boundary = |wire: &Vec<CompressedEdgeIndex>| {
let wire_iter = wire.iter().filter_map(create_edge);
PolyBoundaryPiece::try_new(surface, wire_iter, &sp)
};
let preboundary: Option<Vec<_>> = boundaries.iter().map(create_boundary).collect();
let polygon: Option<PolygonMesh> = (|| {
let boundary = PolyBoundary::new(preboundary?, &surface, tol);
Some(trimming_tessellation(&surface, &boundary, tol))
})();
CompressedFace {
boundaries,
orientation: face.orientation,
surface: polygon,
}
};
#[cfg(not(target_arch = "wasm32"))]
let faces = shell.faces.par_iter().map(tessellate_face).collect();
#[cfg(target_arch = "wasm32")]
let faces = shell.faces.iter().map(tessellate_face).collect();
MeshedCShell {
vertices,
edges,
faces,
}
}
fn shell_create_polygon<S: PreMeshableSurface>(
surface: &S,
wires: Vec<Wire<Point3, PolylineCurve>>,
orientation: bool,
tol: f64,
sp: impl SP<S>,
) -> Face<Point3, PolylineCurve, Option<PolygonMesh>> {
let preboundary = wires
.iter()
.map(|wire: &Wire<_, _>| {
let wire_iter = wire.iter().map(Edge::oriented_curve);
PolyBoundaryPiece::try_new(surface, wire_iter, &sp)
})
.collect::<Option<Vec<_>>>();
let polygon: Option<PolygonMesh> = (|| {
let boundary = PolyBoundary::new(preboundary?, &surface, tol);
Some(trimming_tessellation(surface, &boundary, tol))
})();
let mut new_face = Face::debug_new(wires, polygon);
if !orientation {
new_face.invert();
}
new_face
}
#[derive(Debug, Default, Clone)]
struct PolyBoundaryPiece(Vec<Point2>);
impl PolyBoundaryPiece {
fn try_new<S: PreMeshableSurface>(
surface: &S,
wire: impl Iterator<Item = PolylineCurve>,
sp: impl SP<S>,
) -> Option<Self> {
let (up, vp) = (surface.u_period(), surface.v_period());
let (urange, vrange) = surface.try_range_tuple();
let mut bdry3d: Vec<Point3> = wire
.flat_map(|poly_edge| {
let n = poly_edge.len() - 1;
poly_edge.into_iter().take(n)
})
.collect();
bdry3d.push(bdry3d[0]);
let mut previous = None;
let mut vec = bdry3d
.into_iter()
.flat_map(|pt| {
let (mut u, mut v) = match sp(surface, pt, previous) {
Some(hint) => hint,
None => return vec![None],
};
if let (Some(up), Some((u0, _))) = (up, previous) {
u = get_mindiff(u, u0, up);
}
if let (Some(vp), Some((_, v0))) = (vp, previous) {
v = get_mindiff(v, v0, vp);
}
let res = (|| {
if let Some((u0, v0)) = previous {
if !u0.near(&u) && surface.uder(u0, v0).so_small() {
return vec![Some(Point2::new(u, v0)), Some(Point2::new(u, v))];
} else if !v0.near(&v) && surface.vder(u0, v0).so_small() {
return vec![Some(Point2::new(u0, v)), Some(Point2::new(u, v))];
}
}
vec![Some(Point2::new(u, v))]
})();
previous = Some((u, v));
res
})
.collect::<Option<Vec<Point2>>>()?;
let grav = vec.iter().fold(Point2::origin(), |g, p| g + p.to_vec()) / vec.len() as f64;
if let (Some(up), Some((u0, _))) = (up, urange) {
let quot = f64::floor((grav.x - u0) / up);
vec.iter_mut().for_each(|p| p.x -= quot * up);
}
if let (Some(vp), Some((v0, _))) = (vp, vrange) {
let quot = f64::floor((grav.y - v0) / vp);
vec.iter_mut().for_each(|p| p.y -= quot * vp);
}
let last = *vec.last().unwrap();
if !vec[0].near(&last) {
let Point2 { x: u0, y: v0 } = last;
if surface.uder(u0, v0).so_small() || surface.vder(u0, v0).so_small() {
vec.push(vec[0]);
}
}
Some(Self(vec))
}
}
fn abs_diff(previous: f64) -> impl Fn(&f64, &f64) -> std::cmp::Ordering {
let f = move |x: &f64| f64::abs(x - previous);
move |x: &f64, y: &f64| f(x).partial_cmp(&f(y)).unwrap()
}
fn get_mindiff(u: f64, u0: f64, up: f64) -> f64 {
let closure = |i| u + i as f64 * up;
(-2..=2).map(closure).min_by(abs_diff(u0)).unwrap()
}
#[derive(Debug, Default, Clone)]
struct PolyBoundary(Vec<Vec<Point2>>);
fn normalize_range(curve: &mut Vec<Point2>, compidx: usize, (u0, u1): (f64, f64)) {
let p = curve[0];
let q = curve[curve.len() - 1];
let tmp = f64::min(p[compidx], q[compidx]) + TOLERANCE;
let del = f64::floor((tmp - u0) / (u1 - u0)) * (u1 - u0);
curve.iter_mut().for_each(|p| p[compidx] -= del);
let Some(i) = curve
.iter()
.position(|p| (curve[0][compidx] - u1) * (p[compidx] - u1) < 0.0)
else {
return;
};
let mut curve1 = curve.split_off(i + 1);
curve1.pop();
curve1.insert(0, curve[i]);
match curve[0][compidx] < curve[curve.len() - 1][compidx] {
true => curve1.iter_mut(),
false => curve.iter_mut(),
}
.for_each(|p| p[compidx] -= u1 - u0);
curve1.append(curve);
*curve = curve1;
}
fn loop_orientation(curve: &[Point2]) -> bool {
curve
.iter()
.circular_tuple_windows()
.fold(0.0, |sum, (p, q)| sum + (q.x + p.x) * (q.y - p.y))
> 0.0
}
impl PolyBoundary {
fn new(pieces: Vec<PolyBoundaryPiece>, surface: &impl PreMeshableSurface, tol: f64) -> Self {
let (mut closed, mut open) = (Vec::new(), Vec::new());
pieces.into_iter().for_each(|PolyBoundaryPiece(mut vec)| {
match vec[0].distance(vec[vec.len() - 1]) < 1.0e-3 {
true => {
vec.pop();
closed.push(vec)
}
false => open.push(vec),
}
});
fn connect_edges(vecs: impl IntoIterator<Item = Vec<Point2>>) -> Vec<Point2> {
let closure = |vec: Vec<Point2>| {
let len = vec.len();
vec.into_iter().take(len - 1)
};
vecs.into_iter().flat_map(closure).collect()
}
match open.len() {
1 => {
let mut curve = open.pop().unwrap();
let p = curve[0];
let q = curve[curve.len() - 1];
if let (Some((u0, u1)), Some((v0, v1))) = surface.try_range_tuple() {
if p.x < q.x - TOLERANCE {
normalize_range(&mut curve, 0, (u0, u1));
let p = curve[0];
let q = curve[curve.len() - 1];
let x = Point2::new(u0, v1);
let y = Point2::new(u1, v1);
let vec0 = polyline_on_surface(surface, q, y, tol);
let vec1 = polyline_on_surface(surface, y, x, tol);
let vec2 = polyline_on_surface(surface, x, p, tol);
closed.push(connect_edges([vec0, vec1, vec2, curve]));
} else if q.x < p.x - TOLERANCE {
normalize_range(&mut curve, 0, (u0, u1));
let p = curve[0];
let q = curve[curve.len() - 1];
let x = Point2::new(u1, v0);
let y = Point2::new(u0, v0);
let vec0 = polyline_on_surface(surface, q, y, tol);
let vec1 = polyline_on_surface(surface, y, x, tol);
let vec2 = polyline_on_surface(surface, x, p, tol);
closed.push(connect_edges([vec0, vec1, vec2, curve]));
} else if p.y < q.y - TOLERANCE {
normalize_range(&mut curve, 1, (v0, v1));
let p = curve[0];
let q = curve[curve.len() - 1];
let x = Point2::new(u0, v0);
let y = Point2::new(u0, v1);
let vec0 = polyline_on_surface(surface, q, y, tol);
let vec1 = polyline_on_surface(surface, y, x, tol);
let vec2 = polyline_on_surface(surface, x, p, tol);
closed.push(connect_edges([vec0, vec1, vec2, curve]));
} else if q.y < p.y - TOLERANCE {
normalize_range(&mut curve, 1, (v0, v1));
let p = curve[0];
let q = curve[curve.len() - 1];
let x = Point2::new(u1, v1);
let y = Point2::new(u1, v0);
let vec0 = polyline_on_surface(surface, q, y, tol);
let vec1 = polyline_on_surface(surface, y, x, tol);
let vec2 = polyline_on_surface(surface, x, p, tol);
closed.push(connect_edges([vec0, vec1, vec2, curve]));
}
}
}
2 => {
let mut curve1 = open.pop().unwrap();
let mut curve0 = open.pop().unwrap();
fn end_pts<T: Copy>(vec: &[T]) -> (T, T) { (vec[0], vec[vec.len() - 1]) }
let ((p0, p1), (q0, q1)) = (end_pts(&curve0), end_pts(&curve1));
if !p0.x.near(&p1.x) && !q0.x.near(&q1.x) {
if let (Some(urange), _) = surface.try_range_tuple() {
normalize_range(&mut curve0, 0, urange);
normalize_range(&mut curve1, 0, urange);
}
} else if !p0.y.near(&p1.y) && !q0.y.near(&q1.y) {
if let (_, Some(vrange)) = surface.try_range_tuple() {
normalize_range(&mut curve0, 1, vrange);
normalize_range(&mut curve1, 1, vrange);
}
}
let ((p0, p1), (q0, q1)) = (end_pts(&curve0), end_pts(&curve1));
let vec0 = polyline_on_surface(surface, p1, q0, tol);
let vec1 = polyline_on_surface(surface, q1, p0, tol);
closed.push(connect_edges([curve0, vec0, curve1, vec1]));
}
_ => {}
}
if !closed.iter().any(|curve| loop_orientation(curve)) {
if let (Some((u0, u1)), Some((v0, v1))) = surface.try_range_tuple() {
let p = [
Point2::new(u0, v0),
Point2::new(u1, v0),
Point2::new(u1, v1),
Point2::new(u0, v1),
];
let vec0 = polyline_on_surface(surface, p[0], p[1], tol);
let vec1 = polyline_on_surface(surface, p[1], p[2], tol);
let vec2 = polyline_on_surface(surface, p[2], p[3], tol);
let vec3 = polyline_on_surface(surface, p[3], p[0], tol);
closed.push(connect_edges([vec0, vec1, vec2, vec3]));
}
}
Self(closed)
}
/// whether `c` is included in the domain with boundary = `self`.
fn include(&self, c: Point2) -> bool {
let t = 2.0 * std::f64::consts::PI * HashGen::hash1(c);
let r = Vector2::new(f64::cos(t), f64::sin(t));
self.0
.iter()
.flat_map(|vec| vec.iter().circular_tuple_windows())
.try_fold(0_i32, move |counter, (p0, p1)| {
let a = p0 - c;
let b = p1 - c;
let s0 = r.x * a.y - r.y * a.x; // v times a
let s1 = r.x * b.y - r.y * b.x; // v times b
let s2 = a.x * b.y - a.y * b.x; // a times b
let x = s2 / (s1 - s0);
if x.so_small() && s0 * s1 < 0.0 {
None
} else if x > 0.0 && s0 <= 0.0 && s1 > 0.0 {
Some(counter + 1)
} else if x > 0.0 && s0 >= 0.0 && s1 < 0.0 {
Some(counter - 1)
} else {
Some(counter)
}
})
.map(|counter| counter > 0)
.unwrap_or(false)
}
/// Inserts points and adds constraint into triangulation.
fn insert_to(&self, triangulation: &mut Cdt) {
let poly2tri: Vec<_> = self
.0
.iter()
.flatten()
.map(|pt| {
let p = [spade_round(pt.x), spade_round(pt.y)];
triangulation.insert(SPoint2::from(p)).ok()
})
.collect();
let mut prev: Option<usize> = None;
let mut counter = 0;
self.0
.iter()
.map(Vec::len)
.flat_map(|len| {
let range = counter..counter + len;
counter += len;
range.circular_tuple_windows()
})
.for_each(|(i, j)| {
let Some(vj) = poly2tri[j] else { return };
if let Some(p) = prev {
let Some(v) = poly2tri[p] else { return };
if triangulation.can_add_constraint(v, vj) {
triangulation.add_constraint(v, vj);
prev = None;
}
} else {
let Some(vi) = poly2tri[i] else { return };
if triangulation.can_add_constraint(vi, vj) {
triangulation.add_constraint(vi, vj);
} else {
prev = Some(i);
}
}
});
}
}
fn spade_round(x: f64) -> f64 {
match f64::abs(x) < MIN_ALLOWED_VALUE {
true => 0.0,
false => x,
}
}
/// Tessellates one surface trimmed by polyline.
fn trimming_tessellation<S>(surface: &S, polyboundary: &PolyBoundary, tol: f64) -> PolygonMesh
where S: PreMeshableSurface {
let mut triangulation = Cdt::new();
polyboundary.insert_to(&mut triangulation);
insert_surface(&mut triangulation, surface, polyboundary, tol);
let mut mesh = triangulation_into_polymesh(
triangulation.vertices(),
triangulation.inner_faces(),
surface,
polyboundary,
);
mesh.make_face_compatible_to_normal();
mesh
}
/// Inserts parameter divisions into triangulation.
fn insert_surface(
triangulation: &mut Cdt,
surface: impl PreMeshableSurface,
polyline: &PolyBoundary,
tol: f64,
) {
let bdb: BoundingBox<Point2> = polyline.0.iter().flatten().collect();
let range = ((bdb.min()[0], bdb.max()[0]), (bdb.min()[1], bdb.max()[1]));
let (udiv, vdiv) = surface.parameter_division(range, tol);
let insert_res: Vec<Vec<Option<_>>> = udiv
.into_iter()
.map(|u| {
vdiv.iter()
.map(|v| match polyline.include(Point2::new(u, *v)) {
true => triangulation.insert(SPoint2::new(u, *v)).ok(),
false => None,
})
.collect()
})
.collect();
insert_res.windows(2).for_each(|vec| {
vec[0].windows(2).zip(&vec[1]).for_each(|(a, z)| {
if let Some(x) = a[0] {
if let Some(y) = a[1] {
if triangulation.can_add_constraint(x, y) {
triangulation.add_constraint(x, y);
}
}
if let Some(z) = z {
if triangulation.can_add_constraint(x, *z) {
triangulation.add_constraint(x, *z);
}
}
}
});
let idx = vec[0].len() - 1;
if let (Some(x), Some(y)) = (vec[0][idx], vec[1][idx]) {
if triangulation.can_add_constraint(x, y) {
triangulation.add_constraint(x, y);
}
}
});
}
/// Converts triangulation into `PolygonMesh`.
fn triangulation_into_polymesh<'a>(
vertices: VertexIterator<'a, SPoint2, (), CdtEdge<()>, ()>,
triangles: InnerFaceIterator<'a, SPoint2, (), CdtEdge<()>, ()>,
surface: &impl ParametricSurface3D,
polyline: &PolyBoundary,
) -> PolygonMesh {
let mut positions = Vec::<Point3>::new();
let mut uv_coords = Vec::<Vector2>::new();
let mut normals = Vec::<Vector3>::new();
let vmap: HashMap<_, _> = vertices
.enumerate()
.map(|(i, v)| {
let p = *v.as_ref();
let uv = Vector2::new(p.x, p.y);
positions.push(surface.subs(uv[0], uv[1]));
uv_coords.push(uv);
normals.push(surface.normal(uv[0], uv[1]));
(v.fix(), i)
})
.collect();
let tri_faces: Vec<[StandardVertex; 3]> = triangles
.map(|tri| tri.vertices())
.filter(|tri| {
fn sp2cg(p: SPoint2) -> Point2 { Point2::new(p.x, p.y) }
let tri = array![i => sp2cg(*tri[i].as_ref()); 3];
let (a, b) = (tri[1] - tri[0], tri[2] - tri[0]);
let c = tri[0] + (a + b) / 3.0;
let area = a.x * b.y - a.y * b.x;
polyline.include(c) && !area.so_small2()
})
.map(|tri| {
let idcs = array![i => vmap[&tri[i].fix()]; 3];
array![i => [idcs[i], idcs[i], idcs[i]].into(); 3]
})
.collect();
PolygonMesh::debug_new(
StandardAttributes {
positions,
uv_coords,
normals,
},
Faces::from_tri_and_quad_faces(tri_faces, Vec::new()),
)
}
fn polyline_on_surface(
surface: impl PreMeshableSurface,
p: Point2,
q: Point2,
tol: f64,
) -> Vec<Point2> {
use truck_geometry::prelude::*;
let line = Line(p, q);
let pcurve = PCurve::new(line, surface);
let (vec, _) = pcurve.parameter_division(pcurve.range_tuple(), tol);
vec.into_iter().map(|t| line.subs(t)).collect()
}
#[test]
#[ignore]
#[cfg(not(target_arch = "wasm32"))]
fn par_bench() {
use std::time::Instant;
use truck_modeling::*;
const JSON: &str = include_str!(concat!(
env!("CARGO_MANIFEST_DIR"),
"/../resources/shape/bottle.json"
));
let solid: Solid = serde_json::from_str(JSON).unwrap();
let shell = solid.into_boundaries().pop().unwrap();
let instant = Instant::now();
(0..100).for_each(|_| {
let _shell = shell_tessellation(&shell, 0.01, by_search_parameter);
});
println!("{}ms", instant.elapsed().as_millis());
let instant = Instant::now();
(0..100).for_each(|_| {
let _shell = shell_tessellation_single_thread(&shell, 0.01, by_search_parameter);
});
println!("{}ms", instant.elapsed().as_millis());
}
+525
View File
@@ -0,0 +1,525 @@
use crate::*;
use rustc_hash::FxHashMap as HashMap;
use std::convert::identity;
use truck_base::tolerance::TOLERANCE;
use truck_topology::{compress::*, Vertex, *};
pub use vtkio;
use vtkio::model::{Attributes, *};
/// Trait for convert to [`DataSet`]
pub trait ToDataSet {
/// Converts to [`DataSet`]
/// # Remarks
/// The `DataSet` is generated assuming output in XML format.
fn to_data_set(&self) -> DataSet;
}
fn to_vertex_numbers(faces: &Faces<usize>) -> VertexNumbers {
let connectivity = faces
.face_iter()
.flat_map(|face| face.iter().map(|idx| *idx as u64))
.collect::<Vec<_>>();
let mut offset: u64 = 0;
let offsets = faces
.face_iter()
.map(|face| {
offset += face.len() as u64;
offset
})
.collect::<Vec<_>>();
VertexNumbers::XML {
connectivity,
offsets,
}
}
impl ToDataSet for PolygonMesh<usize, Vec<Point3>> {
fn to_data_set(&self) -> DataSet {
let flatten_points = self
.attributes()
.iter()
.flat_map(|p| Into::<[f64; 3]>::into(*p))
.collect::<Vec<_>>();
DataSet::PolyData {
meta: None,
pieces: vec![Piece::Inline(Box::new(PolyDataPiece {
points: IOBuffer::F64(flatten_points),
polys: Some(to_vertex_numbers(self.faces())),
..Default::default()
}))],
}
}
}
impl ToDataSet for PolygonMesh<usize, Vec<(Point3, Vector2)>> {
fn to_data_set(&self) -> DataSet {
let flatten_points = self
.attributes()
.iter()
.flat_map(|(p, _)| Into::<[f64; 3]>::into(*p))
.collect::<Vec<_>>();
let flatten_uvs = self
.attributes()
.iter()
.flat_map(|(_, uv)| Into::<[f64; 2]>::into(*uv))
.collect::<Vec<_>>();
let uvs = Attribute::DataArray(DataArray {
name: "TCoords".to_owned(),
elem: ElementType::TCoords(2),
data: IOBuffer::F64(flatten_uvs),
});
DataSet::PolyData {
meta: None,
pieces: vec![Piece::Inline(Box::new(PolyDataPiece {
points: IOBuffer::F64(flatten_points),
polys: Some(to_vertex_numbers(self.faces())),
data: Attributes {
point: vec![uvs],
..Default::default()
},
..Default::default()
}))],
}
}
}
impl ToDataSet for PolygonMesh<usize, Vec<(Point3, Vector3)>> {
fn to_data_set(&self) -> DataSet {
let flatten_points = self
.attributes()
.iter()
.flat_map(|(p, _)| Into::<[f64; 3]>::into(*p))
.collect::<Vec<_>>();
let flatten_normals = self
.attributes()
.iter()
.flat_map(|(_, n)| Into::<[f64; 3]>::into(*n))
.collect::<Vec<_>>();
let normals = Attribute::DataArray(DataArray {
name: "Normals".to_owned(),
elem: ElementType::Normals,
data: IOBuffer::F64(flatten_normals),
});
DataSet::PolyData {
meta: None,
pieces: vec![Piece::Inline(Box::new(PolyDataPiece {
points: IOBuffer::F64(flatten_points),
polys: Some(to_vertex_numbers(self.faces())),
data: Attributes {
point: vec![normals],
..Default::default()
},
..Default::default()
}))],
}
}
}
impl ToDataSet for PolygonMesh<usize, Vec<StandardAttribute>> {
fn to_data_set(&self) -> DataSet {
let flatten_points = self
.attributes()
.iter()
.flat_map(|attr| Into::<[f64; 3]>::into(attr.position))
.collect::<Vec<_>>();
let flatten_uvs = self
.attributes()
.iter()
.flat_map(|attr| match attr.uv_coord {
Some(uv) => Into::<[f64; 2]>::into(uv),
None => [f64::NAN; 2],
})
.collect::<Vec<_>>();
let uvs = Attribute::DataArray(DataArray {
name: "TCoords".to_owned(),
elem: ElementType::TCoords(2),
data: IOBuffer::F64(flatten_uvs),
});
let flatten_normals = self
.attributes()
.iter()
.flat_map(|attr| match attr.normal {
Some(normal) => Into::<[f64; 3]>::into(normal),
None => [f64::NAN; 3],
})
.collect::<Vec<_>>();
let normals = Attribute::DataArray(DataArray {
name: "Normals".to_owned(),
elem: ElementType::Normals,
data: IOBuffer::F64(flatten_normals),
});
DataSet::PolyData {
meta: None,
pieces: vec![Piece::Inline(Box::new(PolyDataPiece {
points: IOBuffer::F64(flatten_points),
polys: Some(to_vertex_numbers(self.faces())),
data: Attributes {
point: vec![uvs, normals],
..Default::default()
},
..Default::default()
}))],
}
}
}
impl ToDataSet for PolygonMesh {
fn to_data_set(&self) -> DataSet {
let attrs = self.attributes();
match (attrs.uv_coords().is_empty(), attrs.normals().is_empty()) {
(true, true) => self.to_positions_mesh().to_data_set(),
(false, true) => self
.expands(|attr| {
let uv = match attr.uv_coord {
Some(uv) => uv,
None => Vector2::new(f64::NAN, f64::NAN),
};
(attr.position, uv)
})
.to_data_set(),
(true, false) => self
.expands(|attr| {
let normal = match attr.normal {
Some(normal) => normal,
None => Vector3::new(f64::NAN, f64::NAN, f64::NAN),
};
(attr.position, normal)
})
.to_data_set(),
(false, false) => self.expands(identity).to_data_set(),
}
}
}
fn hash_point(p: Point3) -> Option<[i64; 3]> {
(p / (TOLERANCE * 50.0) - Vector3::new(0.5, 0.5, 0.5))
.cast()
.map(Into::into)
}
impl ToDataSet for Vertex<Point3> {
fn to_data_set(&self) -> DataSet {
DataSet::UnstructuredGrid {
meta: None,
pieces: vec![Piece::Inline(Box::new(UnstructuredGridPiece {
points: IOBuffer::F64(Into::<[f64; 3]>::into(self.point()).to_vec()),
cells: Cells {
cell_verts: VertexNumbers::XML {
connectivity: vec![0],
offsets: vec![1],
},
types: vec![CellType::Vertex],
},
data: Default::default(),
}))],
}
}
}
impl ToDataSet for Edge<Point3, PolylineCurve<Point3>> {
fn to_data_set(&self) -> DataSet {
let curve = self.oriented_curve();
let points = curve
.0
.iter()
.copied()
.flat_map(Into::<[f64; 3]>::into)
.collect::<Vec<_>>();
let len = curve.0.len() as u64;
let connectivity = (0..len).chain(vec![0, len - 1]).collect();
DataSet::UnstructuredGrid {
meta: None,
pieces: vec![Piece::Inline(Box::new(UnstructuredGridPiece {
points: IOBuffer::F64(points),
cells: Cells {
cell_verts: VertexNumbers::XML {
connectivity,
offsets: vec![len, len + 1, len + 2],
},
types: vec![CellType::PolyLine, CellType::Vertex, CellType::Vertex],
},
data: Default::default(),
}))],
}
}
}
impl ToDataSet for Wire<Point3, PolylineCurve<Point3>> {
fn to_data_set(&self) -> DataSet {
let mut points = Vec::<[f64; 3]>::new();
let mut connectivity = Vec::<u64>::new();
let mut offsets = Vec::<u64>::new();
let mut types = Vec::<CellType>::new();
let mut vmap = truck_base::entry_map::FxEntryMap::new(
|v: &Vertex<Point3>| v.id(),
|v: &Vertex<Point3>| {
points.push(v.point().into());
points.len() as u64 - 1
},
);
let tmp: Vec<_> = self.vertex_iter().collect();
tmp.iter().for_each(|v| {
vmap.entry_or_insert(v);
});
let vmap: HashMap<_, _> = vmap.into();
self.edge_iter().for_each(|edge| {
let pts = edge.oriented_curve().0;
connectivity.push(*vmap.get(&edge.front().id()).unwrap());
connectivity.extend((1..pts.len() - 1).map(|i| (points.len() + i - 1) as u64));
connectivity.push(*vmap.get(&edge.back().id()).unwrap());
offsets.push(connectivity.len() as u64);
points.extend(
pts[1..pts.len() - 1]
.iter()
.copied()
.map(Into::<[f64; 3]>::into),
);
types.push(CellType::PolyLine);
});
self.vertex_iter().enumerate().for_each(|(idx, _)| {
connectivity.push(idx as u64);
offsets.push(connectivity.len() as u64);
types.push(CellType::Vertex);
});
DataSet::UnstructuredGrid {
meta: None,
pieces: vec![Piece::Inline(Box::new(UnstructuredGridPiece {
points: IOBuffer::F64(points.concat().to_vec()),
cells: Cells {
cell_verts: VertexNumbers::XML {
connectivity,
offsets,
},
types,
},
data: Default::default(),
}))],
}
}
}
fn face_piece(
face: &Face<Point3, PolylineCurve<Point3>, PolygonMesh>,
) -> Piece<UnstructuredGridPiece> {
let polygon = face.oriented_surface().expands(identity);
let mut length = 0;
let map: HashMap<Option<[i64; 3]>, usize> = polygon
.attributes()
.iter()
.map(move |attr| {
length += 1;
(hash_point(attr.position), length - 1)
})
.collect();
let flatten_points = polygon
.attributes()
.iter()
.flat_map(|attr| Into::<[f64; 3]>::into(attr.position))
.collect::<Vec<_>>();
let flatten_uvs = polygon
.attributes()
.iter()
.flat_map(|attr| match attr.uv_coord {
Some(uv) => Into::<[f64; 2]>::into(uv),
None => [f64::NAN; 2],
})
.collect::<Vec<_>>();
let uvs = Attribute::DataArray(DataArray {
name: "TCoords".to_owned(),
elem: ElementType::TCoords(2),
data: IOBuffer::F64(flatten_uvs),
});
let flatten_normals = polygon
.attributes()
.iter()
.flat_map(|attr| match attr.normal {
Some(normal) => Into::<[f64; 3]>::into(normal),
None => [f64::NAN; 3],
})
.collect::<Vec<_>>();
let normals = Attribute::DataArray(DataArray {
name: "Normals".to_owned(),
elem: ElementType::Normals,
data: IOBuffer::F64(flatten_normals),
});
let mut cell_verts = to_vertex_numbers(polygon.faces());
let mut types = vec![CellType::Polygon; cell_verts.num_cells()];
face.boundaries().into_iter().flatten().for_each(|edge| {
types.push(CellType::PolyLine);
let curve = edge.oriented_curve().0;
let polyline = curve
.into_iter()
.filter_map(|p| map.get(&hash_point(p)).map(|i| *i as u64));
if let VertexNumbers::XML {
connectivity,
offsets,
} = &mut cell_verts
{
connectivity.extend(polyline);
offsets.push(connectivity.len() as u64);
}
});
face.vertex_iter().try_for_each(|v| {
types.push(CellType::Vertex);
let idx = map.get(&hash_point(v.point())).map(|i| *i as u64)?;
if let VertexNumbers::XML {
connectivity,
offsets,
} = &mut cell_verts
{
connectivity.push(idx);
offsets.push(connectivity.len() as u64)
}
Some(())
});
Piece::Inline(Box::new(UnstructuredGridPiece {
points: IOBuffer::F64(flatten_points),
cells: Cells { cell_verts, types },
data: Attributes {
point: vec![uvs, normals],
..Default::default()
},
}))
}
impl ToDataSet for Face<Point3, PolylineCurve<Point3>, PolygonMesh> {
fn to_data_set(&self) -> DataSet {
DataSet::UnstructuredGrid {
meta: None,
pieces: vec![face_piece(self)],
}
}
}
impl ToDataSet for Shell<Point3, PolylineCurve<Point3>, PolygonMesh> {
fn to_data_set(&self) -> DataSet {
DataSet::UnstructuredGrid {
meta: None,
pieces: self.face_iter().map(face_piece).collect(),
}
}
}
fn cface_piece(
vertices: &[Point3],
edges: &[CompressedEdge<PolylineCurve<Point3>>],
face: &CompressedFace<PolygonMesh>,
) -> Option<Piece<UnstructuredGridPiece>> {
let polygon = match face.orientation {
true => face.surface.expands(identity),
false => face.surface.inverse().expands(identity),
};
let mut length = 0;
let map: HashMap<Option<[i64; 3]>, usize> = polygon
.attributes()
.iter()
.map(move |attr| {
length += 1;
(hash_point(attr.position), length - 1)
})
.collect();
let flatten_points = polygon
.attributes()
.iter()
.flat_map::<[f64; 3], _>(|attr| attr.position.into())
.collect::<Vec<_>>();
let flatten_uvs = polygon
.attributes()
.iter()
.flat_map(|attr| match attr.uv_coord {
Some(uv) => uv.into(),
None => [f64::NAN; 2],
})
.collect::<Vec<_>>();
let uvs = Attribute::DataArray(DataArray {
name: "TCoords".to_owned(),
elem: ElementType::TCoords(2),
data: IOBuffer::F64(flatten_uvs),
});
let flatten_normals = polygon
.attributes()
.iter()
.flat_map(|attr| match attr.normal {
Some(normal) => normal.into(),
None => [f64::NAN; 3],
})
.collect::<Vec<_>>();
let normals = Attribute::DataArray(DataArray {
name: "Normals".to_owned(),
elem: ElementType::Normals,
data: IOBuffer::F64(flatten_normals),
});
let mut cell_verts = to_vertex_numbers(polygon.faces());
let mut types = vec![CellType::Polygon; cell_verts.num_cells()];
let edge_indices: Vec<_> = face
.boundaries
.iter()
.flat_map::<Box<dyn Iterator<Item = &CompressedEdgeIndex>>, _>(|boundary| {
match face.orientation {
true => Box::new(boundary.iter()),
false => Box::new(boundary.iter().rev()),
}
})
.collect();
edge_indices.iter().try_for_each(|edge| {
let curve = match face.orientation == edge.orientation {
true => edges[edge.index].curve.clone(),
false => edges[edge.index].curve.inverse(),
};
let polyline = curve
.iter()
.filter_map(|p| map.get(&hash_point(*p)).map(|i| *i as u64))
.collect::<Vec<_>>();
if polyline.is_empty() {
return None;
}
types.push(CellType::PolyLine);
if let VertexNumbers::XML {
connectivity,
offsets,
} = &mut cell_verts
{
connectivity.extend(polyline);
offsets.push(connectivity.len() as u64);
}
Some(())
})?;
edge_indices.into_iter().try_for_each(|edge| {
let v = match face.orientation == edge.orientation {
true => edges[edge.index].vertices.0,
false => edges[edge.index].vertices.1,
};
let idx = map.get(&hash_point(vertices[v])).map(|i| *i as u64)?;
types.push(CellType::Vertex);
if let VertexNumbers::XML {
connectivity,
offsets,
} = &mut cell_verts
{
connectivity.push(idx);
offsets.push(connectivity.len() as u64)
}
Some(())
})?;
Some(Piece::Inline(Box::new(UnstructuredGridPiece {
points: IOBuffer::F64(flatten_points),
cells: Cells { cell_verts, types },
data: Attributes {
point: vec![uvs, normals],
..Default::default()
},
})))
}
impl ToDataSet for CompressedShell<Point3, PolylineCurve<Point3>, PolygonMesh> {
fn to_data_set(&self) -> DataSet {
DataSet::UnstructuredGrid {
meta: None,
pieces: self
.faces
.iter()
.filter_map(|face| cface_piece(&self.vertices, &self.edges, face))
.collect(),
}
}
}