From 216d3f60dd2f1d790cb17aaac45b947c6d3c78c9 Mon Sep 17 00:00:00 2001 From: Laurent Le Goff Date: Wed, 22 Apr 2015 19:07:03 +0200 Subject: [PATCH] Clean curve package --- curve/arc.go | 4 +- curve/cubic_float64.go | 49 ++- curve/cubic_float64_others.go | 696 ---------------------------------- curve/curve_test.go | 156 +------- curve/quad_float64.go | 18 +- curve/tracer.go | 6 + 6 files changed, 65 insertions(+), 864 deletions(-) delete mode 100644 curve/cubic_float64_others.go create mode 100644 curve/tracer.go diff --git a/curve/arc.go b/curve/arc.go index 583c158..9c8fb82 100644 --- a/curve/arc.go +++ b/curve/arc.go @@ -1,12 +1,14 @@ // Copyright 2010 The draw2d Authors. All rights reserved. // created: 21/11/2010 by Laurent Le Goff + package curve import ( "math" ) -func SegmentArc(t LineTracer, x, y, rx, ry, start, angle, scale float64) { +// TraceArc trace an arc using a LineTracer +func TraceArc(t LineTracer, x, y, rx, ry, start, angle, scale float64) { end := start + angle clockWise := true if angle < 0 { diff --git a/curve/cubic_float64.go b/curve/cubic_float64.go index 64a7ac6..8c19eb2 100644 --- a/curve/cubic_float64.go +++ b/curve/cubic_float64.go @@ -1,5 +1,7 @@ // Copyright 2010 The draw2d Authors. All rights reserved. // created: 17/05/2011 by Laurent Le Goff + +// Package curve implements Bezier Curve Subdivision using De Casteljau's algorithm package curve import ( @@ -10,38 +12,48 @@ const ( CurveRecursionLimit = 32 ) -// X1, Y1, X2, Y2, X3, Y3, X4, Y4 float64 +// x1, y1, cpx1, cpx2, cpx2, cpy2, x2, y2 float64 type CubicCurveFloat64 [8]float64 -type LineTracer interface { - LineTo(x, y float64) -} - -func (c *CubicCurveFloat64) Subdivide(c1, c2 *CubicCurveFloat64) (x23, y23 float64) { - // Calculate all the mid-points of the line segments - //---------------------- +// Subdivide a Bezier cubic curve in 2 equivalents Bezier cubic curves. +// c1 and c2 parameters are the resulting curves +func (c *CubicCurveFloat64) Subdivide(c1, c2 *CubicCurveFloat64) { + // First point of c is the first point of c1 c1[0], c1[1] = c[0], c[1] + // Last point of c is the last point of c2 c2[6], c2[7] = c[6], c[7] + + // Subdivide segment using midpoints c1[2] = (c[0] + c[2]) / 2 c1[3] = (c[1] + c[3]) / 2 - x23 = (c[2] + c[4]) / 2 - y23 = (c[3] + c[5]) / 2 + + midX := (c[2] + c[4]) / 2 + midY := (c[3] + c[5]) / 2 + c2[4] = (c[4] + c[6]) / 2 c2[5] = (c[5] + c[7]) / 2 - c1[4] = (c1[2] + x23) / 2 - c1[5] = (c1[3] + y23) / 2 - c2[2] = (x23 + c2[4]) / 2 - c2[3] = (y23 + c2[5]) / 2 + + c1[4] = (c1[2] + midX) / 2 + c1[5] = (c1[3] + midY) / 2 + + c2[2] = (midX + c2[4]) / 2 + c2[3] = (midY + c2[5]) / 2 + c1[6] = (c1[4] + c2[2]) / 2 c1[7] = (c1[5] + c2[3]) / 2 + + // Last Point of c1 is equal to the first point of c2 c2[0], c2[1] = c1[6], c1[7] - return } -func (curve *CubicCurveFloat64) Segment(t LineTracer, flattening_threshold float64) { +// Trace generate lines subdividing the curve using a LineTracer +// flattening_threshold helps determines the flattening expectation of the curve +func (curve *CubicCurveFloat64) Trace(t LineTracer, flattening_threshold float64) { + // Allocation curves var curves [CurveRecursionLimit]CubicCurveFloat64 curves[0] = *curve i := 0 + // current curve var c *CubicCurveFloat64 @@ -52,9 +64,10 @@ func (curve *CubicCurveFloat64) Segment(t LineTracer, flattening_threshold float dx = c[6] - c[0] dy = c[7] - c[1] - d2 = math.Abs(((c[2]-c[6])*dy - (c[3]-c[7])*dx)) - d3 = math.Abs(((c[4]-c[6])*dy - (c[5]-c[7])*dx)) + d2 = math.Abs((c[2]-c[6])*dy - (c[3]-c[7])*dx) + d3 = math.Abs((c[4]-c[6])*dy - (c[5]-c[7])*dx) + // if it's flat then trace a line if (d2+d3)*(d2+d3) < flattening_threshold*(dx*dx+dy*dy) || i == len(curves)-1 { t.LineTo(c[6], c[7]) i-- diff --git a/curve/cubic_float64_others.go b/curve/cubic_float64_others.go deleted file mode 100644 index 497f1d3..0000000 --- a/curve/cubic_float64_others.go +++ /dev/null @@ -1,696 +0,0 @@ -// Copyright 2010 The draw2d Authors. All rights reserved. -// created: 17/05/2011 by Laurent Le Goff -package curve - -import ( - "math" -) - -const ( - CurveCollinearityEpsilon = 1e-30 - CurveAngleToleranceEpsilon = 0.01 -) - -//mu ranges from 0 to 1, start to end of curve -func (c *CubicCurveFloat64) ArbitraryPoint(mu float64) (x, y float64) { - - mum1 := 1 - mu - mum13 := mum1 * mum1 * mum1 - mu3 := mu * mu * mu - - x = mum13*c[0] + 3*mu*mum1*mum1*c[2] + 3*mu*mu*mum1*c[4] + mu3*c[6] - y = mum13*c[1] + 3*mu*mum1*mum1*c[3] + 3*mu*mu*mum1*c[5] + mu3*c[7] - return -} - -func (c *CubicCurveFloat64) SubdivideAt(c1, c2 *CubicCurveFloat64, t float64) (x23, y23 float64) { - inv_t := (1 - t) - c1[0], c1[1] = c[0], c[1] - c2[6], c2[7] = c[6], c[7] - - c1[2] = inv_t*c[0] + t*c[2] - c1[3] = inv_t*c[1] + t*c[3] - - x23 = inv_t*c[2] + t*c[4] - y23 = inv_t*c[3] + t*c[5] - - c2[4] = inv_t*c[4] + t*c[6] - c2[5] = inv_t*c[5] + t*c[7] - - c1[4] = inv_t*c1[2] + t*x23 - c1[5] = inv_t*c1[3] + t*y23 - - c2[2] = inv_t*x23 + t*c2[4] - c2[3] = inv_t*y23 + t*c2[5] - - c1[6] = inv_t*c1[4] + t*c2[2] - c1[7] = inv_t*c1[5] + t*c2[3] - - c2[0], c2[1] = c1[6], c1[7] - return -} - -func (c *CubicCurveFloat64) EstimateDistance() float64 { - dx1 := c[2] - c[0] - dy1 := c[3] - c[1] - dx2 := c[4] - c[2] - dy2 := c[5] - c[3] - dx3 := c[6] - c[4] - dy3 := c[7] - c[5] - return math.Sqrt(dx1*dx1+dy1*dy1) + math.Sqrt(dx2*dx2+dy2*dy2) + math.Sqrt(dx3*dx3+dy3*dy3) -} - -// subdivide the curve in straight lines using line approximation and Casteljau recursive subdivision -func (c *CubicCurveFloat64) SegmentRec(t LineTracer, flattening_threshold float64) { - c.segmentRec(t, flattening_threshold) - t.LineTo(c[6], c[7]) -} - -func (c *CubicCurveFloat64) segmentRec(t LineTracer, flattening_threshold float64) { - var c1, c2 CubicCurveFloat64 - c.Subdivide(&c1, &c2) - - // Try to approximate the full cubic curve by a single straight line - //------------------ - dx := c[6] - c[0] - dy := c[7] - c[1] - - d2 := math.Abs(((c[2]-c[6])*dy - (c[3]-c[7])*dx)) - d3 := math.Abs(((c[4]-c[6])*dy - (c[5]-c[7])*dx)) - - if (d2+d3)*(d2+d3) < flattening_threshold*(dx*dx+dy*dy) { - t.LineTo(c[6], c[7]) - return - } - // Continue subdivision - //---------------------- - c1.segmentRec(t, flattening_threshold) - c2.segmentRec(t, flattening_threshold) -} - -/* - The function has the following parameters: - approximationScale : - Eventually determines the approximation accuracy. In practice we need to transform points from the World coordinate system to the Screen one. - It always has some scaling coefficient. - The curves are usually processed in the World coordinates, while the approximation accuracy should be eventually in pixels. - Usually it looks as follows: - curved.approximationScale(transform.scale()); - where transform is the affine matrix that includes all the transformations, including viewport and zoom. - angleTolerance : - You set it in radians. - The less this value is the more accurate will be the approximation at sharp turns. - But 0 means that we don't consider angle conditions at all. - cuspLimit : - An angle in radians. - If 0, only the real cusps will have bevel cuts. - If more than 0, it will restrict the sharpness. - The more this value is the less sharp turns will be cut. - Typically it should not exceed 10-15 degrees. -*/ -func (c *CubicCurveFloat64) AdaptiveSegmentRec(t LineTracer, approximationScale, angleTolerance, cuspLimit float64) { - cuspLimit = computeCuspLimit(cuspLimit) - distanceToleranceSquare := 0.5 / approximationScale - distanceToleranceSquare = distanceToleranceSquare * distanceToleranceSquare - c.adaptiveSegmentRec(t, 0, distanceToleranceSquare, angleTolerance, cuspLimit) - t.LineTo(c[6], c[7]) -} - -func computeCuspLimit(v float64) (r float64) { - if v == 0.0 { - r = 0.0 - } else { - r = math.Pi - v - } - return -} - -func squareDistance(x1, y1, x2, y2 float64) float64 { - dx := x2 - x1 - dy := y2 - y1 - return dx*dx + dy*dy -} - -/** - * http://www.antigrain.com/research/adaptive_bezier/index.html - */ -func (c *CubicCurveFloat64) adaptiveSegmentRec(t LineTracer, level int, distanceToleranceSquare, angleTolerance, cuspLimit float64) { - if level > CurveRecursionLimit { - return - } - var c1, c2 CubicCurveFloat64 - x23, y23 := c.Subdivide(&c1, &c2) - - // Try to approximate the full cubic curve by a single straight line - //------------------ - dx := c[6] - c[0] - dy := c[7] - c[1] - - d2 := math.Abs(((c[2]-c[6])*dy - (c[3]-c[7])*dx)) - d3 := math.Abs(((c[4]-c[6])*dy - (c[5]-c[7])*dx)) - switch { - case d2 <= CurveCollinearityEpsilon && d3 <= CurveCollinearityEpsilon: - // All collinear OR p1==p4 - //---------------------- - k := dx*dx + dy*dy - if k == 0 { - d2 = squareDistance(c[0], c[1], c[2], c[3]) - d3 = squareDistance(c[6], c[7], c[4], c[5]) - } else { - k = 1 / k - da1 := c[2] - c[0] - da2 := c[3] - c[1] - d2 = k * (da1*dx + da2*dy) - da1 = c[4] - c[0] - da2 = c[5] - c[1] - d3 = k * (da1*dx + da2*dy) - if d2 > 0 && d2 < 1 && d3 > 0 && d3 < 1 { - // Simple collinear case, 1---2---3---4 - // We can leave just two endpoints - return - } - if d2 <= 0 { - d2 = squareDistance(c[2], c[3], c[0], c[1]) - } else if d2 >= 1 { - d2 = squareDistance(c[2], c[3], c[6], c[7]) - } else { - d2 = squareDistance(c[2], c[3], c[0]+d2*dx, c[1]+d2*dy) - } - - if d3 <= 0 { - d3 = squareDistance(c[4], c[5], c[0], c[1]) - } else if d3 >= 1 { - d3 = squareDistance(c[4], c[5], c[6], c[7]) - } else { - d3 = squareDistance(c[4], c[5], c[0]+d3*dx, c[1]+d3*dy) - } - } - if d2 > d3 { - if d2 < distanceToleranceSquare { - t.LineTo(c[2], c[3]) - return - } - } else { - if d3 < distanceToleranceSquare { - t.LineTo(c[4], c[5]) - return - } - } - - case d2 <= CurveCollinearityEpsilon && d3 > CurveCollinearityEpsilon: - // p1,p2,p4 are collinear, p3 is significant - //---------------------- - if d3*d3 <= distanceToleranceSquare*(dx*dx+dy*dy) { - if angleTolerance < CurveAngleToleranceEpsilon { - t.LineTo(x23, y23) - return - } - - // Angle Condition - //---------------------- - da1 := math.Abs(math.Atan2(c[7]-c[5], c[6]-c[4]) - math.Atan2(c[5]-c[3], c[4]-c[2])) - if da1 >= math.Pi { - da1 = 2*math.Pi - da1 - } - - if da1 < angleTolerance { - t.LineTo(c[2], c[3]) - t.LineTo(c[4], c[5]) - return - } - - if cuspLimit != 0.0 { - if da1 > cuspLimit { - t.LineTo(c[4], c[5]) - return - } - } - } - - case d2 > CurveCollinearityEpsilon && d3 <= CurveCollinearityEpsilon: - // p1,p3,p4 are collinear, p2 is significant - //---------------------- - if d2*d2 <= distanceToleranceSquare*(dx*dx+dy*dy) { - if angleTolerance < CurveAngleToleranceEpsilon { - t.LineTo(x23, y23) - return - } - - // Angle Condition - //---------------------- - da1 := math.Abs(math.Atan2(c[5]-c[3], c[4]-c[2]) - math.Atan2(c[3]-c[1], c[2]-c[0])) - if da1 >= math.Pi { - da1 = 2*math.Pi - da1 - } - - if da1 < angleTolerance { - t.LineTo(c[2], c[3]) - t.LineTo(c[4], c[5]) - return - } - - if cuspLimit != 0.0 { - if da1 > cuspLimit { - t.LineTo(c[2], c[3]) - return - } - } - } - - case d2 > CurveCollinearityEpsilon && d3 > CurveCollinearityEpsilon: - // Regular case - //----------------- - if (d2+d3)*(d2+d3) <= distanceToleranceSquare*(dx*dx+dy*dy) { - // If the curvature doesn't exceed the distanceTolerance value - // we tend to finish subdivisions. - //---------------------- - if angleTolerance < CurveAngleToleranceEpsilon { - t.LineTo(x23, y23) - return - } - - // Angle & Cusp Condition - //---------------------- - k := math.Atan2(c[5]-c[3], c[4]-c[2]) - da1 := math.Abs(k - math.Atan2(c[3]-c[1], c[2]-c[0])) - da2 := math.Abs(math.Atan2(c[7]-c[5], c[6]-c[4]) - k) - if da1 >= math.Pi { - da1 = 2*math.Pi - da1 - } - if da2 >= math.Pi { - da2 = 2*math.Pi - da2 - } - - if da1+da2 < angleTolerance { - // Finally we can stop the recursion - //---------------------- - t.LineTo(x23, y23) - return - } - - if cuspLimit != 0.0 { - if da1 > cuspLimit { - t.LineTo(c[2], c[3]) - return - } - - if da2 > cuspLimit { - t.LineTo(c[4], c[5]) - return - } - } - } - } - - // Continue subdivision - //---------------------- - c1.adaptiveSegmentRec(t, level+1, distanceToleranceSquare, angleTolerance, cuspLimit) - c2.adaptiveSegmentRec(t, level+1, distanceToleranceSquare, angleTolerance, cuspLimit) - -} - -func (curve *CubicCurveFloat64) AdaptiveSegment(t LineTracer, approximationScale, angleTolerance, cuspLimit float64) { - cuspLimit = computeCuspLimit(cuspLimit) - distanceToleranceSquare := 0.5 / approximationScale - distanceToleranceSquare = distanceToleranceSquare * distanceToleranceSquare - - var curves [CurveRecursionLimit]CubicCurveFloat64 - curves[0] = *curve - i := 0 - // current curve - var c *CubicCurveFloat64 - var c1, c2 CubicCurveFloat64 - var dx, dy, d2, d3, k, x23, y23 float64 - for i >= 0 { - c = &curves[i] - x23, y23 = c.Subdivide(&c1, &c2) - - // Try to approximate the full cubic curve by a single straight line - //------------------ - dx = c[6] - c[0] - dy = c[7] - c[1] - - d2 = math.Abs(((c[2]-c[6])*dy - (c[3]-c[7])*dx)) - d3 = math.Abs(((c[4]-c[6])*dy - (c[5]-c[7])*dx)) - switch { - case i == len(curves)-1: - t.LineTo(c[6], c[7]) - i-- - continue - case d2 <= CurveCollinearityEpsilon && d3 <= CurveCollinearityEpsilon: - // All collinear OR p1==p4 - //---------------------- - k = dx*dx + dy*dy - if k == 0 { - d2 = squareDistance(c[0], c[1], c[2], c[3]) - d3 = squareDistance(c[6], c[7], c[4], c[5]) - } else { - k = 1 / k - da1 := c[2] - c[0] - da2 := c[3] - c[1] - d2 = k * (da1*dx + da2*dy) - da1 = c[4] - c[0] - da2 = c[5] - c[1] - d3 = k * (da1*dx + da2*dy) - if d2 > 0 && d2 < 1 && d3 > 0 && d3 < 1 { - // Simple collinear case, 1---2---3---4 - // We can leave just two endpoints - i-- - continue - } - if d2 <= 0 { - d2 = squareDistance(c[2], c[3], c[0], c[1]) - } else if d2 >= 1 { - d2 = squareDistance(c[2], c[3], c[6], c[7]) - } else { - d2 = squareDistance(c[2], c[3], c[0]+d2*dx, c[1]+d2*dy) - } - - if d3 <= 0 { - d3 = squareDistance(c[4], c[5], c[0], c[1]) - } else if d3 >= 1 { - d3 = squareDistance(c[4], c[5], c[6], c[7]) - } else { - d3 = squareDistance(c[4], c[5], c[0]+d3*dx, c[1]+d3*dy) - } - } - if d2 > d3 { - if d2 < distanceToleranceSquare { - t.LineTo(c[2], c[3]) - i-- - continue - } - } else { - if d3 < distanceToleranceSquare { - t.LineTo(c[4], c[5]) - i-- - continue - } - } - - case d2 <= CurveCollinearityEpsilon && d3 > CurveCollinearityEpsilon: - // p1,p2,p4 are collinear, p3 is significant - //---------------------- - if d3*d3 <= distanceToleranceSquare*(dx*dx+dy*dy) { - if angleTolerance < CurveAngleToleranceEpsilon { - t.LineTo(x23, y23) - i-- - continue - } - - // Angle Condition - //---------------------- - da1 := math.Abs(math.Atan2(c[7]-c[5], c[6]-c[4]) - math.Atan2(c[5]-c[3], c[4]-c[2])) - if da1 >= math.Pi { - da1 = 2*math.Pi - da1 - } - - if da1 < angleTolerance { - t.LineTo(c[2], c[3]) - t.LineTo(c[4], c[5]) - i-- - continue - } - - if cuspLimit != 0.0 { - if da1 > cuspLimit { - t.LineTo(c[4], c[5]) - i-- - continue - } - } - } - - case d2 > CurveCollinearityEpsilon && d3 <= CurveCollinearityEpsilon: - // p1,p3,p4 are collinear, p2 is significant - //---------------------- - if d2*d2 <= distanceToleranceSquare*(dx*dx+dy*dy) { - if angleTolerance < CurveAngleToleranceEpsilon { - t.LineTo(x23, y23) - i-- - continue - } - - // Angle Condition - //---------------------- - da1 := math.Abs(math.Atan2(c[5]-c[3], c[4]-c[2]) - math.Atan2(c[3]-c[1], c[2]-c[0])) - if da1 >= math.Pi { - da1 = 2*math.Pi - da1 - } - - if da1 < angleTolerance { - t.LineTo(c[2], c[3]) - t.LineTo(c[4], c[5]) - i-- - continue - } - - if cuspLimit != 0.0 { - if da1 > cuspLimit { - t.LineTo(c[2], c[3]) - i-- - continue - } - } - } - - case d2 > CurveCollinearityEpsilon && d3 > CurveCollinearityEpsilon: - // Regular case - //----------------- - if (d2+d3)*(d2+d3) <= distanceToleranceSquare*(dx*dx+dy*dy) { - // If the curvature doesn't exceed the distanceTolerance value - // we tend to finish subdivisions. - //---------------------- - if angleTolerance < CurveAngleToleranceEpsilon { - t.LineTo(x23, y23) - i-- - continue - } - - // Angle & Cusp Condition - //---------------------- - k := math.Atan2(c[5]-c[3], c[4]-c[2]) - da1 := math.Abs(k - math.Atan2(c[3]-c[1], c[2]-c[0])) - da2 := math.Abs(math.Atan2(c[7]-c[5], c[6]-c[4]) - k) - if da1 >= math.Pi { - da1 = 2*math.Pi - da1 - } - if da2 >= math.Pi { - da2 = 2*math.Pi - da2 - } - - if da1+da2 < angleTolerance { - // Finally we can stop the recursion - //---------------------- - t.LineTo(x23, y23) - i-- - continue - } - - if cuspLimit != 0.0 { - if da1 > cuspLimit { - t.LineTo(c[2], c[3]) - i-- - continue - } - - if da2 > cuspLimit { - t.LineTo(c[4], c[5]) - i-- - continue - } - } - } - } - - // Continue subdivision - //---------------------- - curves[i+1], curves[i] = c1, c2 - i++ - } - t.LineTo(curve[6], curve[7]) -} - -/********************** Ahmad thesis *******************/ - -/************************************************************************************** -* This code is the implementation of the Parabolic Approximation (PA). Although * -* it uses recursive subdivision as a safe net for the failing cases, this is an * -* iterative routine and reduces considerably the number of vertices (point) * -* generation. * -**************************************************************************************/ - -func (c *CubicCurveFloat64) ParabolicSegment(t LineTracer, flattening_threshold float64) { - estimatedIFP := c.numberOfInflectionPoints() - if estimatedIFP == 0 { - // If no inflection points then apply PA on the full Bezier segment. - c.doParabolicApproximation(t, flattening_threshold) - return - } - // If one or more inflection point then we will have to subdivide the curve - numOfIfP, t1, t2 := c.findInflectionPoints() - if numOfIfP == 2 { - // Case when 2 inflection points then divide at the smallest one first - var sub1, tmp1, sub2, sub3 CubicCurveFloat64 - c.SubdivideAt(&sub1, &tmp1, t1) - // Now find the second inflection point in the second curve an subdivide - numOfIfP, t1, t2 = tmp1.findInflectionPoints() - if numOfIfP == 2 { - tmp1.SubdivideAt(&sub2, &sub3, t2) - } else if numOfIfP == 1 { - tmp1.SubdivideAt(&sub2, &sub3, t1) - } else { - return - } - // Use PA for first subsegment - sub1.doParabolicApproximation(t, flattening_threshold) - // Use RS for the second (middle) subsegment - sub2.Segment(t, flattening_threshold) - // Drop the last point in the array will be added by the PA in third subsegment - //noOfPoints--; - // Use PA for the third curve - sub3.doParabolicApproximation(t, flattening_threshold) - } else if numOfIfP == 1 { - // Case where there is one inflection point, subdivide once and use PA on - // both subsegments - var sub1, sub2 CubicCurveFloat64 - c.SubdivideAt(&sub1, &sub2, t1) - sub1.doParabolicApproximation(t, flattening_threshold) - //noOfPoints--; - sub2.doParabolicApproximation(t, flattening_threshold) - } else { - // Case where there is no inflection USA PA directly - c.doParabolicApproximation(t, flattening_threshold) - } -} - -// Find the third control point deviation form the axis -func (c *CubicCurveFloat64) thirdControlPointDeviation() float64 { - dx := c[2] - c[0] - dy := c[3] - c[1] - l2 := dx*dx + dy*dy - if l2 == 0 { - return 0 - } - l := math.Sqrt(l2) - r := (c[3] - c[1]) / l - s := (c[0] - c[2]) / l - u := (c[2]*c[1] - c[0]*c[3]) / l - return math.Abs(r*c[4] + s*c[5] + u) -} - -// Find the number of inflection point -func (c *CubicCurveFloat64) numberOfInflectionPoints() int { - dx21 := (c[2] - c[0]) - dy21 := (c[3] - c[1]) - dx32 := (c[4] - c[2]) - dy32 := (c[5] - c[3]) - dx43 := (c[6] - c[4]) - dy43 := (c[7] - c[5]) - if ((dx21*dy32 - dy21*dx32) * (dx32*dy43 - dy32*dx43)) < 0 { - return 1 // One inflection point - } else if ((dx21*dy32 - dy21*dx32) * (dx21*dy43 - dy21*dx43)) > 0 { - return 0 // No inflection point - } else { - // Most cases no inflection point - b1 := (dx21*dx32 + dy21*dy32) > 0 - b2 := (dx32*dx43 + dy32*dy43) > 0 - if b1 || b2 && !(b1 && b2) { // xor!! - return 0 - } - } - return -1 // cases where there in zero or two inflection points -} - -// This is the main function where all the work is done -func (curve *CubicCurveFloat64) doParabolicApproximation(tracer LineTracer, flattening_threshold float64) { - var c *CubicCurveFloat64 - c = curve - var d, t, dx, dy, d2, d3 float64 - for { - dx = c[6] - c[0] - dy = c[7] - c[1] - - d2 = math.Abs(((c[2]-c[6])*dy - (c[3]-c[7])*dx)) - d3 = math.Abs(((c[4]-c[6])*dy - (c[5]-c[7])*dx)) - - if (d2+d3)*(d2+d3) < flattening_threshold*(dx*dx+dy*dy) { - // If the subsegment deviation satisfy the flatness then store the last - // point and stop - tracer.LineTo(c[6], c[7]) - break - } - // Find the third control point deviation and the t values for subdivision - d = c.thirdControlPointDeviation() - t = 2 * math.Sqrt(flattening_threshold/d/3) - if t > 1 { - // Case where the t value calculated is invalid so using RS - c.Segment(tracer, flattening_threshold) - break - } - // Valid t value to subdivide at that calculated value - var b1, b2 CubicCurveFloat64 - c.SubdivideAt(&b1, &b2, t) - // First subsegment should have its deviation equal to flatness - dx = b1[6] - b1[0] - dy = b1[7] - b1[1] - - d2 = math.Abs(((b1[2]-b1[6])*dy - (b1[3]-b1[7])*dx)) - d3 = math.Abs(((b1[4]-b1[6])*dy - (b1[5]-b1[7])*dx)) - - if (d2+d3)*(d2+d3) > flattening_threshold*(dx*dx+dy*dy) { - // if not then use RS to handle any mathematical errors - b1.Segment(tracer, flattening_threshold) - } else { - tracer.LineTo(b1[6], b1[7]) - } - // repeat the process for the left over subsegment. - c = &b2 - } -} - -// Find the actual inflection points and return the number of inflection points found -// if 2 inflection points found, the first one returned will be with smaller t value. -func (curve *CubicCurveFloat64) findInflectionPoints() (int, firstIfp, secondIfp float64) { - // For Cubic Bezier curve with equation P=a*t^3 + b*t^2 + c*t + d - // slope of the curve dP/dt = 3*a*t^2 + 2*b*t + c - // a = (float)(-bez.p1 + 3*bez.p2 - 3*bez.p3 + bez.p4); - // b = (float)(3*bez.p1 - 6*bez.p2 + 3*bez.p3); - // c = (float)(-3*bez.p1 + 3*bez.p2); - ax := (-curve[0] + 3*curve[2] - 3*curve[4] + curve[6]) - bx := (3*curve[0] - 6*curve[2] + 3*curve[4]) - cx := (-3*curve[0] + 3*curve[2]) - ay := (-curve[1] + 3*curve[3] - 3*curve[5] + curve[7]) - by := (3*curve[1] - 6*curve[3] + 3*curve[5]) - cy := (-3*curve[1] + 3*curve[3]) - a := (3 * (ay*bx - ax*by)) - b := (3 * (ay*cx - ax*cy)) - c := (by*cx - bx*cy) - r2 := (b*b - 4*a*c) - firstIfp = 0.0 - secondIfp = 0.0 - if r2 >= 0.0 && a != 0.0 { - r := math.Sqrt(r2) - firstIfp = ((-b + r) / (2 * a)) - secondIfp = ((-b - r) / (2 * a)) - if (firstIfp > 0.0 && firstIfp < 1.0) && (secondIfp > 0.0 && secondIfp < 1.0) { - if firstIfp > secondIfp { - tmp := firstIfp - firstIfp = secondIfp - secondIfp = tmp - } - if secondIfp-firstIfp > 0.00001 { - return 2, firstIfp, secondIfp - } else { - return 1, firstIfp, secondIfp - } - } else if firstIfp > 0.0 && firstIfp < 1.0 { - return 1, firstIfp, secondIfp - } else if secondIfp > 0.0 && secondIfp < 1.0 { - firstIfp = secondIfp - return 1, firstIfp, secondIfp - } - return 0, firstIfp, secondIfp - } - return 0, firstIfp, secondIfp -} diff --git a/curve/curve_test.go b/curve/curve_test.go index c101e6a..2b4bb6a 100644 --- a/curve/curve_test.go +++ b/curve/curve_test.go @@ -1,16 +1,15 @@ package curve import ( - "bufio" "fmt" "image" "image/color" "image/draw" - "image/png" "log" "os" "testing" + "github.com/llgcode/draw2d" "github.com/llgcode/draw2d/raster" ) @@ -51,7 +50,8 @@ func (p *Path) LineTo(x, y float64) { } func init() { - f, err := os.Create("_test.html") + os.Mkdir("test_results", 0666) + f, err := os.Create("test_results/_test.html") if err != nil { log.Println(err) os.Exit(1) @@ -60,7 +60,7 @@ func init() { log.Printf("Create html viewer") f.Write([]byte("")) for i := 0; i < len(testsCubicFloat64); i++ { - f.Write([]byte(fmt.Sprintf("
\n\n\n\n\n
\n", i, i, i, i, i))) + f.Write([]byte(fmt.Sprintf("
\n", i))) } for i := 0; i < len(testsQuadFloat64); i++ { f.Write([]byte(fmt.Sprintf("
\n
\n", i))) @@ -69,28 +69,8 @@ func init() { } -func savepng(filePath string, m image.Image) { - f, err := os.Create(filePath) - if err != nil { - log.Println(err) - os.Exit(1) - } - defer f.Close() - b := bufio.NewWriter(f) - err = png.Encode(b, m) - if err != nil { - log.Println(err) - os.Exit(1) - } - err = b.Flush() - if err != nil { - log.Println(err) - os.Exit(1) - } -} - func drawPoints(img draw.Image, c color.Color, s ...float64) image.Image { - /*for i := 0; i < len(s); i += 2 { + for i := 0; i < len(s); i += 2 { x, y := int(s[i]+0.5), int(s[i+1]+0.5) img.Set(x, y, c) img.Set(x, y+1, c) @@ -102,85 +82,21 @@ func drawPoints(img draw.Image, c color.Color, s ...float64) image.Image { img.Set(x-1, y+1, c) img.Set(x-1, y-1, c) - }*/ - return img -} - -func TestCubicCurveRec(t *testing.T) { - for i, curve := range testsCubicFloat64 { - var p Path - p.LineTo(curve[0], curve[1]) - curve.SegmentRec(&p, flattening_threshold) - img := image.NewNRGBA(image.Rect(0, 0, 300, 300)) - raster.PolylineBresenham(img, color.NRGBA{0xff, 0, 0, 0xff}, curve[:]...) - raster.PolylineBresenham(img, image.Black, p.points...) - //drawPoints(img, image.NRGBAColor{0, 0, 0, 0xff}, curve[:]...) - drawPoints(img, color.NRGBA{0, 0, 0, 0xff}, p.points...) - savepng(fmt.Sprintf("_testRec%d.png", i), img) - log.Printf("Num of points: %d\n", len(p.points)) } - fmt.Println() + return img } func TestCubicCurve(t *testing.T) { for i, curve := range testsCubicFloat64 { var p Path p.LineTo(curve[0], curve[1]) - curve.Segment(&p, flattening_threshold) + curve.Trace(&p, flattening_threshold) img := image.NewNRGBA(image.Rect(0, 0, 300, 300)) raster.PolylineBresenham(img, color.NRGBA{0xff, 0, 0, 0xff}, curve[:]...) raster.PolylineBresenham(img, image.Black, p.points...) //drawPoints(img, image.NRGBAColor{0, 0, 0, 0xff}, curve[:]...) drawPoints(img, color.NRGBA{0, 0, 0, 0xff}, p.points...) - savepng(fmt.Sprintf("_test%d.png", i), img) - log.Printf("Num of points: %d\n", len(p.points)) - } - fmt.Println() -} - -func TestCubicCurveAdaptiveRec(t *testing.T) { - for i, curve := range testsCubicFloat64 { - var p Path - p.LineTo(curve[0], curve[1]) - curve.AdaptiveSegmentRec(&p, 1, 0, 0) - img := image.NewNRGBA(image.Rect(0, 0, 300, 300)) - raster.PolylineBresenham(img, color.NRGBA{0xff, 0, 0, 0xff}, curve[:]...) - raster.PolylineBresenham(img, image.Black, p.points...) - //drawPoints(img, image.NRGBAColor{0, 0, 0, 0xff}, curve[:]...) - drawPoints(img, color.NRGBA{0, 0, 0, 0xff}, p.points...) - savepng(fmt.Sprintf("_testAdaptiveRec%d.png", i), img) - log.Printf("Num of points: %d\n", len(p.points)) - } - fmt.Println() -} - -func TestCubicCurveAdaptive(t *testing.T) { - for i, curve := range testsCubicFloat64 { - var p Path - p.LineTo(curve[0], curve[1]) - curve.AdaptiveSegment(&p, 1, 0, 0) - img := image.NewNRGBA(image.Rect(0, 0, 300, 300)) - raster.PolylineBresenham(img, color.NRGBA{0xff, 0, 0, 0xff}, curve[:]...) - raster.PolylineBresenham(img, image.Black, p.points...) - //drawPoints(img, image.NRGBAColor{0, 0, 0, 0xff}, curve[:]...) - drawPoints(img, color.NRGBA{0, 0, 0, 0xff}, p.points...) - savepng(fmt.Sprintf("_testAdaptive%d.png", i), img) - log.Printf("Num of points: %d\n", len(p.points)) - } - fmt.Println() -} - -func TestCubicCurveParabolic(t *testing.T) { - for i, curve := range testsCubicFloat64 { - var p Path - p.LineTo(curve[0], curve[1]) - curve.ParabolicSegment(&p, flattening_threshold) - img := image.NewNRGBA(image.Rect(0, 0, 300, 300)) - raster.PolylineBresenham(img, color.NRGBA{0xff, 0, 0, 0xff}, curve[:]...) - raster.PolylineBresenham(img, image.Black, p.points...) - //drawPoints(img, image.NRGBAColor{0, 0, 0, 0xff}, curve[:]...) - drawPoints(img, color.NRGBA{0, 0, 0, 0xff}, p.points...) - savepng(fmt.Sprintf("_testParabolic%d.png", i), img) + draw2d.SaveToPngFile(fmt.Sprintf("test_results/_test%d.png", i), img) log.Printf("Num of points: %d\n", len(p.points)) } fmt.Println() @@ -190,74 +106,24 @@ func TestQuadCurve(t *testing.T) { for i, curve := range testsQuadFloat64 { var p Path p.LineTo(curve[0], curve[1]) - curve.Segment(&p, flattening_threshold) + curve.Trace(&p, flattening_threshold) img := image.NewNRGBA(image.Rect(0, 0, 300, 300)) raster.PolylineBresenham(img, color.NRGBA{0xff, 0, 0, 0xff}, curve[:]...) raster.PolylineBresenham(img, image.Black, p.points...) //drawPoints(img, image.NRGBAColor{0, 0, 0, 0xff}, curve[:]...) drawPoints(img, color.NRGBA{0, 0, 0, 0xff}, p.points...) - savepng(fmt.Sprintf("_testQuad%d.png", i), img) + draw2d.SaveToPngFile(fmt.Sprintf("test_results/_testQuad%d.png", i), img) log.Printf("Num of points: %d\n", len(p.points)) } fmt.Println() } -func BenchmarkCubicCurveRec(b *testing.B) { - for i := 0; i < b.N; i++ { - for _, curve := range testsCubicFloat64 { - p := Path{make([]float64, 0, 32)} - p.LineTo(curve[0], curve[1]) - curve.SegmentRec(&p, flattening_threshold) - } - } -} - func BenchmarkCubicCurve(b *testing.B) { for i := 0; i < b.N; i++ { for _, curve := range testsCubicFloat64 { p := Path{make([]float64, 0, 32)} p.LineTo(curve[0], curve[1]) - curve.Segment(&p, flattening_threshold) - } - } -} - -func BenchmarkCubicCurveAdaptiveRec(b *testing.B) { - for i := 0; i < b.N; i++ { - for _, curve := range testsCubicFloat64 { - p := Path{make([]float64, 0, 32)} - p.LineTo(curve[0], curve[1]) - curve.AdaptiveSegmentRec(&p, 1, 0, 0) - } - } -} - -func BenchmarkCubicCurveAdaptive(b *testing.B) { - for i := 0; i < b.N; i++ { - for _, curve := range testsCubicFloat64 { - p := Path{make([]float64, 0, 32)} - p.LineTo(curve[0], curve[1]) - curve.AdaptiveSegment(&p, 1, 0, 0) - } - } -} - -func BenchmarkCubicCurveParabolic(b *testing.B) { - for i := 0; i < b.N; i++ { - for _, curve := range testsCubicFloat64 { - p := Path{make([]float64, 0, 32)} - p.LineTo(curve[0], curve[1]) - curve.ParabolicSegment(&p, flattening_threshold) - } - } -} - -func BenchmarkQuadCurve(b *testing.B) { - for i := 0; i < b.N; i++ { - for _, curve := range testsQuadFloat64 { - p := Path{make([]float64, 0, 32)} - p.LineTo(curve[0], curve[1]) - curve.Segment(&p, flattening_threshold) + curve.Trace(&p, flattening_threshold) } } } diff --git a/curve/quad_float64.go b/curve/quad_float64.go index bd72aff..6cb2891 100644 --- a/curve/quad_float64.go +++ b/curve/quad_float64.go @@ -1,19 +1,25 @@ // Copyright 2010 The draw2d Authors. All rights reserved. // created: 17/05/2011 by Laurent Le Goff + package curve import ( "math" ) -//X1, Y1, X2, Y2, X3, Y3 float64 +//x1, y1, cpx1, cpy2, x2, y2 float64 type QuadCurveFloat64 [6]float64 +// Subdivide a Bezier quad curve in 2 equivalents Bezier quad curves. +// c1 and c2 parameters are the resulting curves func (c *QuadCurveFloat64) Subdivide(c1, c2 *QuadCurveFloat64) { - // Calculate all the mid-points of the line segments - //---------------------- + + // First point of c is the first point of c1 c1[0], c1[1] = c[0], c[1] + // Last point of c is the last point of c2 c2[4], c2[5] = c[4], c[5] + + // Subdivide segment using midpoints c1[2] = (c[0] + c[2]) / 2 c1[3] = (c[1] + c[3]) / 2 c2[2] = (c[2] + c[4]) / 2 @@ -24,7 +30,10 @@ func (c *QuadCurveFloat64) Subdivide(c1, c2 *QuadCurveFloat64) { return } -func (curve *QuadCurveFloat64) Segment(t LineTracer, flattening_threshold float64) { +// Trace generate lines subdividing the curve using a LineTracer +// flattening_threshold helps determines the flattening expectation of the curve +func (curve *QuadCurveFloat64) Trace(t LineTracer, flattening_threshold float64) { + // Allocates curves stack var curves [CurveRecursionLimit]QuadCurveFloat64 curves[0] = *curve i := 0 @@ -39,6 +48,7 @@ func (curve *QuadCurveFloat64) Segment(t LineTracer, flattening_threshold float6 d = math.Abs(((c[2]-c[4])*dy - (c[3]-c[5])*dx)) + // if it's flat then trace a line if (d*d) < flattening_threshold*(dx*dx+dy*dy) || i == len(curves)-1 { t.LineTo(c[4], c[5]) i-- diff --git a/curve/tracer.go b/curve/tracer.go new file mode 100644 index 0000000..e80ce1d --- /dev/null +++ b/curve/tracer.go @@ -0,0 +1,6 @@ +package curve + +// LineTracer is an interface that help segmenting curve into small lines +type LineTracer interface { + LineTo(x, y float64) +}