// -------------------------------------------------------------------- // Ipe geometry primitives // -------------------------------------------------------------------- /* This file is part of the extensible drawing editor Ipe. Copyright (C) 1993-2004 Otfried Cheong Ipe is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. As a special exception, you have permission to link Ipe with the CGAL library and distribute executables, as long as you follow the requirements of the Gnu General Public License in regard to all of the software in the executable aside from CGAL. Ipe is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with Ipe; if not, you can find it at "http://www.gnu.org/copyleft/gpl.html", or write to the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */ /*! \defgroup geo Ipe Geometry \brief Geometric primitives for Ipe. The IpeGeo module provides a few classes for constant-size geometric primitives, such as vector, axis-aligned rectangles, lines, rays, line segments, etc. */ #include "ipegeo.h" inline double Sq(double x) { return x * x; } // -------------------------------------------------------------------- /*! \class IpeAngle \ingroup geo \brief A double that's an angle. An IpeAngle is really nothing more than a double. Having a separate type is sometimes useful, for instance in the IpeVector constructor, and this class serves as the right place for a few utility functions. It also makes it clear whether a value is in radians or in degrees. */ double IpeAngle::Degrees() const { return (iAlpha / IpePi * 180.0); } //! Normalize the value to the range lowlimit .. lowlimit + 2 pi. /*! This IpeAngle object is modified, a copy is returned. */ IpeAngle IpeAngle::Normalize(double lowlimit) { while (iAlpha >= lowlimit + 2.0*IpePi) iAlpha -= 2.0 * IpePi; while (iAlpha < lowlimit) iAlpha += 2.0 * IpePi; return *this; } /*! When considering the positively oriented circle arc from angle \a small to \a large, does it cover this angle? */ bool IpeAngle::LiesBetween(IpeAngle small, IpeAngle large) const { large.Normalize(iAlpha); small.Normalize(large.iAlpha - 2.0 * IpePi); return (iAlpha >= small.iAlpha); } // -------------------------------------------------------------------- /*! \class IpeVector \ingroup geo \brief Two-dimensional vector. Unlike some other libraries, I don't make a difference between points and vectors. */ //! Construct a unit vector with this direction. IpeVector::IpeVector(IpeAngle alpha) { iX = cos(alpha); iY = sin(alpha); } //! Return angle of the vector (with positive x-direction). /*! Returns zero for the zero vector. */ IpeAngle IpeVector::Angle() const { if (iX == 0.0 && iY == 0.0) return IpeAngle(0.0); else return IpeAngle(atan2(iY, iX)); } //! The origin (zero vector). IpeVector IpeVector::Zero = IpeVector(0.0, 0.0); double IpeVector::Len() const { return sqrt(SqLen()); } //! Return this vector normalized (with length one). IpeVector IpeVector::Normalized() const { double len = SqLen(); if (len == 0.0 || len == 1.0) return *this; return (1.0/sqrt(len)) * (*this); } //! Return this vector turned 90 degrees to the left. IpeVector IpeVector::Orthogonal() const { return IpeVector(-iY, iX); } /*! Normalizes this vector into \a unit and returns length. If this is the zero vector, \a unit is set to (1,0). */ double IpeVector::Factorize(IpeVector &unit) const { double len = SqLen(); if (len == 0.0) { unit = IpeVector(1,0); return len; } if (len == 1.0) { unit = *this; return len; } len = sqrt(len); unit = (1.0 / len) * (*this); return len; } //! Output operator for IpeVector. IpeStream &operator<<(IpeStream &stream, const IpeVector &rhs) { return stream << rhs.iX << " " << rhs.iY; } // -------------------------------------------------------------------- /*! \class IpeRect \ingroup geo \brief Axis-parallel rectangle (which can be empty) */ //! Create rectangle containing points \a c1 and \a c2. IpeRect::IpeRect(const IpeVector &c1, const IpeVector &c2) : iMin(1,0), iMax(-1,0) { AddPoint(c1); AddPoint(c2); } //! Does (closed) rectangle contain the point? bool IpeRect::Contains(const IpeVector &rhs) const { // this correctly handles empty this return (iMin.iX <= rhs.iX && rhs.iX <= iMax.iX && iMin.iY <= rhs.iY && rhs.iY <= iMax.iY); } //! Does rectangle contain other rectangle? bool IpeRect::Contains(const IpeRect &rhs) const { if (rhs.IsEmpty()) return true; if (IsEmpty()) return false; return (iMin.iX <= rhs.iMin.iX && rhs.iMax.iX <= iMax.iX && iMin.iY <= rhs.iMin.iY && rhs.iMax.iY <= iMax.iY); } //! Does rectangle intersect other rectangle? bool IpeRect::Intersects(const IpeRect &rhs) const { if (IsEmpty() || rhs.IsEmpty()) return false; return (iMin.iX <= rhs.iMax.iX && rhs.iMin.iX <= iMax.iX && iMin.iY <= rhs.iMax.iY && rhs.iMin.iY <= iMax.iY); } //! Enlarge rectangle to contain point. void IpeRect::AddPoint(const IpeVector &rhs) { if (IsEmpty()) { iMin = rhs; iMax = rhs; } else { if (rhs.iX > iMax.iX) iMax.iX = rhs.iX; else if (rhs.iX < iMin.iX) iMin.iX = rhs.iX; if (rhs.iY > iMax.iY) iMax.iY = rhs.iY; else if (rhs.iY < iMin.iY) iMin.iY = rhs.iY; } } //! Enlarge rectangle to contain rhs rectangle. void IpeRect::AddRect(const IpeRect &rhs) { if (IsEmpty()) { iMin = rhs.iMin; iMax = rhs.iMax; } else if (!rhs.IsEmpty()) { if (rhs.iMax.iX > iMax.iX) iMax.iX = rhs.iMax.iX; if (rhs.iMin.iX < iMin.iX) iMin.iX = rhs.iMin.iX; if (rhs.iMax.iY > iMax.iY) iMax.iY = rhs.iMax.iY; if (rhs.iMin.iY < iMin.iY) iMin.iY = rhs.iMin.iY; } } /*! Returns false if the distance between the box and v is smaller than \a bound. Often returns true if their distance is larger than \a bound. */ bool IpeRect::CertainClearance(const IpeVector &v, double bound) const { return ((iMin.iX - v.iX) >= bound || (v.iX - iMax.iX) >= bound || (iMin.iY - v.iY) >= bound || (v.iY - iMax.iY) >= bound); } //! Output operator for IpeRect. IpeStream &operator<<(IpeStream &stream, const IpeRect &rhs) { return stream << rhs.Min() << " " << rhs.Max(); } // -------------------------------------------------------------------- /*! \class IpeLine \ingroup geo \brief A directed line. */ //! Construct a line from \a p with direction \a dir. /*! Asserts unit length of \a dir. */ IpeLine::IpeLine(const IpeVector &p, const IpeVector &dir) { assert(Sq(dir.SqLen() - 1.0) < 1e-10); iP = p; iDir = dir; } //! Construct a line through two points. IpeLine IpeLine::Through(const IpeVector &p, const IpeVector &q) { assert(q != p); return IpeLine(p, (q - p).Normalized()); } //! Result is > 0, = 0, < 0 if point lies to the left, on, to the right. double IpeLine::Side(const IpeVector &p) const { return Dot(Normal(), p - iP); } //! Returns distance between line and \a v. double IpeLine::Distance(const IpeVector &v) const { IpeVector diff = v - iP; return (diff - Dot(diff, iDir) * iDir).Len(); } inline double cross(const IpeVector &v1, const IpeVector &v2) { return v1.iX * v2.iY - v1.iY * v2.iX; } static bool Intersect(double &lambda, const IpeLine &l, const IpeLine &m) { double denom = cross(m.Dir(), l.Dir()); if (denom == 0.0) return false; lambda = cross(l.iP - m.iP, m.Dir()) / denom; return true; } //! Does this line intersect \a line? If so, computes intersection point. bool IpeLine::Intersects(const IpeLine &line, IpeVector &pt) { double lambda; if (Intersect(lambda, *this, line)) { pt = iP + lambda * iDir; return true; } return false; } //! Orthogonally project point \a v onto the line. IpeVector IpeLine::Project(const IpeVector &v) const { double dx = Dot(iDir, v - iP); return iP + dx * iDir; } // -------------------------------------------------------------------- /*! \class IpeSegment \ingroup geo \brief A directed line segment. */ /*! Returns distance between segment and point \a v, but may just return \a bound when its larger than \a bound. */ double IpeSegment::Distance(const IpeVector &v, double bound) const { if (IpeRect(iP, iQ).CertainClearance(v, bound)) return bound; return Distance(v); } /*! Returns distance between segment and point \a v */ double IpeSegment::Distance(const IpeVector &v) const { IpeVector dir = iQ - iP; IpeVector udir; double len = dir.Factorize(udir); double dx = Dot(udir, v - iP); if (dx <= 0) return (v - iP).Len(); if (dx >= len) return (v - iQ).Len(); return (v - (iP + dx * udir)).Len(); } /*! Project point \a v orthogonally on segment. Returns false if the point falls outside the segment. */ bool IpeSegment::Project(const IpeVector &v, IpeVector &projection) const { IpeVector dir = iQ - iP; IpeVector udir; double len = dir.Factorize(udir); double dx = Dot(udir, v - iP); if (dx <= 0 || dx >= len) return false; projection = iP + dx * udir; return true; } //! Compute intersection point. Return \c false if segs don't intersect. bool IpeSegment::Intersects(const IpeSegment &seg, IpeVector &pt) const { if (iP == iQ || seg.iP == seg.iQ) return false; if (!IpeRect(iP, iQ).Intersects(IpeRect(seg.iP, seg.iQ))) return false; if (!Line().Intersects(seg.Line(), pt)) return false; // have intersection point, check whether it's on both segments. IpeVector dir = iQ - iP; IpeVector dir1 = seg.iQ - seg.iP; return (Dot(pt - iP, dir) >= 0 && Dot(pt - iQ, dir) <= 0 && Dot(pt - seg.iP, dir1) >= 0 && Dot(pt - seg.iQ, dir1) <= 0); } //! Compute intersection point. Return \c false if no intersection. bool IpeSegment::Intersects(const IpeLine &l, IpeVector &pt) const { if (!Line().Intersects(l, pt)) return false; // have intersection point, check whether it's on both segments. IpeVector dir = iQ - iP; return (Dot(pt - iP, dir) >= 0 && Dot(pt - iQ, dir) <= 0); } // -------------------------------------------------------------------- /*! \class IpeLinear \ingroup geo \brief Linear transformation in the plane (2x2 matrix). */ //! Create matrix representing a rotation by angle. IpeLinear::IpeLinear(IpeAngle angle) { iA[0] = cos(angle); iA[1] = sin(angle); iA[2] = -iA[1]; iA[3] = iA[0]; } //! Parse string. IpeLinear::IpeLinear(IpeString str) { IpeLex lex(str); lex >> iA[0] >> iA[1] >> iA[2] >> iA[3]; } //! Return inverse. IpeLinear IpeLinear::Inverse() const { double t = 1.0/(iA[0]*iA[3]-iA[1]*iA[2]); return IpeLinear(iA[3]*t, -iA[1]*t, -iA[2]*t, iA[0]*t); } //! Output operator for IpeLinear. IpeStream &operator<<(IpeStream &stream, const IpeLinear &rhs) { return stream << rhs.iA[0] << " " << rhs.iA[1] << " " << rhs.iA[2] << " " << rhs.iA[3]; } // -------------------------------------------------------------------- /*! \class IpeMatrix \ingroup geo \brief Homogeneous transformation in the plane. */ //! Parse string. IpeMatrix::IpeMatrix(IpeString str) { IpeLex lex(str); lex >> iA[0] >> iA[1] >> iA[2] >> iA[3] >> iA[4] >> iA[5]; } //! Return inverse. IpeMatrix IpeMatrix::Inverse() const { double t = iA[0]*iA[3]-iA[1]*iA[2]; assert(t != 0); t = 1.0/t; return IpeMatrix(iA[3]*t, -iA[1]*t, -iA[2]*t, iA[0]*t, (iA[2]*iA[5]-iA[3]*iA[4])*t, -(iA[0]*iA[5]-iA[1]*iA[4])*t); } //! Output operator for IpeMatrix. IpeStream &operator<<(IpeStream &stream, const IpeMatrix &rhs) { return stream << rhs.iA[0] << " " << rhs.iA[1] << " " << rhs.iA[2] << " " << rhs.iA[3] << " " << rhs.iA[4] << " " << rhs.iA[5]; } // -------------------------------------------------------------------- /*! \class IpeBezier \ingroup geo \brief A cubic Bezier spline. */ inline IpeVector midpoint(const IpeVector& p, const IpeVector& q) { return 0.5 * (p + q); } inline IpeVector thirdpoint(const IpeVector& p, const IpeVector& q) { return (1.0/3.0) * ((2 * p) + q); } //! Return point on curve with parameter \a t (from 0.0 to 1.0). IpeVector IpeBezier::Point(double t) const { double t1 = 1.0 - t; return t1 * t1 * t1 * iV[0] + 3 * t * t1 * t1 * iV[1] + 3 * t * t * t1 * iV[2] + t * t * t * iV[3]; } /*! Returns true if the Bezier curve is nearly identical to the line segment iV[0]..iV[3]. */ bool IpeBezier::Straight(double precision) const { if (iV[0] == iV[3]) { return ((iV[1] - iV[0]).Len() < precision && (iV[1] - iV[0]).Len() < precision); } else { IpeLine l = IpeLine::Through(iV[0], iV[3]); double d1 = l.Distance(iV[1]); double d2 = l.Distance(iV[2]); return (d1 < precision && d2 < precision); } } //! Subdivide this Bezier curve in the middle. void IpeBezier::Subdivide(IpeBezier &l, IpeBezier &r) const { IpeVector h; l.iV[0] = iV[0]; l.iV[1] = 0.5 * (iV[0] + iV[1]); h = 0.5 * (iV[1] + iV[2]); l.iV[2] = 0.5 * (l.iV[1] + h); r.iV[2] = 0.5 * (iV[2] + iV[3]); r.iV[1] = 0.5 * (h + r.iV[2]); r.iV[0] = 0.5 * (l.iV[2] + r.iV[1]); l.iV[3] = r.iV[0]; r.iV[3] = iV[3]; } //! Approximate by a polygonal chain. /*! \a result must be empty when calling this. */ void IpeBezier::Approximate(double precision, std::vector &result) const { if (Straight(precision)) { result.push_back(iV[3]); } else { IpeBezier l, r; Subdivide(l, r); l.Approximate(precision, result); r.Approximate(precision, result); } } //! Convert a quadratic Bezier-spline to a cubic one. /*! The quadratic Bezier-spline with control points p0, p1, p2 is identical to the cubic Bezier-spline with control points q0 = p0, q1 = (2p1 + p0)/3, q2 = (2p1 + p2)/3, q3 = p2. */ IpeBezier IpeBezier::QuadBezier(const IpeVector &p0, const IpeVector &p1, const IpeVector &p2) { IpeVector q1 = thirdpoint(p1, p0); IpeVector q2 = thirdpoint(p1, p2); return IpeBezier(p0, q1, q2, p2); } //! Convert a uniform cubic B-spline to a series of Bezier splines. /*! First and last control point are given multiplicity 3. Bezier splines are appended to \a result. */ void IpeBezier::Spline(int n, const IpeVector *v, std::vector &result) { IpeVector p0, p1, p2, p3, q0, q1, q2, q3; // First segment (p1 = p2 = p0 => q1 = q2 = q0 = p0) p0 = v[0]; p3 = v[1]; q3 = midpoint(thirdpoint(p0, p3), p0); result.push_back(IpeBezier(p0, p0, p0, q3)); if (n > 2) { // Second segment p1 = v[0]; p2 = v[1]; p3 = v[2]; q0 = q3; // from previous q1 = thirdpoint(p1, p2); q2 = thirdpoint(p2, p1); q3 = midpoint(thirdpoint(p2, p3), q2); result.push_back(IpeBezier(q0, q1, q2, q3)); // create n - 3 segments for (int i = 0; i < n - 3; ++i) { p0 = v[i]; p1 = v[i + 1]; p2 = v[i + 2]; p3 = v[i + 3]; q0 = q3; // from previous // q0 = midpoint(thirdpoint(p1, p0), q1); // the real formula q1 = thirdpoint(p1, p2); q2 = thirdpoint(p2, p1); q3 = midpoint(thirdpoint(p2, p3), q2); result.push_back(IpeBezier(q0, q1, q2, q3)); } } // Second to last segment p1 = v[n-2]; p2 = v[n-1]; p3 = v[n-1]; q0 = q3; // from previous q1 = thirdpoint(p1, p2); q2 = thirdpoint(p2, p1); q3 = midpoint(p3, q2); result.push_back(IpeBezier(q0, q1, q2, q3)); // Last segment (p1 = p2 = p3 => q1 = q2 = q3 = p3) result.push_back(IpeBezier(q3, p3, p3, p3)); } //! Convert a closed uniform cubic B-spline to a series of Bezier splines. /*! Bezier splines are appended to \a result. */ void IpeBezier::ClosedSpline(int n, const IpeVector *v, std::vector &result) { for (int i = 0; i < n; ++i) { IpeVector p0 = v[i]; IpeVector p1 = v[(i + 1) % n]; IpeVector p2 = v[(i + 2) % n]; IpeVector p3 = v[(i + 3) % n]; IpeVector q1 = thirdpoint(p1, p2); IpeVector q2 = thirdpoint(p2, p1); IpeVector q0 = midpoint(thirdpoint(p1, p0), q1); IpeVector q3 = midpoint(thirdpoint(p2, p3), q2); result.push_back(IpeBezier(q0, q1, q2, q3)); } } //! Return distance to Bezier spline. /*! But may just return \a bound if actual distance is larger. */ double IpeBezier::Distance(const IpeVector &v, double bound) { IpeRect box; box.AddPoint(iV[0]); box.AddPoint(iV[1]); box.AddPoint(iV[2]); box.AddPoint(iV[3]); if (box.CertainClearance(v, bound)) return bound; std::vector approx; Approximate(3.0, approx); IpeVector cur = iV[0]; double d = bound; double d1; for (std::vector::const_iterator it = approx.begin(); it != approx.end(); ++it) { if ((d1 = IpeSegment(cur, *it).Distance(v, d)) < d) d = d1; cur = *it; } return d; } //! Return a tight bounding box (accurate to within 0.5). IpeRect IpeBezier::BBox() const { IpeRect box(iV[0]); std::vector approx; Approximate(0.5, approx); for (std::vector::const_iterator it = approx.begin(); it != approx.end(); ++it) { box.AddPoint(*it); } return IpeRect(box.Min() - IpeVector(0.5, 0.5), box.Max() + IpeVector(0.5, 0.5)); } bool IpeBezier::Intersects(const IpeLine &, IpeVector &) const { // XXX return false; } bool IpeBezier::Intersects(const IpeSegment &, IpeVector &) const { // XXX return false; } // -------------------------------------------------------------------- /*! \class IpeArc \ingroup geo \brief An arc of an ellipse. This object can also describe a full ellipse if iAlpha = iBeta = 0. */ //! Construct arc for ellipse defined by m, from begp to endp. IpeArc::IpeArc(const IpeMatrix &m, const IpeVector &begp, const IpeVector &endp) { IpeMatrix inv = m.Inverse(); iAlpha = (inv * begp).Angle(); iBeta = (inv * endp).Angle(); iM = m; } //! This doesn't really compute the distance, but a reasonable approximation. double IpeArc::Distance(const IpeVector &v, double bound) { IpeVector pos; IpeAngle angle; return Distance(v, bound, pos, angle); } /*! Like Distance(), but sets pos to point on arc and angle to its angle in arc coordinates. \a angle and \a pos are not modified if result is larger than bound. */ double IpeArc::Distance(const IpeVector &v, double bound, IpeVector &pos, IpeAngle &angle) { IpeMatrix inv1 = iM.Inverse(); IpeVector v1 = inv1 * v; IpeVector pos1 = iM * v1.Normalized(); double d = (v - pos1).Len(); if (iAlpha != 0.0 || iBeta != 0.0) { if (d < bound && v1.Angle().LiesBetween(iAlpha, iBeta)) { bound = d; pos = pos1; angle = v1.Angle(); } pos1 = iM * IpeVector(iAlpha); d = (v - pos1).Len(); if (d < bound) { bound = d; pos = pos1; angle = iAlpha; } pos1 = iM * IpeVector(iBeta); d = (v - pos1).Len(); if (d < bound) { bound = d; pos = pos1; angle = iBeta; } } else { // full ellipse if (d < bound) { bound = d; pos = pos1; angle = v1.Angle(); } } return bound; } //! Return a tight bounding box. IpeRect IpeArc::BBox() const { IpeRect box; // add begin and end point box.AddPoint(iM * IpeVector(iAlpha)); box.AddPoint(iM * IpeVector(iBeta)); IpeLinear inv = iM.Linear().Inverse(); bool ell = (iAlpha == 0.0 && iBeta == 0.0); IpeAngle alpha = (inv * IpeVector(0,1)).Angle() - IpeHalfPi; if (ell || alpha.LiesBetween(iAlpha, iBeta)) box.AddPoint(iM * IpeVector(alpha)); alpha = (inv * IpeVector(0,-1)).Angle() - IpeHalfPi; if (ell || alpha.LiesBetween(iAlpha, iBeta)) box.AddPoint(iM * IpeVector(alpha)); alpha = (inv * IpeVector(1,0)).Angle() - IpeHalfPi; if (ell || alpha.LiesBetween(iAlpha, iBeta)) box.AddPoint(iM * IpeVector(alpha)); alpha = (inv * IpeVector(-1,0)).Angle() - IpeHalfPi; if (ell || alpha.LiesBetween(iAlpha, iBeta)) box.AddPoint(iM * IpeVector(alpha)); return box; } bool IpeArc::Intersects(const IpeLine &, IpeVector &) const { /* IpeMatrix inv = iM.Inverse(); IpeVector p = inv * l.iP; IpeVector q = inv * (l.iP + l.Dir()); IpeLine l1 = IpeLine::Through(p, q); */ // XXX return false; } bool IpeArc::Intersects(const IpeSegment &, IpeVector &) const { // XXX return false; } // --------------------------------------------------------------------