diff --git a/src/LNLib/Algorithm/KnotVectorUtils.cpp b/src/LNLib/Algorithm/KnotVectorUtils.cpp index c279bfe..9236fc4 100644 --- a/src/LNLib/Algorithm/KnotVectorUtils.cpp +++ b/src/LNLib/Algorithm/KnotVectorUtils.cpp @@ -16,7 +16,36 @@ #include "LNLibExceptions.h" #include -using namespace LNLib; +namespace LNLib +{ + void InsertMidKnotCore(std::vector& unqiueKnotVector, std::vector& insert, int limitNumber) + { + if (insert.size() == limitNumber) + { + return; + } + else + { + double standard = Constants::DoubleEpsilon; + int index = -1; + for (int i = 0; i < unqiueKnotVector.size() - 1; i++) + { + double delta = unqiueKnotVector[i + 1] - unqiueKnotVector[i]; + if (MathUtils::IsGreaterThan(delta, standard)) + { + standard = delta; + index = i; + } + } + double current = unqiueKnotVector[index] + standard / 2.0; + unqiueKnotVector.emplace_back(current); + std::sort(unqiueKnotVector.begin(), unqiueKnotVector.end()); + insert.emplace_back(current); + + InsertMidKnotCore(unqiueKnotVector, insert, limitNumber); + } + } +} int LNLib::KnotVectorUtils::GetContinuity(int degree, const std::vector& knotVector, double knot) { @@ -218,6 +247,17 @@ std::vector> LNLib::KnotVectorUtils::GetInsertedKnotElements return result; } +std::vector LNLib::KnotVectorUtils::GetInsertedKnotElements(int insertKnotsNumber, const std::vector& knotVector) +{ + std::vector unqiueKnotVector = knotVector; + unqiueKnotVector.erase(std::unique(unqiueKnotVector.begin(), unqiueKnotVector.end()), unqiueKnotVector.end()); + + std::vector insert; + InsertMidKnotCore(unqiueKnotVector, insert, insertKnotsNumber); + std::sort(insert.begin(), insert.end()); + return insert; +} + bool LNLib::KnotVectorUtils::IsUniform(const std::vector& knotVector) { auto map = GetKnotMultiplicityMap(knotVector); diff --git a/src/LNLib/Geometry/Surface/NurbsSurface.cpp b/src/LNLib/Geometry/Surface/NurbsSurface.cpp index a3fb4d1..8402689 100644 --- a/src/LNLib/Geometry/Surface/NurbsSurface.cpp +++ b/src/LNLib/Geometry/Surface/NurbsSurface.cpp @@ -90,6 +90,23 @@ namespace LNLib } return ind; } + + XYZ LocalToWorld(const XYZ& localOrigin, const XYZ& localXdir, const XYZ& localYdir, const XYZ& localZdir, const XYZ& worldPoint) + { + double x = localXdir.GetX() * worldPoint.GetX() + localYdir.GetX() * worldPoint.GetY() + localZdir.GetX() * worldPoint.GetZ(); + double y = localXdir.GetY() * worldPoint.GetX() + localYdir.GetY() * worldPoint.GetY() + localZdir.GetY() * worldPoint.GetZ(); + double z = localXdir.GetZ() * worldPoint.GetX() + localYdir.GetZ() * worldPoint.GetY() + localZdir.GetZ() * worldPoint.GetZ(); + return XYZ(x, y, z) + localOrigin; + } + + XYZ WorldToLocal(const XYZ& localOrigin, const XYZ& localXdir, const XYZ& localYdir, const XYZ& localZdir, const XYZ& worldPoint) + { + XYZ traslation = (worldPoint - localOrigin); + double x = localXdir.GetX() * traslation.GetX() + localXdir.GetY() * traslation.GetY() + localXdir.GetZ() * traslation.GetZ(); + double y = localYdir.GetX() * traslation.GetX() + localYdir.GetY() * traslation.GetY() + localYdir.GetZ() * traslation.GetZ(); + double z = localZdir.GetX() * traslation.GetX() + localZdir.GetY() * traslation.GetY() + localZdir.GetZ() * traslation.GetZ(); + return XYZ(x, y, z); + } } void LNLib::NurbsSurface::Check(const LN_NurbsSurface& surface) @@ -2226,6 +2243,7 @@ void LNLib::NurbsSurface::CreateLoftSurface(const std::vector& se int degreeU = degree_max; int degreeV = degree_max; + int K = size - 1; std::vector vl(size); vl[0] = 0; @@ -2239,7 +2257,7 @@ void LNLib::NurbsSurface::CreateLoftSurface(const std::vector& se for (int i = 0; i <= n; i++) { double delta = curvesControlPoints[k][i].ToXYZ(true).Distance(curvesControlPoints[k - 1][i].ToXYZ(true)); - + std::vector tempPoints(size); for (int j = 0; j <= K; j++) { @@ -2249,7 +2267,7 @@ void LNLib::NurbsSurface::CreateLoftSurface(const std::vector& se double di = Interpolation::GetTotalChordLength(tempPoints); sum += delta / di; } - + vl[k] = vl[k - 1] + (1.0 / (n + 1)) * sum; } std::vector knotVectorV = Interpolation::AverageKnotVector(degreeV, vl); @@ -2315,33 +2333,86 @@ void LNLib::NurbsSurface::CreateGeneralizedTranslationalSweepSurface(const LN_Nu surface.ControlPoints = controlPoints; } -void LNLib::NurbsSurface::CreateSweepSurface(const LN_NurbsCurve& path, const std::vector& profiles, LN_NurbsSurface& surface) +void LNLib::NurbsSurface::CreateSweepSurface(const LN_NurbsCurve& profile, const LN_NurbsCurve& trajectory, int minimumProfiles, LN_NurbsSurface& surface) { - int profilesSize = profiles.size(); - double path_min = path.KnotVector[0]; - double path_max = path.KnotVector[path.KnotVector.size() - 1]; - double delta = (path_max - path_min) / profilesSize; + // Determine number of sections for lofting. + int q = trajectory.Degree; + const std::vector& trajectoryKnotVector = trajectory.KnotVector; + int ktv = trajectoryKnotVector.size(); + int K = minimumProfiles; + int nsect = K + 1; + + LN_NurbsCurve trajectoryCopy = trajectory; + if (ktv <= nsect + q) + { + // Insert knots to widest spans + // to reach the required minimum number of sections. + int m = nsect + q - ktv + 1; + std::vector insertedElements = KnotVectorUtils::GetInsertedKnotElements(m, trajectoryKnotVector); + NurbsCurve::RefineKnotVector(trajectory, insertedElements, trajectoryCopy); + } + else + { + if (ktv > nsect + q + 1) + { + nsect = ktv - q - 1; + } + } + + // Determine lofting profile positions. + std::vector knotVectorV = trajectoryCopy.KnotVector; + auto minmax = std::minmax_element(knotVectorV.begin(), knotVectorV.end()); + + std::vector vk(nsect); + vk[0] = *minmax.first; + vk[nsect - 1] = *minmax.second; + + for (int k = 1; k < nsect - 1; k++) + { + double sum = 0.0; + for (int i = 1; i <= q; i++) + { + sum += knotVectorV[k + i]; + } + vk[k] = sum / (double)q; + } + + // Compute trajectory normals. + std::vector Bv = NurbsCurve::ProjectNormal(trajectoryCopy); + const std::vector& profileControlPoints = profile.ControlPoints; + int profileCpSize = profileControlPoints.size(); - std::vector sections(profilesSize); - for (int i = 0; i < profilesSize; i++) + // for each lofting section + std::vector sections(nsect); + for (int k = 0; k < nsect; k++) { - double param = path_min + i * delta; - XYZ point = NurbsCurve::GetPointOnCurve(path, param); - std::vector ders = NurbsCurve::ComputeRationalCurveDerivatives(profiles[i], 1, param); - XYZ tangent = ders[1]; - Matrix4d transform = Matrix4d::CreateTranslation(point); - if (!tangent.IsZero()) + // Build local coordinate system for the section. + XYZ zdir = NurbsCurve::ComputeRationalCurveDerivatives(trajectoryCopy, 1, vk[k])[1].Normalize(); + int spanIndex = Polynomials::GetKnotSpanIndex(trajectoryCopy.Degree, trajectoryCopy.KnotVector, vk[k]); + const XYZ& xdir = Bv[spanIndex]; + XYZ ydir = zdir.CrossProduct(xdir).Normalize(); + + // Transform the input profile to the lofting position. + std::vector transformedControlPoints(profileCpSize); + for (int i = 0; i < profileCpSize; i++) { - XYZ binormal = NurbsCurve::Normal(profiles[i], CurveNormal::Binormal, param); - double rad = binormal.AngleTo(tangent); - Matrix4d rotation = Matrix4d::CreateRotationAtPoint(point, binormal, rad); - transform = transform.Multiply(rotation); + XYZW wp = profileControlPoints[i]; + double w = wp.GetW(); + XYZ p = wp.ToXYZ(true); + + XYZ transformed = WorldToLocal(NurbsCurve::GetPointOnCurve(trajectoryCopy, vk[k]), xdir, ydir, zdir, p); + transformedControlPoints[i] = XYZW(transformed, w); } - LN_NurbsCurve newProfile; - NurbsCurve::CreateTransformed(profiles[i], transform, newProfile); - sections.emplace_back(newProfile); + + // Make section curve. + // note: std::vector::swap relinks pointers, which is faster than copying. + LN_NurbsCurve& section = sections[k]; + section.Degree = profile.Degree; + section.KnotVector = profile.KnotVector; + section.ControlPoints.swap(transformedControlPoints); } + // Do the lofting. CreateLoftSurface(sections, surface); } diff --git a/src/LNLib/include/KnotVectorUtils.h b/src/LNLib/include/KnotVectorUtils.h index 9fe9f69..29f421a 100644 --- a/src/LNLib/include/KnotVectorUtils.h +++ b/src/LNLib/include/KnotVectorUtils.h @@ -60,6 +60,13 @@ namespace LNLib /// static std::vector> GetInsertedKnotElements(const std::vector>& knotVectors); + /// + /// The NURBS Book 2nd Edition Page476 + /// Insert [insertKnotsNumber] knots into [knotVector] + /// by insert at the midpoint of the longest span will do. + /// + static std::vector GetInsertedKnotElements(int insertKnotsNumber, const std::vector& knotVector); + /// /// The NURBS Book 2nd Edition Page572 /// diff --git a/src/LNLib/include/NurbsSurface.h b/src/LNLib/include/NurbsSurface.h index acfe150..323e14f 100644 --- a/src/LNLib/include/NurbsSurface.h +++ b/src/LNLib/include/NurbsSurface.h @@ -216,11 +216,12 @@ namespace LNLib static void CreateGeneralizedTranslationalSweepSurface(const LN_NurbsCurve& profile, const LN_NurbsCurve& trajectory, LN_NurbsSurface& surface); /// - /// The NURBS Book 2nd Edition Page472 - /// Algorithm A10.2 - /// Create Sweep Surface. + /// The NURBS Book 2nd Edition Page475 + /// Algorithm A10.1 + /// Create sweep surface by trajectory interpolated. + /// Profile must lie on XOY plane. /// - static void CreateSweepSurface(const LN_NurbsCurve& path, const std::vector& profiles, LN_NurbsSurface& surface); + static void CreateSweepSurface(const LN_NurbsCurve& profile, const LN_NurbsCurve& trajectory, int minimumProfiles, LN_NurbsSurface& surface); /// /// The NURBS Book 2nd Edition Page494 @@ -241,7 +242,6 @@ namespace LNLib /// Create Coons Surface. /// The difference between Coons and Gordon is that Coons Surface is created by 4 anti-clock curves. /// The coons surface is the special case of Gordon Surface. - /// /// static void CreateCoonsSurface(const LN_NurbsCurve& leftCurve, const LN_NurbsCurve& bottomCurve, const LN_NurbsCurve& rightCurve, const LN_NurbsCurve& topCurve, LN_NurbsSurface& surface); diff --git a/tests/T_AdvancedSurface.cpp b/tests/T_AdvancedSurface.cpp index 3eae9b2..84d8ebb 100644 --- a/tests/T_AdvancedSurface.cpp +++ b/tests/T_AdvancedSurface.cpp @@ -181,4 +181,9 @@ TEST(Test_AdvancedSurface, CreateGeneralizedTranslationalSweepSurface) double expectedArea = (endRad - startRad) * radius * height; double area = NurbsSurface::ApproximateArea(surface); EXPECT_NEAR(area, expectedArea, 1e-4); + + + NurbsSurface::CreateSweepSurface(profile, trajectory, 5, surface); + area = NurbsSurface::ApproximateArea(surface); + EXPECT_NEAR(area, expectedArea, 1e-4); }