Divide intersection edges from their neighbours' parameters (#1038) (#1412)

Thanks @jijinbei.
2026-10-16 00:27:55 -04:00
1 parent 42d4ebeb16
commit 3333f8a62b
@@ -0,0 +1,161 @@
+use super::*;
+use std::f64::consts::PI;
+
+impl<P> UnitCircle<P> {
+    /// constructor
+    #[inline]
+    pub const fn new() -> Self { Self(std::marker::PhantomData) }
+}
+
+impl ParametricCurve for UnitCircle<Point2> {
+    type Point = Point2;
+    type Vector = Vector2;
+    #[inline]
+    fn subs(&self, t: f64) -> Self::Point { Point2::new(f64::cos(t), f64::sin(t)) }
+    #[inline]
+    fn der(&self, t: f64) -> Self::Vector { Vector2::new(-f64::sin(t), f64::cos(t)) }
+    #[inline]
+    fn der2(&self, t: f64) -> Self::Vector { Vector2::new(-f64::cos(t), -f64::sin(t)) }
+    #[inline]
+    fn parameter_range(&self) -> ParameterRange {
+        (Bound::Included(0.0), Bound::Excluded(2.0 * PI))
+    }
+}
+
+impl BoundedCurve for UnitCircle<Point2> {}
+
+impl ParametricCurve for UnitCircle<Point3> {
+    type Point = Point3;
+    type Vector = Vector3;
+    #[inline]
+    fn subs(&self, t: f64) -> Self::Point { Point3::new(f64::cos(t), f64::sin(t), 0.0) }
+    #[inline]
+    fn der(&self, t: f64) -> Self::Vector { Vector3::new(-f64::sin(t), f64::cos(t), 0.0) }
+    #[inline]
+    fn der2(&self, t: f64) -> Self::Vector { Vector3::new(-f64::cos(t), -f64::sin(t), 0.0) }
+    #[inline]
+    fn period(&self) -> Option<f64> { Some(2.0 * PI) }
+    #[inline]
+    fn parameter_range(&self) -> ParameterRange {
+        (Bound::Included(0.0), Bound::Excluded(2.0 * PI))
+    }
+}
+
+impl BoundedCurve for UnitCircle<Point3> {}
+
+impl<P> ParameterDivision1D for UnitCircle<P>
+where UnitCircle<P>: ParametricCurve<Point = P>
+{
+    type Point = P;
+    fn parameter_division(&self, range: (f64, f64), tol: f64) -> (Vec<f64>, Vec<P>) {
+        nonpositive_tolerance!(tol);
+        let tol = f64::min(tol, 0.8);
+        let delta = 2.0 * f64::acos(1.0 - tol);
+        let n = 1 + ((range.1 - range.0) / delta) as usize;
+        let params = (0..=n)
+            .map(|i| {
+                let t = i as f64 / n as f64;
+                range.0 * (1.0 - t) + range.1 * t
+            })
+            .collect::<Vec<_>>();
+        let pts = params.iter().map(|t| self.subs(*t)).collect();
+        (params, pts)
+    }
+}
+
+impl SearchNearestParameter<D1> for UnitCircle<Point2> {
+    type Point = Point2;
+    fn search_nearest_parameter<H: Into<SPHint1D>>(
+        &self,
+        pt: Point2,
+        _: H,
+        _: usize,
+    ) -> Option<f64> {
+        let v = pt.to_vec();
+        if v.magnitude().so_small() {
+            return None;
+        }
+        let v = v.normalize();
+        let theta = f64::acos(f64::clamp(v.x, -1.0, 1.0));
+        let theta = match v.y > 0.0 {
+            true => theta,
+            false => 2.0 * PI - theta,
+        };
+        Some(theta)
+    }
+}
+
+#[test]
+fn search_nearest_parameter() {
+    const N: usize = 100;
+    let circle = UnitCircle::<Point2>::new();
+    for i in 0..N {
+        let t = 2.0 * PI * i as f64 / N as f64;
+        let a = 5.0 * rand::random::<f64>() + 0.1;
+        let p = a * circle.subs(t);
+        let s = circle.search_nearest_parameter(p, None, 1).unwrap();
+        let q = a * circle.subs(s);
+        assert_near!(p, q);
+    }
+}
+
+impl SearchParameter<D1> for UnitCircle<Point2> {
+    type Point = Point2;
+    fn search_parameter<H: Into<SPHint1D>>(&self, pt: Point2, _: H, _: usize) -> Option<f64> {
+        let v = pt.to_vec();
+        if !v.magnitude().near(&1.0) {
+            return None;
+        }
+        let v = v.normalize();
+        let theta = f64::acos(f64::clamp(v.x, -1.0, 1.0));
+        let theta = match v.y > 0.0 {
+            true => theta,
+            false => 2.0 * PI - theta,
+        };
+        Some(theta)
+    }
+}
+
+#[test]
+fn search_parameter() {
+    const N: usize = 100;
+    let circle = UnitCircle::<Point2>::new();
+    for i in 1..N {
+        let t = 2.0 * PI * i as f64 / N as f64;
+        let p = circle.subs(t);
+        let s = circle.search_parameter(p, None, 1).unwrap();
+        assert_near!(s, t);
+    }
+}
+
+impl SearchNearestParameter<D1> for UnitCircle<Point3> {
+    type Point = Point3;
+    fn search_nearest_parameter<H: Into<SPHint1D>>(
+        &self,
+        pt: Point3,
+        _: H,
+        _: usize,
+    ) -> Option<f64> {
+        UnitCircle::<Point2>::new().search_nearest_parameter(Point2::new(pt.x, pt.y), None, 0)
+    }
+}
+
+impl SearchParameter<D1> for UnitCircle<Point3> {
+    type Point = Point3;
+    fn search_parameter<H: Into<SPHint1D>>(&self, pt: Point3, _: H, _: usize) -> Option<f64> {
+        if !f64::abs(pt.z).so_small() {
+            return None;
+        }
+        UnitCircle::<Point2>::new().search_parameter(Point2::new(pt.x, pt.y), None, 0)
+    }
+}
+
+#[test]
+fn parameter_division() {
+    let c = UnitCircle::<Point2>::new();
+    let (_div, pts) = c.parameter_division(c.range_tuple(), 0.05);
+    for a in pts.windows(2) {
+        let p = a[0].midpoint(a[1]);
+        assert!(p.to_vec().magnitude() > 0.95);
+    }
+}
@@ -0,0 +1,125 @@
+use super::*;
+
+impl<P> UnitHyperbola<P> {
+    /// constructor
+    #[inline]
+    pub const fn new() -> UnitHyperbola<P> { UnitHyperbola(std::marker::PhantomData) }
+}
+
+impl ParametricCurve for UnitHyperbola<Point2> {
+    type Point = Point2;
+    type Vector = Vector2;
+    #[inline]
+    fn subs(&self, t: f64) -> Self::Point { Point2::new(f64::cosh(t), f64::sinh(t)) }
+    #[inline]
+    fn der(&self, t: f64) -> Self::Vector { Vector2::new(f64::sinh(t), f64::cosh(t)) }
+    #[inline]
+    fn der2(&self, t: f64) -> Self::Vector { Vector2::new(f64::cosh(t), f64::sinh(t)) }
+}
+
+impl ParametricCurve for UnitHyperbola<Point3> {
+    type Point = Point3;
+    type Vector = Vector3;
+    #[inline]
+    fn subs(&self, t: f64) -> Self::Point { Point3::new(f64::cosh(t), f64::sinh(t), 0.0) }
+    #[inline]
+    fn der(&self, t: f64) -> Self::Vector { Vector3::new(f64::sinh(t), f64::cosh(t), 0.0) }
+    #[inline]
+    fn der2(&self, t: f64) -> Self::Vector { Vector3::new(f64::cosh(t), f64::sinh(t), 0.0) }
+}
+
+impl<P> ParameterDivision1D for UnitHyperbola<P>
+where
+    UnitHyperbola<P>: ParametricCurve<Point = P>,
+    P: EuclideanSpace<Scalar = f64> + MetricSpace<Metric = f64> + HashGen<f64>,
+{
+    type Point = P;
+    fn parameter_division(&self, range: (f64, f64), tol: f64) -> (Vec<f64>, Vec<P>) {
+        algo::curve::parameter_division(self, range, tol)
+    }
+}
+
+impl SearchNearestParameter<D1> for UnitHyperbola<Point2> {
+    type Point = Point2;
+    fn search_nearest_parameter<H: Into<SPHint1D>>(
+        &self,
+        p: Point2,
+        _: H,
+        _: usize,
+    ) -> Option<f64> {
+        let a = -p.y;
+        let b = (p.y * p.y - p.x * p.x) / 4.0 + 1.0;
+        let c = -p.y;
+        let d = p.y * p.y / 4.0;
+        let y = solver::solve_quartic(a, b, c, d)
+            .into_iter()
+            .filter_map(|z| match z.im.so_small() {
+                true => Some(z.re),
+                false => None,
+            })
+            .min_by(|s, t| {
+                p.distance2(self.subs(*s))
+                    .partial_cmp(&p.distance2(self.subs(*t)))
+                    .unwrap()
+            })?;
+        Some(f64::asinh(y))
+    }
+}
+
+impl SearchNearestParameter<D1> for UnitHyperbola<Point3> {
+    type Point = Point3;
+    fn search_nearest_parameter<H: Into<SPHint1D>>(
+        &self,
+        p: Point3,
+        _: H,
+        _trials: usize,
+    ) -> Option<f64> {
+        UnitHyperbola::<Point2>::new().search_nearest_parameter(
+            Point2::new(p.x, p.y),
+            None,
+            _trials,
+        )
+    }
+}
+
+impl SearchParameter<D1> for UnitHyperbola<Point2> {
+    type Point = Point2;
+    fn search_parameter<H: Into<SPHint1D>>(&self, p: Point2, _: H, _: usize) -> Option<f64> {
+        let t = f64::asinh(p.y);
+        match p.near(&self.subs(t)) {
+            true => Some(t),
+            false => None,
+        }
+    }
+}
+
+impl SearchParameter<D1> for UnitHyperbola<Point3> {
+    type Point = Point3;
+    fn search_parameter<H: Into<SPHint1D>>(&self, p: Point3, _: H, _: usize) -> Option<f64> {
+        let t = f64::asinh(p.y);
+        match p.near(&self.subs(t)) {
+            true => Some(t),
+            false => None,
+        }
+    }
+}
+
+#[test]
+fn snp_test() {
+    let curve = UnitHyperbola::<Point2>::new();
+    let p = curve.subs(2.0);
+    let q = p + Vector2::new(-p.x, p.y);
+    let t = curve.search_nearest_parameter(q, None, 0).unwrap();
+    assert_near!(t, 2.0);
+}
+
+#[test]
+fn sp_test() {
+    let curve = UnitHyperbola::<Point2>::new();
+    let t = 100.0 * rand::random::<f64>() - 50.0;
+    let p = curve.subs(t);
+    assert_near!(curve.search_parameter(p, None, 0).unwrap(), t);
+
+    let q = Point2::new(-1.0, 0.0);
+    assert!(curve.search_parameter(q, None, 0).is_none());
+}
@@ -0,0 +1,217 @@
+use super::*;
+use truck_base::cgmath64::control_point::ControlPoint;
+
+impl<P: Copy> Line<P> {
+    /// initialize line from vector
+    #[inline]
+    pub fn from_origin_direction<V>(origin: P, direction: V) -> Self
+    where P: std::ops::Add<V, Output = P> {
+        Self(origin, origin + direction)
+    }
+    /// to a bspline curve
+    #[inline]
+    pub fn to_bspline(self) -> BSplineCurve<P> {
+        BSplineCurve::new(KnotVec::bezier_knot(1), vec![self.0, self.1])
+    }
+}
+
+impl<P> Line<P>
+where
+    P: EuclideanSpace<Scalar = f64>,
+    P::Diff: InnerSpace<Scalar = f64>,
+{
+    /// Returns the projected point to the line.
+    /// # Examples
+    /// ```
+    /// use truck_geometry::prelude::*;
+    /// let line = Line(Point2::new(0.0, 0.0), Point2::new(1.0, 2.0));
+    /// let pt = Point2::new(0.0, 1.0);
+    /// assert_near!(line.projection(pt), Point2::new(0.4, 0.8));
+    /// ```
+    pub fn projection(self, point: P) -> P {
+        let (u, v) = (point - self.0, self.1 - self.0);
+        self.0 + v * u.dot(v) / v.dot(v)
+    }
+    /// Returns the distance between the line and the point `point`.
+    /// # Examples
+    /// ```
+    /// use truck_geometry::prelude::*;
+    /// let line = Line(Point2::new(0.0, 0.0), Point2::new(3.0, 4.0));
+    ///
+    /// // The foot of the perpendicular line is on a line segment
+    /// let pt = Point2::new(3.0, 0.0);
+    /// assert_near!(line.distance_to_point(pt), 2.4);
+    ///
+    /// // The foot of the perpendicular line is not on a line segment
+    /// let pt = Point2::new(0.0, -4.0);
+    /// assert_near!(line.distance_to_point(pt), 2.4);
+    /// ```
+    pub fn distance_to_point(self, point: P) -> f64 {
+        let (u, v) = (point - self.0, self.1 - self.0);
+        (u - v * u.dot(v) / v.dot(v)).magnitude()
+    }
+    /// Returns the distance between the sengment and the point `point`.
+    /// # Examples
+    /// ```
+    /// use truck_geometry::prelude::*;
+    /// let line = Line(Point2::new(0.0, 0.0), Point2::new(3.0, 4.0));
+    ///
+    /// // The foot of the perpendicular line is on a line segment
+    /// let pt = Point2::new(3.0, 0.0);
+    /// assert_near!(line.distance_to_point_as_segment(pt), 2.4);
+    ///
+    /// // The foot of the perpendicular line is not on a line segment
+    /// let pt = Point2::new(0.0, -4.0);
+    /// assert_near!(line.distance_to_point_as_segment(pt), 4.0);
+    /// ```
+    pub fn distance_to_point_as_segment(self, point: P) -> f64 {
+        let (u, v) = (point - self.0, self.1 - self.0);
+        let t = f64::clamp(u.dot(v) / v.dot(v), 0.0, 1.0);
+        (u - v * t).magnitude()
+    }
+}
+
+impl Line<Point2> {
+    /// Returns the intersection of two lines and its parameters.
+    /// # Examples
+    /// ```
+    /// use truck_geometry::prelude::*;
+    /// let line0 = Line(Point2::new(0.0, 0.0), Point2::new(9.0, 3.0));
+    /// let line1 = Line(Point2::new(0.0, 6.0), Point2::new(9.0, 0.0));
+    /// let (s, t, p) = line0.intersection(line1).unwrap();
+    /// assert_near!(line0.subs(s), Point2::new(6.0, 2.0));
+    /// assert_near!(line1.subs(t), Point2::new(6.0, 2.0));
+    /// assert_near!(p, Point2::new(6.0, 2.0));
+    /// ```
+    /// # Failures
+    /// Returns `None` if two lines are parallel.
+    /// ```
+    /// use truck_geometry::prelude::*;
+    /// let line0 = Line(Point2::new(0.0, 0.0), Point2::new(9.0, 3.0));
+    /// let line1 = Line(Point2::new(-1.0, 0.0), Point2::new(8.0, 3.0));
+    /// assert!(line0.intersection(line1).is_none());
+    /// ```
+    pub fn intersection(self, other: Line<Point2>) -> Option<(f64, f64, Point2)> {
+        let mat = Matrix2::from_cols(self.1 - self.0, other.0 - other.1);
+        let v = other.0 - self.0;
+        let params = mat.invert().map(|inv| inv * v)?;
+        Some((params.x, params.y, self.subs(params.x)))
+    }
+}
+
+impl<P: ControlPoint<f64>> ParametricCurve for Line<P> {
+    type Point = P;
+    type Vector = P::Diff;
+    #[inline]
+    fn subs(&self, t: f64) -> Self::Point { self.0 + (self.1 - self.0) * t }
+    #[inline]
+    fn der(&self, _: f64) -> Self::Vector { self.1 - self.0 }
+    #[inline]
+    fn der2(&self, _: f64) -> Self::Vector { Self::Vector::zero() }
+    /// Return `0.0..=1.0` i.e. we regard it as a segment
+    #[inline]
+    fn parameter_range(&self) -> ParameterRange { (Bound::Included(0.0), Bound::Included(1.0)) }
+}
+
+impl<P: ControlPoint<f64>> BoundedCurve for Line<P> {}
+
+impl<P: ControlPoint<f64>> Cut for Line<P> {
+    #[inline]
+    fn cut(&mut self, t: f64) -> Self {
+        let r = self.subs(t);
+        let res = Self(r, self.1);
+        self.1 = r;
+        res
+    }
+}
+
+impl<P: ControlPoint<f64>> ParameterDivision1D for Line<P> {
+    type Point = P;
+    #[inline]
+    fn parameter_division(&self, range: (f64, f64), _: f64) -> (Vec<f64>, Vec<P>) {
+        (
+            vec![range.0, range.1],
+            vec![self.subs(range.0), self.subs(range.1)],
+        )
+    }
+}
+
+impl<P: Copy> Invertible for Line<P> {
+    #[inline]
+    fn invert(&mut self) { *self = Self(self.1, self.0); }
+    #[inline]
+    fn inverse(&self) -> Self { Self(self.1, self.0) }
+}
+
+impl<P> SearchNearestParameter<D1> for Line<P>
+where
+    P: ControlPoint<f64>,
+    P::Diff: InnerSpace<Scalar = f64>,
+{
+    type Point = P;
+    #[inline]
+    fn search_nearest_parameter<H: Into<SPHint1D>>(&self, pt: P, _: H, _: usize) -> Option<f64> {
+        let b = self.1 - self.0;
+        Some((pt - self.0).dot(b) / b.dot(b))
+    }
+}
+
+impl<P> SearchParameter<D1> for Line<P>
+where
+    P: ControlPoint<f64> + Tolerance,
+    P::Diff: InnerSpace<Scalar = f64>,
+{
+    type Point = P;
+    #[inline]
+    fn search_parameter<H: Into<SPHint1D>>(&self, pt: P, _: H, _: usize) -> Option<f64> {
+        let b = self.1 - self.0;
+        let t = (pt - self.0).dot(b) / b.dot(b);
+        match self.subs(t).near(&pt) {
+            true => Some(t),
+            false => None,
+        }
+    }
+}
+
+impl<P: EuclideanSpace, M: Transform<P>> Transformed<M> for Line<P> {
+    #[inline]
+    fn transform_by(&mut self, trans: M) {
+        self.0 = trans.transform_point(self.0);
+        self.1 = trans.transform_point(self.1);
+    }
+    #[inline]
+    fn transformed(&self, trans: M) -> Self {
+        Line(trans.transform_point(self.0), trans.transform_point(self.1))
+    }
+}
+
+#[test]
+fn line() {
+    let line = Line(Point2::new(1.0, 0.0), Point2::new(0.0, 1.0));
+
+    // subs
+    assert_near!(line.subs(0.4), Point2::new(0.6, 0.4));
+
+    // inverse
+    let line_inverse = line.inverse();
+    assert_eq!(line.0, line_inverse.1);
+    assert_eq!(line.1, line_inverse.0);
+
+    // cut
+    let mut line0 = line;
+    let line1 = line0.cut(0.4);
+    assert_eq!(line.0, line0.0);
+    assert_near!(line0.1, line.subs(0.4));
+    assert_eq!(line0.1, line1.0);
+    assert_eq!(line1.1, line.1);
+
+    // SNP
+    assert_near!(
+        line.search_nearest_parameter(Point2::new(1.0, 1.0), None, 0)
+            .unwrap(),
+        0.5
+    );
+    assert!(line
+        .search_parameter(Point2::new(1.0, 1.0), None, 0)
+        .is_none());
+}
@@ -0,0 +1,115 @@
+use crate::{prelude::*, *};
+
+/// line
+/// # Example
+/// ```
+/// use truck_geometry::prelude::*;
+/// let line = Line(Point2::new(0.0, 0.0), Point2::new(1.0, 1.0));
+/// assert_near!(line.subs(0.5), Point2::new(0.5, 0.5));
+/// ```
+#[derive(Clone, Copy, Debug, PartialEq, Eq, Serialize, Deserialize)]
+pub struct Line<P>(pub P, pub P);
+
+/// unit circle
+#[derive(Clone, Copy, Debug, Default, Serialize, Deserialize)]
+pub struct UnitCircle<P>(std::marker::PhantomData<P>);
+
+/// unit hyperbola
+#[derive(Clone, Copy, Debug, Default, Serialize, Deserialize)]
+pub struct UnitHyperbola<P>(std::marker::PhantomData<P>);
+
+/// parabola whose apex is the origin.
+#[derive(Clone, Copy, Debug, Default, Serialize, Deserialize)]
+pub struct UnitParabola<P>(std::marker::PhantomData<P>);
+
+/// plane
+/// # Example
+/// ```
+/// use truck_geometry::prelude::*;
+///
+/// // arbitrary three points
+/// let pt0 = Point3::new(0.0, 1.0, 2.0);
+/// let pt1 = Point3::new(1.0, 1.0, 3.0);
+/// let pt2 = Point3::new(0.0, 2.0, 3.0);
+///
+/// // Creates a plane
+/// let plane: Plane = Plane::new(pt0, pt1, pt2);
+/// // The origin of the plane is pt0.
+/// assert_near!(plane.origin(), pt0);
+/// // The u-axis of the plane is the vector from pt0 to pt1.
+/// assert_near!(plane.u_axis(), pt1 - pt0);
+/// // The v-axis of the plane is the vector from pt0 to pt2.
+/// assert_near!(plane.v_axis(), pt2 - pt0);
+/// // The normal is the normalized u-axis × v-axis
+/// assert_near!(plane.normal(), (pt1 - pt0).cross(pt2 - pt0).normalize());
+/// ```
+#[derive(Clone, Copy, Debug, PartialEq, Serialize, Deserialize)]
+pub struct Plane {
+    o: Point3,
+    p: Point3,
+    q: Point3,
+}
+
+/// sphere
+/// # Examples
+/// ```
+/// use truck_geometry::prelude::*;
+/// use std::f64::consts::PI;
+///
+/// let center = Point3::new(1.0, 2.0, 3.0);
+/// let radius = 4.56;
+///
+/// let sphere = Sphere::new(center, radius);
+/// const N: usize = 100;
+/// for i in 0..=N {
+///     for j in 0..=N {
+///         // the parameter u is latitude
+///         let u = PI * i as f64 / N as f64;
+///         // the parameter v is longitude
+///         let v = 2.0 * PI * j as f64 / N as f64;
+///
+///         // simple relation between a point and its normal.
+///         let pt = sphere.subs(u, v);
+///         let n = sphere.normal(u, v);
+///         assert_near!(pt - center, n * radius);
+///
+///         // the proof of u is latitude and v is longitude
+///         assert!((PI / 2.0 - u) * (pt.z - center.z) >= 0.0);
+///         assert!((PI - v) * (pt.y - center.y) >= 0.0);
+///     }
+/// }
+/// ```
+#[derive(Clone, Copy, Debug, PartialEq, Serialize, Deserialize)]
+pub struct Sphere {
+    center: Point3,
+    radius: f64,
+}
+
+/// torus
+#[derive(Clone, Copy, Debug, PartialEq, Serialize, Deserialize)]
+pub struct Torus {
+    center: Point3,
+    large_radius: f64,
+    small_radius: f64,
+}
+
+mod circle;
+mod hyperbola;
+mod line;
+mod parabola;
+mod plane;
+mod sphere;
+mod torus;
+
+macro_rules! always_true {
+    ($ty: tt) => {
+        impl<P> PartialEq for $ty<P> {
+            fn eq(&self, _: &Self) -> bool { true }
+        }
+        impl<P> Eq for $ty<P> {}
+    };
+}
+
+always_true!(UnitCircle);
+always_true!(UnitParabola);
+always_true!(UnitHyperbola);
@@ -0,0 +1,137 @@
+use super::*;
+
+impl<P> UnitParabola<P> {
+    /// constructor
+    #[inline]
+    pub const fn new() -> Self { Self(std::marker::PhantomData) }
+}
+
+impl ParametricCurve for UnitParabola<Point2> {
+    type Point = Point2;
+    type Vector = Vector2;
+    #[inline]
+    fn subs(&self, t: f64) -> Self::Point { Point2::new(t * t, 2.0 * t) }
+    #[inline]
+    fn der(&self, t: f64) -> Self::Vector { Vector2::new(2.0 * t, 2.0) }
+    #[inline]
+    fn der2(&self, _: f64) -> Self::Vector { Vector2::new(2.0, 0.0) }
+}
+
+impl ParametricCurve for UnitParabola<Point3> {
+    type Point = Point3;
+    type Vector = Vector3;
+    #[inline]
+    fn subs(&self, t: f64) -> Self::Point { Point3::new(t * t, 2.0 * t, 0.0) }
+    #[inline]
+    fn der(&self, t: f64) -> Self::Vector { Vector3::new(2.0 * t, 2.0, 0.0) }
+    #[inline]
+    fn der2(&self, _: f64) -> Self::Vector { Vector3::new(2.0, 0.0, 0.0) }
+}
+
+impl<P> ParameterDivision1D for UnitParabola<P>
+where
+    UnitParabola<P>: ParametricCurve<Point = P>,
+    P: EuclideanSpace<Scalar = f64> + MetricSpace<Metric = f64> + HashGen<f64>,
+{
+    type Point = P;
+    fn parameter_division(&self, range: (f64, f64), tol: f64) -> (Vec<f64>, Vec<P>) {
+        algo::curve::parameter_division(self, range, tol)
+    }
+}
+
+impl SearchNearestParameter<D1> for UnitParabola<Point2> {
+    type Point = Point2;
+    #[inline]
+    fn search_nearest_parameter<H: Into<SPHint1D>>(
+        &self,
+        pt: Point2,
+        _: H,
+        _: usize,
+    ) -> Option<f64> {
+        let p = 2.0 - pt.x;
+        let q = -pt.y;
+        solver::pre_solve_cubic(p, q)
+            .into_iter()
+            .filter_map(|x| match x.im.so_small() {
+                true => Some(x.re),
+                false => None,
+            })
+            .min_by(|s, t| {
+                pt.distance2(self.subs(*s))
+                    .partial_cmp(&pt.distance2(self.subs(*t)))
+                    .unwrap()
+            })
+    }
+}
+
+impl SearchNearestParameter<D1> for UnitParabola<Point3> {
+    type Point = Point3;
+    #[inline]
+    fn search_nearest_parameter<H: Into<SPHint1D>>(
+        &self,
+        pt: Point3,
+        _hint: H,
+        _trials: usize,
+    ) -> Option<f64> {
+        UnitParabola::<Point2>::new().search_nearest_parameter(
+            Point2::new(pt.x, pt.y),
+            _hint,
+            _trials,
+        )
+    }
+}
+
+impl SearchParameter<D1> for UnitParabola<Point2> {
+    type Point = Point2;
+    #[inline]
+    fn search_parameter<H: Into<SPHint1D>>(&self, pt: Point2, _: H, _: usize) -> Option<f64> {
+        let t = pt.y / 2.0;
+        let pt0 = self.subs(t);
+        match pt.near(&pt0) {
+            true => Some(t),
+            false => None,
+        }
+    }
+}
+
+impl SearchParameter<D1> for UnitParabola<Point3> {
+    type Point = Point3;
+    #[inline]
+    fn search_parameter<H: Into<SPHint1D>>(
+        &self,
+        pt: Point3,
+        _hint: H,
+        _trials: usize,
+    ) -> Option<f64> {
+        match pt.z.so_small() {
+            true => UnitParabola::<Point2>::new().search_parameter(
+                Point2::new(pt.x, pt.y),
+                _hint,
+                _trials,
+            ),
+            false => None,
+        }
+    }
+}
+
+#[test]
+fn snp_test() {
+    let curve = UnitParabola::<Point2>::new();
+
+    let p = Point2::new(-3.0, 0.0);
+    assert_near!(curve.search_nearest_parameter(p, None, 0).unwrap(), 0.0);
+    let p = Point2::new(-3.0, 6.0);
+    assert_near!(curve.search_nearest_parameter(p, None, 0).unwrap(), 1.0);
+    let p = Point2::new(1.5, 1.5);
+    assert_near!(curve.search_nearest_parameter(p, None, 0).unwrap(), 1.0);
+}
+
+#[test]
+fn sp_test() {
+    let curve = UnitParabola::<Point2>::new();
+
+    let p = Point2::new(4.0, -4.0);
+    assert_near!(curve.search_parameter(p, None, 0).unwrap(), -2.0);
+    let p = Point2::new(-3.0, 6.0);
+    assert!(curve.search_parameter(p, None, 0).is_none());
+}
@@ -0,0 +1,256 @@
+use super::*;
+
+impl Plane {
+    /// Creates a new plane from three points.
+    #[inline(always)]
+    pub const fn new(origin: Point3, one: Point3, another: Point3) -> Plane {
+        Plane {
+            o: origin,
+            p: one,
+            q: another,
+        }
+    }
+    /// Returns the origin
+    #[inline(always)]
+    pub const fn origin(&self) -> Point3 { self.o }
+    /// Returns the u-axis
+    #[inline(always)]
+    pub fn u_axis(&self) -> Vector3 { self.p - self.o }
+    /// Returns the v-axis
+    #[inline(always)]
+    pub fn v_axis(&self) -> Vector3 { self.q - self.o }
+    /// Returns the normal
+    /// # Examples
+    /// ```
+    /// use truck_geometry::prelude::*;
+    /// let plane = Plane::new(
+    ///     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_near!(plane.normal(), Vector3::unit_z());
+    /// ```
+    #[inline(always)]
+    pub fn normal(&self) -> Vector3 { self.u_axis().cross(self.v_axis()).normalize() }
+    /// Gets the parameter of `pt` in plane's matrix.
+    /// # Examples
+    /// ```
+    /// use truck_geometry::prelude::*;
+    /// let plane = Plane::new(
+    ///     Point3::new(1.0, 2.0, 3.0),
+    ///     Point3::new(2.0, 1.0, 3.0),
+    ///     Point3::new(3.0, 4.0, -1.0),
+    /// );
+    ///
+    /// let pt = Point3::new(2.1, -6.5, 4.7);
+    /// let prm = plane.get_parameter(pt);
+    /// let rev = plane.origin()
+    ///     + prm[0] * plane.u_axis()
+    ///     + prm[1] * plane.v_axis()
+    ///     + prm[2] * plane.normal();
+    /// assert_near!(pt, rev);
+    /// ```
+    #[inline(always)]
+    pub fn get_parameter(&self, pt: Point3) -> Vector3 {
+        let a = self.u_axis();
+        let b = self.v_axis();
+        let c = self.normal();
+        let mat = Matrix3::from_cols(a, b, c).invert().unwrap();
+        mat * (pt - self.o)
+    }
+    /// into B-spline surface
+    /// # Examples
+    /// ```
+    /// use truck_geometry::prelude::*;
+    /// let pt0 = Point3::new(0.0, 1.0, 2.0);
+    /// let pt1 = Point3::new(1.0, 1.0, 3.0);
+    /// let pt2 = Point3::new(0.0, 2.0, 3.0);
+    /// let plane: Plane = Plane::new(pt0, pt1, pt2);
+    /// let surface: BSplineSurface<Point3> = plane.into_bspline();
+    /// assert_eq!(surface.range_tuple(), ((0.0, 1.0), (0.0, 1.0)));
+    ///
+    /// const N: usize = 100;
+    /// for i in 0..=N {
+    ///     for j in 0..=N {
+    ///         let u = i as f64 / N as f64;
+    ///         let v = j as f64 / N as f64;
+    ///         let res = surface.subs(u, v);
+    ///         let ans = plane.subs(u, v);
+    ///         assert_near!(ans, res);
+    ///     }
+    /// }
+    /// ```
+    #[inline(always)]
+    pub fn into_bspline(&self) -> BSplineSurface<Point3> {
+        let o = self.o;
+        let p = self.p;
+        let q = self.q;
+        BSplineSurface::debug_new(
+            (KnotVec::bezier_knot(1), KnotVec::bezier_knot(1)),
+            vec![vec![o, q], vec![p, p + (q - o)]],
+        )
+    }
+    /// into NURBS surface
+    /// # Examples
+    /// ```
+    /// use truck_geometry::prelude::*;
+    /// let pt0 = Point3::new(0.0, 1.0, 2.0);
+    /// let pt1 = Point3::new(1.0, 1.0, 3.0);
+    /// let pt2 = Point3::new(0.0, 2.0, 3.0);
+    /// let plane: Plane = Plane::new(pt0, pt1, pt2);
+    /// let surface: NurbsSurface<Vector4> = plane.into_nurbs();
+    /// assert_eq!(surface.range_tuple(), ((0.0, 1.0), (0.0, 1.0)));
+    ///
+    /// const N: usize = 100;
+    /// for i in 0..=N {
+    ///     for j in 0..=N {
+    ///         let u = i as f64 / N as f64;
+    ///         let v = j as f64 / N as f64;
+    ///         let res = surface.subs(u, v);
+    ///         let ans = plane.subs(u, v);
+    ///         assert_near!(ans, res);
+    ///     }
+    /// }
+    /// ```
+    #[inline(always)]
+    pub fn into_nurbs(&self) -> NurbsSurface<Vector4> {
+        let o = self.o.to_homogeneous();
+        let p = self.p.to_homogeneous();
+        let q = self.q.to_homogeneous();
+        NurbsSurface::new(BSplineSurface::debug_new(
+            (KnotVec::bezier_knot(1), KnotVec::bezier_knot(1)),
+            vec![vec![o, q], vec![p, p + q - o]],
+        ))
+    }
+}
+
+impl ParametricSurface for Plane {
+    type Point = Point3;
+    type Vector = Vector3;
+    #[inline(always)]
+    fn subs(&self, u: f64, v: f64) -> Point3 {
+        self.o + u * (self.p - self.o) + v * (self.q - self.o)
+    }
+    #[inline(always)]
+    fn uder(&self, _: f64, _: f64) -> Vector3 { self.p - self.o }
+    #[inline(always)]
+    fn vder(&self, _: f64, _: f64) -> Vector3 { self.q - self.o }
+    #[inline(always)]
+    fn uuder(&self, _: f64, _: f64) -> Vector3 { Vector3::zero() }
+    #[inline(always)]
+    fn uvder(&self, _: f64, _: f64) -> Vector3 { Vector3::zero() }
+    #[inline(always)]
+    fn vvder(&self, _: f64, _: f64) -> Vector3 { Vector3::zero() }
+    /// as square
+    #[inline(always)]
+    fn parameter_range(&self) -> (ParameterRange, ParameterRange) {
+        let range = (Bound::Included(0.0), Bound::Included(1.0));
+        (range, range)
+    }
+}
+
+impl ParametricSurface3D for Plane {
+    #[inline(always)]
+    fn normal(&self, _: f64, _: f64) -> Vector3 { self.normal() }
+}
+
+impl BoundedSurface for Plane {}
+
+impl Invertible for Plane {
+    #[inline(always)]
+    fn inverse(&self) -> Self {
+        Plane {
+            o: self.o,
+            p: self.q,
+            q: self.p,
+        }
+    }
+    #[inline(always)]
+    fn invert(&mut self) { *self = self.inverse(); }
+}
+
+impl IncludeCurve<BSplineCurve<Point3>> for Plane {
+    #[inline(always)]
+    fn include(&self, curve: &BSplineCurve<Point3>) -> bool {
+        let origin = self.origin();
+        let normal = self.normal();
+        curve
+            .control_points()
+            .iter()
+            .all(|pt| (pt - origin).dot(normal).so_small())
+    }
+}
+
+impl IncludeCurve<NurbsCurve<Vector4>> for Plane {
+    fn include(&self, curve: &NurbsCurve<Vector4>) -> bool {
+        let origin = self.origin();
+        let normal = self.normal();
+        let (s, e) = (curve.front(), curve.back());
+        if !(s - origin).dot(normal).so_small() || !(e - origin).dot(normal).so_small() {
+            return false;
+        }
+        curve.non_rationalized().control_points().iter().all(|pt| {
+            if pt[3].so_small() {
+                true
+            } else {
+                let pt = Point3::from_homogeneous(*pt);
+                (pt - origin).dot(normal).so_small()
+            }
+        })
+    }
+}
+
+impl ParameterDivision2D for Plane {
+    #[inline(always)]
+    fn parameter_division(&self, range: ((f64, f64), (f64, f64)), _: f64) -> (Vec<f64>, Vec<f64>) {
+        (vec![range.0 .0, range.0 .1], vec![range.1 .0, range.1 .1])
+    }
+}
+
+impl<T: Transform3<Scalar = f64>> Transformed<T> for Plane {
+    #[inline(always)]
+    fn transform_by(&mut self, trans: T) {
+        self.o = trans.transform_point(self.o);
+        self.p = trans.transform_point(self.p);
+        self.q = trans.transform_point(self.q);
+    }
+    #[inline(always)]
+    fn transformed(&self, trans: T) -> Self {
+        Plane {
+            o: trans.transform_point(self.o),
+            p: trans.transform_point(self.p),
+            q: trans.transform_point(self.q),
+        }
+    }
+}
+
+impl SearchParameter<D2> for Plane {
+    type Point = Point3;
+    #[inline(always)]
+    fn search_parameter<H: Into<SPHint2D>>(
+        &self,
+        point: Point3,
+        _: H,
+        _: usize,
+    ) -> Option<(f64, f64)> {
+        let v = self.get_parameter(point);
+        match v[2].so_small() {
+            true => Some((v[0], v[1])),
+            false => None,
+        }
+    }
+}
+
+impl SearchNearestParameter<D2> for Plane {
+    type Point = Point3;
+    #[inline(always)]
+    fn search_nearest_parameter<H: Into<SPHint2D>>(
+        &self,
+        point: Point3,
+        _: H,
+        _: usize,
+    ) -> Option<(f64, f64)> {
+        let v = self.get_parameter(point);
+        Some((v[0], v[1]))
+    }
+}
@@ -0,0 +1,219 @@
+use super::*;
+use std::f64::consts::PI;
+
+impl Sphere {
+    /// Creates a sphere
+    #[inline(always)]
+    pub const fn new(center: Point3, radius: f64) -> Sphere { Sphere { center, radius } }
+    /// Returns the center
+    #[inline(always)]
+    pub const fn center(&self) -> Point3 { self.center }
+    /// Returns the radius
+    #[inline(always)]
+    pub const fn radius(&self) -> f64 { self.radius }
+    /// Returns whether the point `pt` is on sphere
+    #[inline(always)]
+    pub fn include(&self, pt: Point3) -> bool { self.center.distance(pt).near(&self.radius) }
+}
+
+impl ParametricSurface for Sphere {
+    type Point = Point3;
+    type Vector = Vector3;
+    #[inline(always)]
+    fn subs(&self, u: f64, v: f64) -> Point3 { self.center() + self.radius * self.normal(u, v) }
+    #[inline(always)]
+    fn uder(&self, u: f64, v: f64) -> Vector3 {
+        self.radius
+            * Vector3::new(
+                f64::cos(u) * f64::cos(v),
+                f64::cos(u) * f64::sin(v),
+                -f64::sin(u),
+            )
+    }
+    #[inline(always)]
+    fn vder(&self, u: f64, v: f64) -> Vector3 {
+        self.radius * f64::sin(u) * Vector3::new(-f64::sin(v), f64::cos(v), 0.0)
+    }
+    #[inline(always)]
+    fn uuder(&self, u: f64, v: f64) -> Vector3 { -self.radius * self.normal(u, v) }
+    #[inline(always)]
+    fn uvder(&self, u: f64, v: f64) -> Vector3 {
+        self.radius * f64::cos(u) * Vector3::new(-f64::sin(v), f64::cos(v), 0.0)
+    }
+    #[inline(always)]
+    fn vvder(&self, u: f64, v: f64) -> Vector3 {
+        -self.radius * f64::sin(u) * Vector3::new(f64::cos(v), f64::sin(v), 0.0)
+    }
+    #[inline(always)]
+    fn parameter_range(&self) -> (ParameterRange, ParameterRange) {
+        (
+            (Bound::Included(0.0), Bound::Included(PI)),
+            (Bound::Included(0.0), Bound::Excluded(2.0 * PI)),
+        )
+    }
+    #[inline(always)]
+    fn v_period(&self) -> Option<f64> { Some(2.0 * PI) }
+}
+
+impl ParametricSurface3D for Sphere {
+    #[inline(always)]
+    fn normal(&self, u: f64, v: f64) -> Vector3 {
+        Vector3::new(
+            f64::sin(u) * f64::cos(v),
+            f64::sin(u) * f64::sin(v),
+            f64::cos(u),
+        )
+    }
+}
+
+#[test]
+fn sphere_derivation_test() {
+    let center = Point3::new(1.0, 2.0, 3.0);
+    let radius = 4.56;
+    let sphere = Sphere::new(center, radius);
+    const N: usize = 100;
+    for i in 0..N {
+        for j in 0..N {
+            let u = PI * i as f64 / N as f64;
+            let v = 2.0 * PI * j as f64 / N as f64;
+            let normal = sphere.normal(u, v);
+            assert!(normal.dot(sphere.uder(u, v)).so_small());
+            assert!(normal.dot(sphere.vder(u, v)).so_small());
+        }
+    }
+}
+
+impl BoundedSurface for Sphere {}
+
+impl IncludeCurve<BSplineCurve<Point3>> for Sphere {
+    #[inline(always)]
+    fn include(&self, curve: &BSplineCurve<Point3>) -> bool {
+        curve.is_const() && self.include(curve.front())
+    }
+}
+
+impl IncludeCurve<NurbsCurve<Vector4>> for Sphere {
+    fn include(&self, curve: &NurbsCurve<Vector4>) -> bool {
+        let (knots, _) = curve.knot_vec().to_single_multi();
+        let degree = curve.degree() * 2;
+        knots
+            .windows(2)
+            .flat_map(move |window| (1..degree).map(move |i| (window, i)))
+            .all(move |(window, i)| {
+                let t = i as f64 / degree as f64;
+                let t = window[0] * (1.0 - t) + window[1] * t;
+                self.include(curve.subs(t))
+            })
+    }
+}
+
+impl ParameterDivision2D for Sphere {
+    #[inline(always)]
+    fn parameter_division(
+        &self,
+        (urange, vrange): ((f64, f64), (f64, f64)),
+        tol: f64,
+    ) -> (Vec<f64>, Vec<f64>) {
+        nonpositive_tolerance!(tol);
+        assert!(
+            tol < self.radius,
+            "Tolerance is larger than the radius of sphere."
+        );
+        let delta = 2.0 * f64::acos(1.0 - tol / self.radius);
+        let u_div = 1 + ((urange.1 - urange.0) / delta).floor() as usize;
+        let v_div = 1 + ((vrange.1 - vrange.0) / delta).floor() as usize;
+        (
+            (0..=u_div)
+                .map(|i| urange.0 + (urange.1 - urange.0) * i as f64 / u_div as f64)
+                .collect(),
+            (0..=v_div)
+                .map(|j| vrange.0 + (vrange.1 - vrange.0) * j as f64 / v_div as f64)
+                .collect(),
+        )
+    }
+}
+
+impl SearchParameter<D2> for Sphere {
+    type Point = Point3;
+    #[inline(always)]
+    fn search_parameter<H: Into<SPHint2D>>(
+        &self,
+        point: Point3,
+        hint: H,
+        _: usize,
+    ) -> Option<(f64, f64)> {
+        let radius = point - self.center;
+        if (self.radius * self.radius).near(&radius.magnitude2()) {
+            let radius = radius.normalize();
+            let u = f64::acos(radius[2]);
+            let sinu = f64::sqrt(1.0 - radius[2] * radius[2]);
+            let cosv = f64::clamp(radius[0] / sinu, -1.0, 1.0);
+            let v = if sinu.so_small() {
+                match hint.into() {
+                    SPHint2D::Parameter(_, hint) => hint,
+                    _ => 0.0,
+                }
+            } else if radius[1] > 0.0 {
+                f64::acos(cosv)
+            } else {
+                2.0 * PI - f64::acos(cosv)
+            };
+            Some((u, v))
+        } else {
+            None
+        }
+    }
+}
+
+impl SearchNearestParameter<D2> for Sphere {
+    type Point = Point3;
+    #[inline(always)]
+    fn search_nearest_parameter<H: Into<SPHint2D>>(
+        &self,
+        point: Point3,
+        _: H,
+        _: usize,
+    ) -> Option<(f64, f64)> {
+        let radius = (point - self.center).normalize();
+        let u = f64::acos(radius[2]);
+        let sinu = f64::sqrt(1.0 - radius[2] * radius[2]);
+        let cosv = radius[0] / sinu;
+        let v = if radius[1] > 0.0 {
+            f64::acos(cosv)
+        } else {
+            2.0 * PI - f64::acos(cosv)
+        };
+        Some((u, v))
+    }
+}
+
+#[cfg(test)]
+fn exec_search_parameter_test() {
+    let center = Point3::new(
+        100.0 * rand::random::<f64>() - 50.0,
+        100.0 * rand::random::<f64>() - 50.0,
+        100.0 * rand::random::<f64>() - 50.0,
+    );
+    let radius = 100.0 * rand::random::<f64>();
+    let sphere = Sphere::new(center, radius);
+    let u = PI * rand::random::<f64>();
+    let v = 2.0 * PI * rand::random::<f64>();
+    let pt = sphere.subs(u, v);
+    let (u0, v0) = sphere.search_parameter(pt, None, 100).unwrap();
+    assert_near!(Vector2::new(u, v), Vector2::new(u0, v0));
+    let pt = pt
+        + Vector3::new(
+            (0.1 * rand::random::<f64>() + 0.01) * f64::signum(rand::random::<f64>() - 0.5),
+            (0.1 * rand::random::<f64>() + 0.01) * f64::signum(rand::random::<f64>() - 0.5),
+            (0.1 * rand::random::<f64>() + 0.01) * f64::signum(rand::random::<f64>() - 0.5),
+        );
+    assert!(sphere.search_parameter(pt, None, 100).is_none());
+    let (u, v) = sphere.search_nearest_parameter(pt, None, 100).unwrap();
+    assert_near!(
+        sphere.subs(u, v),
+        center + (pt - center).normalize() * radius
+    );
+}
+
+#[test]
+fn search_parameter_test() { (0..10).for_each(|_| exec_search_parameter_test()) }
@@ -0,0 +1,238 @@
+use super::*;
+use std::f64::consts::PI;
+
+impl Torus {
+    /// constructor
+    #[inline(always)]
+    pub fn new(center: Point3, large_radius: f64, small_radius: f64) -> Self {
+        if large_radius <= 0.0 || small_radius <= 0.0 {
+            panic!("radius must be larger than 0");
+        }
+        Self {
+            center,
+            large_radius,
+            small_radius,
+        }
+    }
+
+    /// get center
+    #[inline(always)]
+    pub const fn center(&self) -> Point3 { self.center }
+
+    /// get large radius
+    #[inline(always)]
+    pub const fn large_radius(&self) -> f64 { self.large_radius }
+
+    /// get small radius
+    #[inline(always)]
+    pub const fn small_radius(&self) -> f64 { self.small_radius }
+}
+
+impl ParametricSurface for Torus {
+    type Point = Point3;
+    type Vector = Vector3;
+    #[inline(always)]
+    fn subs(&self, u: f64, v: f64) -> Point3 {
+        let sr = self.small_radius() * Vector2::new(f64::cos(v), f64::sin(v));
+        let lr = (self.large_radius() + sr.x) * Vector2::new(f64::cos(u), f64::sin(u));
+        self.center() + Vector3::new(lr.x, lr.y, sr.y)
+    }
+    #[inline(always)]
+    fn uder(&self, u: f64, v: f64) -> Vector3 {
+        let sr = self.small_radius() * f64::cos(v);
+        let lr = (self.large_radius() + sr) * Vector2::new(f64::cos(u), f64::sin(u));
+        Vector3::new(-lr.y, lr.x, 0.0)
+    }
+    #[inline(always)]
+    fn vder(&self, u: f64, v: f64) -> Vector3 {
+        let sv = self.small_radius() * Vector2::new(-f64::sin(v), f64::cos(v));
+        Vector3::new(sv.x * f64::cos(u), sv.x * f64::sin(u), sv.y)
+    }
+    #[inline(always)]
+    fn uuder(&self, u: f64, v: f64) -> Vector3 {
+        let sr = self.small_radius() * f64::cos(v);
+        let lr = (self.large_radius() + sr) * Vector2::new(f64::cos(u), f64::sin(u));
+        Vector3::new(-lr.x, -lr.y, 0.0)
+    }
+    #[inline(always)]
+    fn uvder(&self, u: f64, v: f64) -> Vector3 {
+        let sr = -self.small_radius() * f64::sin(v);
+        let lr = sr * Vector2::new(f64::cos(u), f64::sin(u));
+        Vector3::new(-lr.y, lr.x, 0.0)
+    }
+    #[inline(always)]
+    fn vvder(&self, u: f64, v: f64) -> Vector3 {
+        let sv = -self.small_radius() * Vector2::new(f64::cos(v), f64::sin(v));
+        Vector3::new(sv.x * f64::cos(u), sv.x * f64::sin(u), sv.y)
+    }
+    #[inline(always)]
+    fn parameter_range(&self) -> (ParameterRange, ParameterRange) {
+        const RANGE: (Bound<f64>, Bound<f64>) = (Bound::Included(0.0), Bound::Excluded(2.0 * PI));
+        (RANGE, RANGE)
+    }
+    #[inline(always)]
+    fn u_period(&self) -> Option<f64> { Some(2.0 * PI) }
+    #[inline(always)]
+    fn v_period(&self) -> Option<f64> { Some(2.0 * PI) }
+}
+
+impl ParametricSurface3D for Torus {
+    fn normal(&self, u: f64, v: f64) -> Vector3 {
+        let sv = Vector2::new(f64::cos(v), f64::sin(v));
+        Vector3::new(sv.x * f64::cos(u), sv.x * f64::sin(u), sv.y)
+    }
+}
+
+impl BoundedSurface for Torus {}
+
+impl SearchParameter<D2> for Torus {
+    type Point = Point3;
+    fn search_parameter<H: Into<SPHint2D>>(
+        &self,
+        point: Point3,
+        _: H,
+        _: usize,
+    ) -> Option<(f64, f64)> {
+        let r = point - self.center();
+        let rxy = Vector2::new(r.x, r.y);
+        let v = f64::asin(f64::clamp(r.z / self.small_radius(), -1.0, 1.0));
+        let minus = rxy.magnitude2() < self.large_radius() * self.large_radius();
+        let v = match (minus, v < 0.0) {
+            (true, _) => PI - v,
+            (false, false) => v,
+            (false, true) => 2.0 * PI + v,
+        };
+        let rxy_n = rxy.normalize();
+        let u = f64::acos(f64::clamp(rxy_n.x, -1.0, 1.0));
+        let u = match rxy_n.y < 0.0 {
+            true => 2.0 * PI - u,
+            false => u,
+        };
+        match self.subs(u, v).near(&point) {
+            true => Some((u, v)),
+            false => None,
+        }
+    }
+}
+
+impl SearchNearestParameter<D2> for Torus {
+    type Point = Point3;
+    fn search_nearest_parameter<H: Into<SPHint2D>>(
+        &self,
+        point: Point3,
+        _: H,
+        _: usize,
+    ) -> Option<(f64, f64)> {
+        let r = point - self.center();
+        let rxy = Vector2::new(r.x, r.y);
+        if rxy.so_small() {
+            return None;
+        }
+        let rxy_n = rxy.normalize();
+        let large_r = self.large_radius() * rxy_n.extend(0.0);
+        let diff = r - large_r;
+        if diff.so_small() {
+            return None;
+        }
+        let small_r = diff.normalize();
+
+        let u = f64::acos(f64::clamp(rxy_n.x, -1.0, 1.0));
+        let u = match rxy_n.y < 0.0 {
+            true => 2.0 * PI - u,
+            false => u,
+        };
+        let v = f64::asin(f64::clamp(small_r.z, -1.0, 1.0));
+        let v = match (small_r.dot(large_r) < 0.0, v < 0.0) {
+            (true, _) => PI - v,
+            (false, false) => v,
+            (false, true) => 2.0 * PI + v,
+        };
+        Some((u, v))
+    }
+}
+
+impl ParameterDivision2D for Torus {
+    fn parameter_division(
+        &self,
+        (urange, vrange): ((f64, f64), (f64, f64)),
+        tol: f64,
+    ) -> (Vec<f64>, Vec<f64>) {
+        let circle = UnitCircle::<Point2>::new();
+        let utol = tol / (self.small_radius() + self.large_radius());
+        let (udiv, _) = circle.parameter_division(urange, utol);
+        let vtol = tol / self.small_radius();
+        let (vdiv, _) = circle.parameter_division(vrange, vtol);
+        (udiv, vdiv)
+    }
+}
+
+#[test]
+fn surface() {
+    use std::ops::RangeBounds;
+    let small_radius = 5.0 * rand::random::<f64>() + 1.0;
+    let large_radius = 2.0 * small_radius + 5.0 * rand::random::<f64>();
+    let torus = Torus::new(Point3::new(0.0, 0.0, 0.0), large_radius, small_radius);
+    const N: usize = 10;
+    const EPS: f64 = 1.0e-2;
+    for i in 0..N {
+        for j in 0..N {
+            let u = 2.0 * PI * i as f64 / N as f64;
+            let v = 2.0 * PI * j as f64 / N as f64;
+            let p = torus.subs(u, v);
+            let r = large_radius * Point3::new(f64::cos(u), f64::sin(u), 0.0);
+            assert_near!(p.distance(r), small_radius);
+            assert_near!(p.z, small_radius * f64::sin(v));
+
+            let uder0 = torus.uder(u, v);
+            let uder1 = (torus.subs(u + EPS, v) - torus.subs(u - EPS, v)) / (2.0 * EPS);
+            assert!((uder0 - uder1).magnitude() < EPS, "{:?} {:?}", uder0, uder1);
+
+            let vder0 = torus.vder(u, v);
+            let vder1 = (torus.subs(u, v + EPS) - torus.subs(u, v - EPS)) / (2.0 * EPS);
+            assert!((vder0 - vder1).magnitude() < EPS, "{:?} {:?}", vder0, vder1);
+
+            let uuder0 = torus.uuder(u, v);
+            let uuder1 = (torus.uder(u + EPS, v) - torus.uder(u - EPS, v)) / (2.0 * EPS);
+            assert!(
+                (uuder0 - uuder1).magnitude() < EPS,
+                "{:?} {:?}",
+                uuder0,
+                uuder1
+            );
+
+            let uvder0 = torus.uvder(u, v);
+            let uvder1 = (torus.vder(u + EPS, v) - torus.vder(u - EPS, v)) / (2.0 * EPS);
+            assert!(
+                (uvder0 - uvder1).magnitude() < EPS,
+                "{:?} {:?}",
+                uvder0,
+                uvder1
+            );
+
+            let vvder0 = torus.vvder(u, v);
+            let vvder1 = (torus.vder(u, v + EPS) - torus.vder(u, v - EPS)) / (2.0 * EPS);
+            assert!(
+                (vvder0 - vvder1).magnitude() < EPS,
+                "{:?} {:?}",
+                vvder0,
+                vvder1
+            );
+
+            let n0 = torus.normal(u, v);
+            let n1 = torus.uder(u, v).cross(torus.vder(u, v)).normalize();
+            assert_near!(n0, n1);
+
+            let (u0, v0) = torus.search_parameter(p, None, 1).unwrap();
+            let (urange, vrange) = torus.parameter_range();
+            assert!(urange.contains(&u0) && vrange.contains(&v0), "{u0}, {v0}");
+            assert_near!(torus.subs(u0, v0), p);
+
+            let deform = (0.8 * rand::random::<f64>() - 0.4) * small_radius;
+            let q = p + deform * n0;
+            let (u0, v0) = torus.search_nearest_parameter(q, None, 1).unwrap();
+            let (urange, vrange) = torus.parameter_range();
+            assert!(urange.contains(&u0) && vrange.contains(&v0), "{u0}, {v0}");
+            assert_near!(torus.subs(u0, v0), p);
+        }
+    }
+}