use crate::projection::{self, SurfaceProjection};
use crate::surface_math::analytic;
use crate::tol;
use crate::{Frame3D, Point3D, Tolerance, Vector3D};
use std::fmt;
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum ConeConstructionError {
NegativeRadius,
BadAngle,
NullAxis,
NullAngle,
NullNormal,
}
impl fmt::Display for ConeConstructionError {
fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
let message = match self {
ConeConstructionError::NegativeRadius => "radius is negative",
ConeConstructionError::BadAngle => "semi-angle is out of range",
ConeConstructionError::NullAxis => "the two points are confused",
ConeConstructionError::NullAngle => "semi-angle is null or a right angle",
ConeConstructionError::NullNormal => "axis direction has zero length",
};
f.write_str(message)
}
}
impl std::error::Error for ConeConstructionError {}
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct Cone {
frame: Frame3D,
semi_angle: f64,
ref_radius: f64,
}
impl Cone {
pub fn from_frame(
frame: Frame3D,
semi_angle: f64,
ref_radius: f64,
) -> Result<Cone, ConeConstructionError> {
if !is_valid_semi_angle(semi_angle) {
return Err(ConeConstructionError::BadAngle);
}
if ref_radius < 0.0 {
return Err(ConeConstructionError::NegativeRadius);
}
Ok(Cone {
frame,
semi_angle,
ref_radius,
})
}
pub fn new(
center: Point3D,
axis_direction: Vector3D,
semi_angle: f64,
ref_radius: f64,
) -> Result<Cone, ConeConstructionError> {
let frame = Frame3D::from_z(center, axis_direction)
.map_err(|_| ConeConstructionError::NullNormal)?;
Cone::from_frame(frame, semi_angle, ref_radius)
}
pub fn from_two_points_and_radii(
p1: Point3D,
p2: Point3D,
r1: f64,
r2: f64,
) -> Result<Cone, ConeConstructionError> {
let delta = p2 - p1;
let dist = delta.magnitude();
if dist < tol::CONFUSION {
return Err(ConeConstructionError::NullAxis);
}
if r1 < 0.0 || r2 < 0.0 {
return Err(ConeConstructionError::NegativeRadius);
}
let semi_angle = ((r2 - r1) / dist).atan();
if !is_valid_semi_angle(semi_angle) {
return Err(ConeConstructionError::NullAngle);
}
let z_dir = delta * (1.0 / dist);
let frame = arbitrary_perpendicular_frame(p1, z_dir);
Ok(Cone {
frame,
semi_angle,
ref_radius: r1,
})
}
pub fn frame(&self) -> Frame3D {
self.frame
}
pub fn semi_angle(&self) -> f64 {
self.semi_angle
}
pub fn ref_radius(&self) -> f64 {
self.ref_radius
}
pub fn apex(&self) -> Point3D {
let v = -self.ref_radius / self.semi_angle.tan();
self.frame.point_at(0.0, 0.0, v)
}
pub fn eval_point(&self, u: f64, v: f64) -> Point3D {
analytic::cone_d0(&self.frame, self.ref_radius, self.semi_angle, u, v)
}
pub fn eval_points(&self, uvs: &[(f64, f64)]) -> Vec<Point3D> {
uvs.iter().map(|&(u, v)| self.eval_point(u, v)).collect()
}
pub fn eval_derivative(&self, u: f64, v: f64, du: u32, dv: u32) -> Vector3D {
match du + dv {
0 => panic!(
"eval_derivative: du + dv must be >= 1 (use eval_point for the (0, 0) order)"
),
1..=2 => analytic::cone_derivative(
&self.frame,
self.ref_radius,
self.semi_angle,
u,
v,
du,
dv,
),
_ => panic!(
"eval_derivative: order du={du}, dv={dv} is not supported (du + dv must be <= 2)"
),
}
}
pub fn parameters_of(&self, point: Point3D) -> (f64, f64) {
analytic::cone_parameters(&self.frame, self.ref_radius, self.semi_angle, point)
}
pub fn contains(&self, point: Point3D, tol: Tolerance) -> bool {
let (u, v) = self.parameters_of(point);
self.eval_point(u, v).distance(point) <= tol.confusion
}
pub fn project_point(&self, point: Point3D, tol: Tolerance) -> SurfaceProjection {
let (u, v) = self.parameters_of(point);
let distance = self.eval_point(u, v).distance(point);
SurfaceProjection {
u,
v,
distance: projection::snap_distance(distance, tol.confusion),
}
}
pub fn project_points(&self, points: &[Point3D], tol: Tolerance) -> Vec<SurfaceProjection> {
points.iter().map(|&p| self.project_point(p, tol)).collect()
}
}
fn is_valid_semi_angle(semi_angle: f64) -> bool {
let angle = semi_angle.abs();
angle >= tol::ANGULAR && (std::f64::consts::FRAC_PI_2 - angle) >= tol::ANGULAR
}
fn arbitrary_perpendicular_frame(origin: Point3D, axis_direction: Vector3D) -> Frame3D {
let z = axis_direction;
let (ax, ay, az) = (z.x.abs(), z.y.abs(), z.z.abs());
let x_dir = if ax <= ay && ax <= az {
Vector3D::new(0.0, -z.z, z.y)
} else if ay <= ax && ay <= az {
Vector3D::new(-z.z, 0.0, z.x)
} else {
Vector3D::new(-z.y, z.x, 0.0)
}
.normalized()
.expect("z is a unit vector and nonzero on at least two axes, so the swap is nonzero");
Frame3D::new(origin, z, x_dir).expect("x_dir constructed perpendicular to z by design")
}
#[cfg(test)]
mod tests {
use super::*;
use crate::Vector3D;
use std::f64::consts::{FRAC_PI_2, FRAC_PI_4};
#[test]
fn test_from_frame_ok() {
let c = Cone::from_frame(Frame3D::WORLD, 0.4, 2.0).unwrap();
assert_eq!(c.semi_angle(), 0.4);
assert_eq!(c.ref_radius(), 2.0);
}
#[test]
fn test_from_frame_negative_angle_ok() {
let c = Cone::from_frame(Frame3D::WORLD, -0.4, 2.0).unwrap();
assert_eq!(c.semi_angle(), -0.4);
}
#[test]
fn test_from_frame_negative_radius_errors() {
assert_eq!(
Cone::from_frame(Frame3D::WORLD, 0.4, -1.0),
Err(ConeConstructionError::NegativeRadius)
);
}
#[test]
fn test_from_frame_zero_angle_errors() {
assert_eq!(
Cone::from_frame(Frame3D::WORLD, 0.0, 1.0),
Err(ConeConstructionError::BadAngle)
);
}
#[test]
fn test_from_frame_right_angle_errors() {
assert_eq!(
Cone::from_frame(Frame3D::WORLD, FRAC_PI_2, 1.0),
Err(ConeConstructionError::BadAngle)
);
assert_eq!(
Cone::from_frame(Frame3D::WORLD, -FRAC_PI_2, 1.0),
Err(ConeConstructionError::BadAngle)
);
}
#[test]
fn test_from_frame_valid_boundary_ok() {
assert!(Cone::from_frame(Frame3D::WORLD, FRAC_PI_4, 1.0).is_ok());
}
#[test]
fn test_new_ok() {
let c = Cone::new(Point3D::ORIGIN, Vector3D::Z, 0.4, 2.0).unwrap();
assert_eq!(c.frame().z_direction(), Vector3D::Z);
}
#[test]
fn test_new_null_normal_errors() {
assert_eq!(
Cone::new(Point3D::ORIGIN, Vector3D::ZERO, 0.4, 2.0),
Err(ConeConstructionError::NullNormal)
);
}
#[test]
fn test_from_two_points_and_radii_null_axis_errors() {
let p = Point3D::new(1.0, 2.0, 3.0);
assert_eq!(
Cone::from_two_points_and_radii(p, p, 1.0, 2.0),
Err(ConeConstructionError::NullAxis)
);
}
#[test]
fn test_from_two_points_and_radii_negative_radius_errors() {
assert_eq!(
Cone::from_two_points_and_radii(
Point3D::ORIGIN,
Point3D::new(0.0, 0.0, 1.0),
-1.0,
2.0
),
Err(ConeConstructionError::NegativeRadius)
);
assert_eq!(
Cone::from_two_points_and_radii(
Point3D::ORIGIN,
Point3D::new(0.0, 0.0, 1.0),
1.0,
-2.0
),
Err(ConeConstructionError::NegativeRadius)
);
}
#[test]
fn test_from_two_points_and_radii_null_angle_errors_when_radii_equal() {
assert_eq!(
Cone::from_two_points_and_radii(Point3D::ORIGIN, Point3D::new(0.0, 0.0, 1.0), 2.0, 2.0),
Err(ConeConstructionError::NullAngle)
);
}
#[test]
fn test_from_two_points_and_radii_matches_golden_case_shrinking() {
let c =
Cone::from_two_points_and_radii(Point3D::ORIGIN, Point3D::new(0.0, 0.0, 4.0), 2.0, 1.0)
.unwrap();
assert!((c.semi_angle() - (-0.24497866312686414)).abs() < 1e-9);
assert_eq!(c.ref_radius(), 2.0);
let frame = c.frame();
assert_eq!(frame.origin(), Point3D::ORIGIN);
assert_vector3_close(frame.x_direction(), Vector3D::new(0.0, -1.0, 0.0));
assert_vector3_close(frame.y_direction(), Vector3D::new(1.0, 0.0, 0.0));
assert_vector3_close(frame.z_direction(), Vector3D::Z);
}
#[test]
fn test_from_two_points_and_radii_matches_golden_case_growing() {
let c = Cone::from_two_points_and_radii(
Point3D::new(1.0, 1.0, 1.0),
Point3D::new(2.0, 3.0, 1.0),
1.0,
2.5,
)
.unwrap();
assert!((c.semi_angle() - 0.5908727501454191).abs() < 1e-9);
assert_eq!(c.ref_radius(), 1.0);
let frame = c.frame();
assert_eq!(frame.origin(), Point3D::new(1.0, 1.0, 1.0));
assert_vector3_close(
frame.x_direction(),
Vector3D::new(-0.894427190999916, 0.447213595499958, 0.0),
);
assert_vector3_close(frame.y_direction(), Vector3D::new(0.0, 0.0, 1.0));
assert_vector3_close(
frame.z_direction(),
Vector3D::new(0.4472135954999579, 0.8944271909999159, 0.0),
);
}
fn assert_vector3_close(actual: Vector3D, expected: Vector3D) {
assert!(
(actual.x - expected.x).abs() < 1e-9,
"x: {actual:?} vs {expected:?}"
);
assert!(
(actual.y - expected.y).abs() < 1e-9,
"y: {actual:?} vs {expected:?}"
);
assert!(
(actual.z - expected.z).abs() < 1e-9,
"z: {actual:?} vs {expected:?}"
);
}
#[test]
fn test_eval_point() {
let c = Cone::new(Point3D::ORIGIN, Vector3D::Z, 0.4, 2.0).unwrap();
assert_eq!(c.eval_point(0.0, 0.0), Point3D::new(2.0, 0.0, 0.0));
}
#[test]
fn test_eval_points_matches_loop() {
let c = Cone::new(Point3D::ORIGIN, Vector3D::Z, 0.4, 2.0).unwrap();
let uvs = [(0.0, 0.0), (0.5, 1.0)];
let expected: Vec<Point3D> = uvs.iter().map(|&(u, v)| c.eval_point(u, v)).collect();
assert_eq!(c.eval_points(&uvs), expected);
}
#[test]
fn test_apex() {
let c = Cone::new(Point3D::ORIGIN, Vector3D::Z, FRAC_PI_4, 1.0).unwrap();
let apex = c.apex();
assert!(apex.x.abs() < 1e-9);
assert!(apex.y.abs() < 1e-9);
assert!((apex.z + 1.0).abs() < 1e-9);
}
#[test]
#[should_panic(expected = "du + dv must be >= 1")]
fn test_eval_derivative_zero_order_panics() {
let c = Cone::new(Point3D::ORIGIN, Vector3D::Z, 0.4, 2.0).unwrap();
c.eval_derivative(0.0, 0.0, 0, 0);
}
#[test]
#[should_panic(expected = "du + dv must be <= 2")]
fn test_eval_derivative_order_too_high_panics() {
let c = Cone::new(Point3D::ORIGIN, Vector3D::Z, 0.4, 2.0).unwrap();
c.eval_derivative(0.0, 0.0, 2, 1);
}
#[test]
fn test_parameters_of_round_trip() {
let c = Cone::new(Point3D::ORIGIN, Vector3D::Z, 0.4, 2.0).unwrap();
let (u, v) = c.parameters_of(Point3D::new(2.0, 0.0, 0.0));
assert_eq!((u, v), (0.0, 0.0));
}
#[test]
fn test_error_display() {
assert_eq!(
ConeConstructionError::NegativeRadius.to_string(),
"radius is negative"
);
assert_eq!(
ConeConstructionError::BadAngle.to_string(),
"semi-angle is out of range"
);
assert_eq!(
ConeConstructionError::NullAxis.to_string(),
"the two points are confused"
);
assert_eq!(
ConeConstructionError::NullAngle.to_string(),
"semi-angle is null or a right angle"
);
assert_eq!(
ConeConstructionError::NullNormal.to_string(),
"axis direction has zero length"
);
}
#[test]
fn test_error_is_std_error() {
fn takes_error(_e: &dyn std::error::Error) {}
takes_error(&ConeConstructionError::NegativeRadius);
}
#[test]
fn test_cone_contains() {
let cone = Cone::from_frame(Frame3D::WORLD, 0.4, 2.0).unwrap();
let tol = Tolerance::DEFAULT;
assert!(cone.contains(cone.eval_point(0.0, 0.0), tol));
assert!(cone.contains(cone.eval_point(1.0, 2.0), tol));
assert!(!cone.contains(Point3D::ORIGIN, tol));
assert!(!cone.contains(Point3D::new(0.0, 0.0, 10.0), tol));
}
#[test]
fn test_cone_project_point() {
let cone = Cone::from_frame(Frame3D::WORLD, 0.4, 2.0).unwrap();
let tol = Tolerance::DEFAULT;
let p = cone.eval_point(1.0, 2.0);
let proj = cone.project_point(p, tol);
assert!((proj.u - 1.0).abs() < 1e-9);
assert!((proj.v - 2.0).abs() < 1e-9);
assert_eq!(proj.distance, 0.0);
}
}