Merge pull request BIMCoderLiang#23 from BIMCoderLiang/sweepsurface

imporve: reconstruct create sweep surface.

src/LNLib/Algorithm/KnotVectorUtils.cpp

@@ -16,7 +16,36 @@
 #include "LNLibExceptions.h"
 #include <algorithm>
 
-using namespace LNLib;
+namespace LNLib
+{
+	void InsertMidKnotCore(std::vector<double>& unqiueKnotVector, std::vector<double>& 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<double>& knotVector, double knot)
 {
@@ -218,6 +247,17 @@ std::vector<std::vector<double>> LNLib::KnotVectorUtils::GetInsertedKnotElements
 	return result;
 }
 
+std::vector<double> LNLib::KnotVectorUtils::GetInsertedKnotElements(int insertKnotsNumber, const std::vector<double>& knotVector)
+{
+	std::vector<double> unqiueKnotVector = knotVector;
+	unqiueKnotVector.erase(std::unique(unqiueKnotVector.begin(), unqiueKnotVector.end()), unqiueKnotVector.end());
+
+	std::vector<double> insert;
+	InsertMidKnotCore(unqiueKnotVector, insert, insertKnotsNumber);
+	std::sort(insert.begin(), insert.end());
+	return insert;
+}
+
 bool LNLib::KnotVectorUtils::IsUniform(const std::vector<double>& knotVector)
 {
 	auto map = GetKnotMultiplicityMap(knotVector);

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<LN_NurbsCurve>& se
 	int degreeU = degree_max;
 	int degreeV = degree_max;
 
+
 	int K = size - 1;
 	std::vector<double> vl(size);
 	vl[0] = 0;
@@ -2239,7 +2257,7 @@ void LNLib::NurbsSurface::CreateLoftSurface(const std::vector<LN_NurbsCurve>& se
 		for (int i = 0; i <= n; i++)
 		{
 			double delta = curvesControlPoints[k][i].ToXYZ(true).Distance(curvesControlPoints[k - 1][i].ToXYZ(true));
-			
+
 			std::vector<XYZ> tempPoints(size);
 			for (int j = 0; j <= K; j++)
 			{
@@ -2249,7 +2267,7 @@ void LNLib::NurbsSurface::CreateLoftSurface(const std::vector<LN_NurbsCurve>& se
 			double di = Interpolation::GetTotalChordLength(tempPoints);
 			sum += delta / di;
 		}
-		
+
 		vl[k] = vl[k - 1] + (1.0 / (n + 1)) * sum;
 	}
 	std::vector<double> 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<LN_NurbsCurve>& 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<double>& 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<double> 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<double> knotVectorV = trajectoryCopy.KnotVector;
+	auto minmax = std::minmax_element(knotVectorV.begin(), knotVectorV.end());
+
+	std::vector<double> 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<XYZ> Bv =  NurbsCurve::ProjectNormal(trajectoryCopy);
+	const std::vector<XYZW>& profileControlPoints = profile.ControlPoints;
+	int profileCpSize = profileControlPoints.size();
 
-	std::vector<LN_NurbsCurve> sections(profilesSize);
-	for (int i = 0; i < profilesSize; i++)
+	// for each lofting section
+	std::vector<LN_NurbsCurve> sections(nsect);
+	for (int k = 0; k < nsect; k++)
 	{
-		double param = path_min + i * delta;
-		XYZ point = NurbsCurve::GetPointOnCurve(path, param);
-		std::vector<XYZ> 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<XYZW> 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);
 }
 

src/LNLib/include/KnotVectorUtils.h

@@ -60,6 +60,13 @@ namespace LNLib
 		/// </summary>
 		static std::vector<std::vector<double>> GetInsertedKnotElements(const std::vector<std::vector<double>>& knotVectors);
 
+		/// <summary>
+		/// The NURBS Book 2nd Edition Page476
+		/// Insert [insertKnotsNumber] knots into [knotVector] 
+		/// by insert at the midpoint of the longest span will do.
+		/// </summary>
+		static std::vector<double> GetInsertedKnotElements(int insertKnotsNumber, const std::vector<double>& knotVector);
+
 		/// <summary>
 		/// The NURBS Book 2nd Edition Page572
 		/// </summary>

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);
 
 		/// <summary>
-		/// 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.
 		/// </summary>
-		static void CreateSweepSurface(const LN_NurbsCurve& path, const std::vector<LN_NurbsCurve>& profiles, LN_NurbsSurface& surface);
+		static void CreateSweepSurface(const LN_NurbsCurve& profile, const LN_NurbsCurve& trajectory, int minimumProfiles, LN_NurbsSurface& surface);
 
 		/// <summary>
 		/// 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.
-		/// 
 		/// </summary>
 		static void CreateCoonsSurface(const LN_NurbsCurve& leftCurve, const LN_NurbsCurve& bottomCurve, const LN_NurbsCurve& rightCurve, const LN_NurbsCurve& topCurve, LN_NurbsSurface& surface);
 

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);
 }