Divide Close Curve(Polyline) By List Of Area
Hi, I wrote c# code to split a close curve (multiple lines) based on area list
My Component has several options:
divide 2 part to Target Area = 0
divide in Equal Area = 1
divide by list of Area = 2 Explanation that this code can divide a closed curve in the direction of the y axis from top to bottom based on target areas
*i share my Code * If anyone can optimize this code or have a better suggestion to make it faster, please help to share divide Area new -----.gh (37.2 KB)
Indeed there’s a thing or two that require some attention more. But although I have various C#'s for that matter … they work differently so I can’t cut and paste Methods.
I’ll examine the thing in depth and I’ll be back soon (but all my similar things are in the practice that reopens April 18).
In general one suggested ability (more pragmatinc for facades and the likes - if you are after that kind of stuff) more would be using BrepFaces. Like:
Hello
for me a simple way is to move the cutter from the beginning to the end with a specified number of steps. For each step calculate the area of the intersection.
Also all the calculation could be done with parallel calculations. So it take 0.2 s for 1000 calculations
private void RunScript(Curve cutter, Curve cutted, ref object dys, ref object areas)
{
BoundingBox bb_cutted = cutted.GetBoundingBox(true);
BoundingBox bb_cutter = cutter.GetBoundingBox(true);
double maxYtranslation = bb_cutter.Center.Y - bb_cutter.Diagonal.Y / 2 - (bb_cutted.Center.Y - bb_cutted.Diagonal.Y / 2 ) * 1.2;
int nSteps = 1000;
double tolerance = 0.01;
//Mouvement in -Y
double[] lst_dy = new double[nSteps];
//Area of each intersection
double[] lst_area = new double[nSteps];
Polyline pl_cutted;
ClipperLib.Polyline3D.ConvertCurveToPolyline(cutted, tolerance, out pl_cutted);
Parallel.For(0, nSteps, i =>
{
double dy = maxYtranslation * (double) i / (double) (nSteps - 1);
Vector3d vY = new Vector3d(0, -dy, 0);
lst_dy[i] = dy;
Curve cutterTranslated = cutter.DuplicateCurve();
cutterTranslated.Transform(Transform.Translation(vY));
Polyline pl_cutter;
ClipperLib.Polyline3D.ConvertCurveToPolyline(cutterTranslated, tolerance, out pl_cutter);
List<Polyline> pls = ClipperLib.Polyline3D.Boolean(ClipperLib.ClipType.ctIntersection, pl_cutter, pl_cutted, Plane.WorldXY, tolerance, true);
double area = 0;
foreach (Polyline pl in pls)
{
area += LDLIB.PrintTools.PolylineArea(pl, Plane.WorldXY);
}
lst_area[i] = area;
}
);
dys = lst_dy;
areas = lst_area;
}
}
it is a kind of brute force
you can use more optimized tool like this one that was used recently by @DanielPiker
the function f(x) is the area depending on the movement in Y.
Similar concept to what already said by others, Curve-Curve boolean intersection + area calc done inside a c# in parallel as Laurent said.
Pseudo-geometric solution, pseudo brute force.
Intersection with the “intergral graph”, which have kinks where the section curve intercepts a discontinuity/corner on the area curve to split.
10 minutes rework
Here a solution with only Rhinocommon and no plugin
The inputs must be closed curves. But it is possible to have a cutter not using a closed curve and close it in the code. It is also possible to the other things inside.
private void RunScript(Curve cutter, Curve cutted, double area, ref object vector)
{
//Laurent Delrieu 21/04/2023
double tolerance = 0.001;
sCC = new sliceClosedCurve(cutter, cutted, tolerance);
double dy = RootFinding.Brent(new FunctionOfOneVariable(f), 0, sCC.maxYtranslation, tolerance, area);
//Mouvement of cutter in order to have the area
vector = new Vector3d(0, -dy, 0);
}
// <Custom additional code>
static sliceClosedCurve sCC;
/// <summary>
/// Daniel Piker tool transformed
///https://discourse.mcneel.com/t/sweep-plane-until-lower-section-has-a-certain-volume/128708/36
/// </summary>
public class sliceClosedCurve
{
Curve cutter;
Curve cutted;
public double maxYtranslation;
double tolerance;
public sliceClosedCurve(Curve c_cutter, Curve c_cutted, double tol)
{
cutter = c_cutter.DuplicateCurve();
cutted = c_cutted.DuplicateCurve();
tolerance = tol;
BoundingBox bb_cutted = cutted.GetBoundingBox(true);
BoundingBox bb_cutter = cutter.GetBoundingBox(true);
maxYtranslation = bb_cutter.Center.Y - bb_cutter.Diagonal.Y / 2 - (bb_cutted.Center.Y - bb_cutted.Diagonal.Y / 2 ) * 1.2;
}
public double getArea(double dy)
{
Vector3d vY = new Vector3d(0, -dy, 0);
Curve cutterTranslated = cutter.DuplicateCurve();
cutterTranslated.Transform(Transform.Translation(vY));
Curve[] tab_curves = Curve.CreateBooleanIntersection(cutterTranslated, cutted, tolerance);
double area = 0;
foreach (Curve curve in tab_curves)
{
var s = Rhino.Geometry.AreaMassProperties.Compute(curve);
area += s.Area;
}
return area;
}
}
public static double f(double x)
{
return sCC.getArea(x);
}
/// John D. Cook
/// https://www.codeproject.com/Articles/79541/Three-Methods-for-Root-finding-in-C
/// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
public delegate double FunctionOfOneVariable(double x);
public static class RootFinding
{
const int maxIterations = 50;
/// <summary>
/// John D. Cook
/// https://www.codeproject.com/Articles/79541/Three-Methods-for-Root-finding-in-C
/// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
/// </summary>
/// <param name="f"></param>
/// <param name="left"></param>
/// <param name="right"></param>
/// <param name="tolerance"></param>
/// <param name="target"></param>
/// <returns></returns>
public static double Bisect
(
FunctionOfOneVariable f,
double left,
double right,
double tolerance = 1e-6,
double target = 0.0
)
{
// extra info that callers may not always want
int iterationsUsed;
double errorEstimate;
return Bisect(f, left, right, tolerance, target, out iterationsUsed, out errorEstimate);
}
public static double Bisect
(
FunctionOfOneVariable f,
double left,
double right,
double tolerance,
double target,
out int iterationsUsed,
out double errorEstimate
)
{
if (tolerance <= 0.0)
{
string msg = string.Format("Tolerance must be positive. Recieved {0}.", tolerance);
throw new ArgumentOutOfRangeException(msg);
}
iterationsUsed = 0;
errorEstimate = double.MaxValue;
// Standardize the problem. To solve f(x) = target,
// solve g(x) = 0 where g(x) = f(x) - target.
FunctionOfOneVariable g = delegate(double x) { return f(x) - target; };
double g_left = g(left); // evaluation of f at left end of interval
double g_right = g(right);
double mid;
double g_mid;
if (g_left * g_right >= 0.0)
{
string str = "Invalid starting bracket. Function must be above target on one end and below target on other end.";
string msg = string.Format("{0} Target: {1}. f(left) = {2}. f(right) = {3}", str, g_left + target, g_right + target);
throw new ArgumentException(msg);
}
double intervalWidth = right - left;
for
(
iterationsUsed = 0;
iterationsUsed < maxIterations && intervalWidth > tolerance;
iterationsUsed++
)
{
intervalWidth *= 0.5;
mid = left + intervalWidth;
if ((g_mid = g(mid)) == 0.0)
{
errorEstimate = 0.0;
return mid;
}
if (g_left * g_mid < 0.0) // g changes sign in (left, mid)
g_right = g(right = mid);
else // g changes sign in (mid, right)
g_left = g(left = mid);
}
errorEstimate = right - left;
return left;
}
public static double Brent
(
FunctionOfOneVariable f,
double left,
double right,
double tolerance = 1e-6,
double target = 0.0
)
{
// extra info that callers may not always want
int iterationsUsed;
double errorEstimate;
return Brent(f, left, right, tolerance, target, out iterationsUsed, out errorEstimate);
}
public static double Brent
(
FunctionOfOneVariable g,
double left,
double right,
double tolerance,
double target,
out int iterationsUsed,
out double errorEstimate
)
{
if (tolerance <= 0.0)
{
string msg = string.Format("Tolerance must be positive. Recieved {0}.", tolerance);
throw new ArgumentOutOfRangeException(msg);
}
errorEstimate = double.MaxValue;
// Standardize the problem. To solve g(x) = target,
// solve f(x) = 0 where f(x) = g(x) - target.
FunctionOfOneVariable f = delegate(double x) { return g(x) - target; };
// Implementation and notation based on Chapter 4 in
// "Algorithms for Minimization without Derivatives"
// by Richard Brent.
double c, d, e, fa, fb, fc, tol, m, p, q, r, s;
// set up aliases to match Brent's notation
double a = left; double b = right; double t = tolerance;
iterationsUsed = 0;
fa = f(a);
fb = f(b);
if (fa * fb > 0.0)
{
string str = "Invalid starting bracket. Function must be above target on one end and below target on other end.";
string msg = string.Format("{0} Target: {1}. f(left) = {2}. f(right) = {3}", str, target, fa + target, fb + target);
throw new ArgumentException(msg);
}
label_int:
c = a; fc = fa; d = e = b - a;
label_ext:
if (Math.Abs(fc) < Math.Abs(fb))
{
a = b; b = c; c = a;
fa = fb; fb = fc; fc = fa;
}
iterationsUsed++;
tol = 2.0 * t * Math.Abs(b) + t;
errorEstimate = m = 0.5 * (c - b);
if (Math.Abs(m) > tol && fb != 0.0) // exact comparison with 0 is OK here
{
// See if bisection is forced
if (Math.Abs(e) < tol || Math.Abs(fa) <= Math.Abs(fb))
{
d = e = m;
}
else
{
s = fb / fa;
if (a == c)
{
// linear interpolation
p = 2.0 * m * s; q = 1.0 - s;
}
else
{
// Inverse quadratic interpolation
q = fa / fc; r = fb / fc;
p = s * (2.0 * m * q * (q - r) - (b - a) * (r - 1.0));
q = (q - 1.0) * (r - 1.0) * (s - 1.0);
}
if (p > 0.0)
q = -q;
else
p = -p;
s = e; e = d;
if (2.0 * p < 3.0 * m * q - Math.Abs(tol * q) && p < Math.Abs(0.5 * s * q))
d = p / q;
else
d = e = m;
}
a = b; fa = fb;
if (Math.Abs(d) > tol)
b += d;
else if (m > 0.0)
b += tol;
else
b -= tol;
if (iterationsUsed == maxIterations)
return b;
fb = f(b);
if ((fb > 0.0 && fc > 0.0) || (fb <= 0.0 && fc <= 0.0))
goto label_int;
else
goto label_ext;
}
else
return b;
}
public static double Newton
(
FunctionOfOneVariable f,
FunctionOfOneVariable fprime,
double guess,
double tolerance = 1e-6,
double target = 0.0
)
{
// extra info that callers may not always want
int iterationsUsed;
double errorEstimate;
return Newton(f, fprime, guess, tolerance, target, out iterationsUsed, out errorEstimate);
}
public static double Newton
(
FunctionOfOneVariable f,
FunctionOfOneVariable fprime,
double guess,
double tolerance,
double target,
out int iterationsUsed,
out double errorEstimate
)
{
if (tolerance <= 0)
{
string msg = string.Format("Tolerance must be positive. Recieved {0}.", tolerance);
throw new ArgumentOutOfRangeException(msg);
}
iterationsUsed = 0;
errorEstimate = double.MaxValue;
// Standardize the problem. To solve f(x) = target,
// solve g(x) = 0 where g(x) = f(x) - target.
// Note that f(x) and g(x) have the same derivative.
FunctionOfOneVariable g = delegate(double x) { return f(x) - target; };
double oldX, newX = guess;
for
(
iterationsUsed = 0;
iterationsUsed < maxIterations && errorEstimate > tolerance;
iterationsUsed++
)
{
oldX = newX;
double gx = g(oldX);
double gprimex = fprime(oldX);
double absgprimex = Math.Abs(gprimex);
if (absgprimex > 1.0 || Math.Abs(gx) < double.MaxValue * absgprimex)
{
// The division will not overflow
newX = oldX - gx / gprimex;
errorEstimate = Math.Abs(newX - oldX);
}
else
{
newX = oldX;
errorEstimate = double.MaxValue;
break;
}
}
return newX;
}
}
I didn’t test this one. The only problem I seebis that you will have to sort the splitted region.
Combining the 2 methods is surely the way to go. I made 3 methods you have another one with Grasshopper component plus the one from you. Up to you to decide from the speed reliability simplicity or whatever.
For my tools I often use clipper for reliability. For example if I made intersection of 2 same curve it work better with clipper than RhinoCommon
