From b88f2dc3a35d7c2d026538bd55146b25a42d89f7 Mon Sep 17 00:00:00 2001 From: Laurent Le Goff Date: Fri, 20 May 2011 17:35:40 +0200 Subject: [PATCH] reduce number of points with distance_threshold parameter deduce from flattening_threshold --- draw2d/curve/Makefile | 2 + draw2d/curve/cubic_float64.go | 141 +-- draw2d/curve/cubic_float64_others.go | 1400 +++++++++++++------------- draw2d/curve/curve_test.go | 524 +++++----- draw2d/curve/quad_float64.go | 114 ++- 5 files changed, 1100 insertions(+), 1081 deletions(-) diff --git a/draw2d/curve/Makefile b/draw2d/curve/Makefile index f7939a1..15ceee0 100644 --- a/draw2d/curve/Makefile +++ b/draw2d/curve/Makefile @@ -6,4 +6,6 @@ GOFILES=\ quad_float64.go\ cubic_float64_others.go\ + + include $(GOROOT)/src/Make.pkg diff --git a/draw2d/curve/cubic_float64.go b/draw2d/curve/cubic_float64.go index a707bc3..b972c74 100644 --- a/draw2d/curve/cubic_float64.go +++ b/draw2d/curve/cubic_float64.go @@ -1,66 +1,75 @@ -// Copyright 2010 The draw2d Authors. All rights reserved. -// created: 17/05/2011 by Laurent Le Goff -package curve - -import ( - "math" -) - -const ( - CurveRecursionLimit = 32 -) - -type CubicCurveFloat64 struct { - X1, Y1, X2, Y2, X3, Y3, X4, Y4 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 - //---------------------- - c1.X1, c1.Y1 = c.X1, c.Y1 - c2.X4, c2.Y4 = c.X4, c.Y4 - c1.X2 = (c.X1 + c.X2) / 2 - c1.Y2 = (c.Y1 + c.Y2) / 2 - x23 = (c.X2 + c.X3) / 2 - y23 = (c.Y2 + c.Y3) / 2 - c2.X3 = (c.X3 + c.X4) / 2 - c2.Y3 = (c.Y3 + c.Y4) / 2 - c1.X3 = (c1.X2 + x23) / 2 - c1.Y3 = (c1.Y2 + y23) / 2 - c2.X2 = (x23 + c2.X3) / 2 - c2.Y2 = (y23 + c2.Y3) / 2 - c1.X4 = (c1.X3 + c2.X2) / 2 - c1.Y4 = (c1.Y3 + c2.Y2) / 2 - c2.X1, c2.Y1 = c1.X4, c1.Y4 - return -} - -func (curve *CubicCurveFloat64) Segment(t LineTracer, flattening_threshold float64) { - var curves [CurveRecursionLimit]CubicCurveFloat64 - curves[0] = *curve - i := 0 - // current curve - var c *CubicCurveFloat64 - var dx, dy, d2, d3 float64 - for i >= 0 { - c = &curves[i] - dx = c.X4 - c.X1 - dy = c.Y4 - c.Y1 - - d2 = math.Fabs(((c.X2-c.X4)*dy - (c.Y2-c.Y4)*dx)) - d3 = math.Fabs(((c.X3-c.X4)*dy - (c.Y3-c.Y4)*dx)) - - if (d2+d3)*(d2+d3) < flattening_threshold*(dx*dx+dy*dy) || i == len(curves)-1 { - t.LineTo(c.X4, c.Y4) - i-- - } else { - // second half of bezier go lower onto the stack - c.Subdivide(&curves[i+1], &curves[i]) - i++ - } - } -} +// Copyright 2010 The draw2d Authors. All rights reserved. +// created: 17/05/2011 by Laurent Le Goff +package curve + +import ( + "math" +) + +const ( + CurveRecursionLimit = 32 +) + +type CubicCurveFloat64 struct { + X1, Y1, X2, Y2, X3, Y3, X4, Y4 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 + //---------------------- + c1.X1, c1.Y1 = c.X1, c.Y1 + c2.X4, c2.Y4 = c.X4, c.Y4 + c1.X2 = (c.X1 + c.X2) / 2 + c1.Y2 = (c.Y1 + c.Y2) / 2 + x23 = (c.X2 + c.X3) / 2 + y23 = (c.Y2 + c.Y3) / 2 + c2.X3 = (c.X3 + c.X4) / 2 + c2.Y3 = (c.Y3 + c.Y4) / 2 + c1.X3 = (c1.X2 + x23) / 2 + c1.Y3 = (c1.Y2 + y23) / 2 + c2.X2 = (x23 + c2.X3) / 2 + c2.Y2 = (y23 + c2.Y3) / 2 + c1.X4 = (c1.X3 + c2.X2) / 2 + c1.Y4 = (c1.Y3 + c2.Y2) / 2 + c2.X1, c2.Y1 = c1.X4, c1.Y4 + return +} + +func (curve *CubicCurveFloat64) Segment(t LineTracer, flattening_threshold float64) { + var curves [CurveRecursionLimit]CubicCurveFloat64 + curves[0] = *curve + i := 0 + // current curve + var c *CubicCurveFloat64 + + var dx, dy, d2, d3 float64 + var lx, ly float64 + distance_threshold := flattening_threshold * 5 + lx, ly = curve.X1, curve.Y1 + + for i >= 0 { + c = &curves[i] + dx = c.X4 - c.X1 + dy = c.Y4 - c.Y1 + + d2 = math.Fabs(((c.X2-c.X4)*dy - (c.Y2-c.Y4)*dx)) + d3 = math.Fabs(((c.X3-c.X4)*dy - (c.Y3-c.Y4)*dx)) + + if (d2+d3)*(d2+d3) < flattening_threshold*(dx*dx+dy*dy) || i == len(curves)-1 { + if !(math.Fabs(lx - c.X4) < distance_threshold && math.Fabs(ly - c.Y4)< distance_threshold ) { + t.LineTo(c.X4, c.Y4) + lx, ly = c.X4, c.Y4 + } + + i-- + } else { + // second half of bezier go lower onto the stack + c.Subdivide(&curves[i+1], &curves[i]) + i++ + } + } +} diff --git a/draw2d/curve/cubic_float64_others.go b/draw2d/curve/cubic_float64_others.go index 5537d60..20c8e79 100644 --- a/draw2d/curve/cubic_float64_others.go +++ b/draw2d/curve/cubic_float64_others.go @@ -1,700 +1,700 @@ -// 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.X1 + 3*mu*mum1*mum1*c.X2 + 3*mu*mu*mum1*c.X3 + mu3*c.X4 - y = mum13*c.Y1 + 3*mu*mum1*mum1*c.Y2 + 3*mu*mu*mum1*c.Y3 + mu3*c.Y4 - return -} - -func (c *CubicCurveFloat64) SubdivideAt(c1, c2 *CubicCurveFloat64, t float64) (x23, y23 float64) { - inv_t := (1 - t) - c1.X1, c1.Y1 = c.X1, c.Y1 - c2.X4, c2.Y4 = c.X4, c.Y4 - - c1.X2 = inv_t*c.X1 + t*c.X2 - c1.Y2 = inv_t*c.Y1 + t*c.Y2 - - x23 = inv_t*c.X2 + t*c.X3 - y23 = inv_t*c.Y2 + t*c.Y3 - - c2.X3 = inv_t*c.X3 + t*c.X4 - c2.Y3 = inv_t*c.Y3 + t*c.Y4 - - c1.X3 = inv_t*c1.X2 + t*x23 - c1.Y3 = inv_t*c1.Y2 + t*y23 - - c2.X2 = inv_t*x23 + t*c2.X3 - c2.Y2 = inv_t*y23 + t*c2.Y3 - - c1.X4 = inv_t*c1.X3 + t*c2.X2 - c1.Y4 = inv_t*c1.Y3 + t*c2.Y2 - - c2.X1, c2.Y1 = c1.X4, c1.Y4 - return -} - -func (c *CubicCurveFloat64) EstimateDistance() float64 { - dx1 := c.X2 - c.X1 - dy1 := c.Y2 - c.Y1 - dx2 := c.X3 - c.X2 - dy2 := c.Y3 - c.Y2 - dx3 := c.X4 - c.X3 - dy3 := c.Y4 - c.Y3 - 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.X4, c.Y4) -} - -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.X4 - c.X1 - dy := c.Y4 - c.Y1 - - d2 := math.Fabs(((c.X2-c.X4)*dy - (c.Y2-c.Y4)*dx)) - d3 := math.Fabs(((c.X3-c.X4)*dy - (c.Y3-c.Y4)*dx)) - - if (d2+d3)*(d2+d3) < flattening_threshold*(dx*dx+dy*dy) { - t.LineTo(c.X4, c.Y4) - 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.X4, c.Y4) -} - -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.X4 - c.X1 - dy := c.Y4 - c.Y1 - - d2 := math.Fabs(((c.X2-c.X4)*dy - (c.Y2-c.Y4)*dx)) - d3 := math.Fabs(((c.X3-c.X4)*dy - (c.Y3-c.Y4)*dx)) - switch { - case d2 <= CurveCollinearityEpsilon && d3 <= CurveCollinearityEpsilon: - // All collinear OR p1==p4 - //---------------------- - k := dx*dx + dy*dy - if k == 0 { - d2 = squareDistance(c.X1, c.Y1, c.X2, c.Y2) - d3 = squareDistance(c.X4, c.Y4, c.X3, c.Y3) - } else { - k = 1 / k - da1 := c.X2 - c.X1 - da2 := c.Y2 - c.Y1 - d2 = k * (da1*dx + da2*dy) - da1 = c.X3 - c.X1 - da2 = c.Y3 - c.Y1 - 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.X2, c.Y2, c.X1, c.Y1) - } else if d2 >= 1 { - d2 = squareDistance(c.X2, c.Y2, c.X4, c.Y4) - } else { - d2 = squareDistance(c.X2, c.Y2, c.X1+d2*dx, c.Y1+d2*dy) - } - - if d3 <= 0 { - d3 = squareDistance(c.X3, c.Y3, c.X1, c.Y1) - } else if d3 >= 1 { - d3 = squareDistance(c.X3, c.Y3, c.X4, c.Y4) - } else { - d3 = squareDistance(c.X3, c.Y3, c.X1+d3*dx, c.Y1+d3*dy) - } - } - if d2 > d3 { - if d2 < distanceToleranceSquare { - t.LineTo(c.X2, c.Y2) - return - } - } else { - if d3 < distanceToleranceSquare { - t.LineTo(c.X3, c.Y3) - 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.Fabs(math.Atan2(c.Y4-c.Y3, c.X4-c.X3) - math.Atan2(c.Y3-c.Y2, c.X3-c.X2)) - if da1 >= math.Pi { - da1 = 2*math.Pi - da1 - } - - if da1 < angleTolerance { - t.LineTo(c.X2, c.Y2) - t.LineTo(c.X3, c.Y3) - return - } - - if cuspLimit != 0.0 { - if da1 > cuspLimit { - t.LineTo(c.X3, c.Y3) - 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.Fabs(math.Atan2(c.Y3-c.Y2, c.X3-c.X2) - math.Atan2(c.Y2-c.Y1, c.X2-c.X1)) - if da1 >= math.Pi { - da1 = 2*math.Pi - da1 - } - - if da1 < angleTolerance { - t.LineTo(c.X2, c.Y2) - t.LineTo(c.X3, c.Y3) - return - } - - if cuspLimit != 0.0 { - if da1 > cuspLimit { - t.LineTo(c.X2, c.Y2) - 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.Y3-c.Y2, c.X3-c.X2) - da1 := math.Fabs(k - math.Atan2(c.Y2-c.Y1, c.X2-c.X1)) - da2 := math.Fabs(math.Atan2(c.Y4-c.Y3, c.X4-c.X3) - 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.X2, c.Y2) - return - } - - if da2 > cuspLimit { - t.LineTo(c.X3, c.Y3) - 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.X4 - c.X1 - dy = c.Y4 - c.Y1 - - d2 = math.Fabs(((c.X2-c.X4)*dy - (c.Y2-c.Y4)*dx)) - d3 = math.Fabs(((c.X3-c.X4)*dy - (c.Y3-c.Y4)*dx)) - switch { - case i == len(curves)-1: - t.LineTo(c.X4, c.Y4) - i-- - continue - case d2 <= CurveCollinearityEpsilon && d3 <= CurveCollinearityEpsilon: - // All collinear OR p1==p4 - //---------------------- - k = dx*dx + dy*dy - if k == 0 { - d2 = squareDistance(c.X1, c.Y1, c.X2, c.Y2) - d3 = squareDistance(c.X4, c.Y4, c.X3, c.Y3) - } else { - k = 1 / k - da1 := c.X2 - c.X1 - da2 := c.Y2 - c.Y1 - d2 = k * (da1*dx + da2*dy) - da1 = c.X3 - c.X1 - da2 = c.Y3 - c.Y1 - 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.X2, c.Y2, c.X1, c.Y1) - } else if d2 >= 1 { - d2 = squareDistance(c.X2, c.Y2, c.X4, c.Y4) - } else { - d2 = squareDistance(c.X2, c.Y2, c.X1+d2*dx, c.Y1+d2*dy) - } - - if d3 <= 0 { - d3 = squareDistance(c.X3, c.Y3, c.X1, c.Y1) - } else if d3 >= 1 { - d3 = squareDistance(c.X3, c.Y3, c.X4, c.Y4) - } else { - d3 = squareDistance(c.X3, c.Y3, c.X1+d3*dx, c.Y1+d3*dy) - } - } - if d2 > d3 { - if d2 < distanceToleranceSquare { - t.LineTo(c.X2, c.Y2) - i-- - continue - } - } else { - if d3 < distanceToleranceSquare { - t.LineTo(c.X3, c.Y3) - 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.Fabs(math.Atan2(c.Y4-c.Y3, c.X4-c.X3) - math.Atan2(c.Y3-c.Y2, c.X3-c.X2)) - if da1 >= math.Pi { - da1 = 2*math.Pi - da1 - } - - if da1 < angleTolerance { - t.LineTo(c.X2, c.Y2) - t.LineTo(c.X3, c.Y3) - i-- - continue - } - - if cuspLimit != 0.0 { - if da1 > cuspLimit { - t.LineTo(c.X3, c.Y3) - 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.Fabs(math.Atan2(c.Y3-c.Y2, c.X3-c.X2) - math.Atan2(c.Y2-c.Y1, c.X2-c.X1)) - if da1 >= math.Pi { - da1 = 2*math.Pi - da1 - } - - if da1 < angleTolerance { - t.LineTo(c.X2, c.Y2) - t.LineTo(c.X3, c.Y3) - i-- - continue - } - - if cuspLimit != 0.0 { - if da1 > cuspLimit { - t.LineTo(c.X2, c.Y2) - 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.Y3-c.Y2, c.X3-c.X2) - da1 := math.Fabs(k - math.Atan2(c.Y2-c.Y1, c.X2-c.X1)) - da2 := math.Fabs(math.Atan2(c.Y4-c.Y3, c.X4-c.X3) - 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.X2, c.Y2) - i-- - continue - } - - if da2 > cuspLimit { - t.LineTo(c.X3, c.Y3) - i-- - continue - } - } - } - } - - // Continue subdivision - //---------------------- - curves[i+1], curves[i] = c1, c2 - i++ - } - t.LineTo(curve.X4, curve.Y4) -} - - -/********************** 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.X2 - c.X1 - dy := c.Y2 - c.Y1 - l2 := dx*dx + dy*dy - if l2 == 0 { - return 0 - } - l := math.Sqrt(l2) - r := (c.Y2 - c.Y1) / l - s := (c.X1 - c.X2) / l - u := (c.X2*c.Y1 - c.X1*c.Y2) / l - return math.Fabs(r*c.X3 + s*c.Y3 + u) -} - -// Find the number of inflection point -func (c *CubicCurveFloat64) numberOfInflectionPoints() int { - dx21 := (c.X2 - c.X1) - dy21 := (c.Y2 - c.Y1) - dx32 := (c.X3 - c.X2) - dy32 := (c.Y3 - c.Y2) - dx43 := (c.X4 - c.X3) - dy43 := (c.Y4 - c.Y3) - 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.X4 - c.X1 - dy = c.Y4 - c.Y1 - - d2 = math.Fabs(((c.X2-c.X4)*dy - (c.Y2-c.Y4)*dx)) - d3 = math.Fabs(((c.X3-c.X4)*dy - (c.Y3-c.Y4)*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.X4, c.Y4) - 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.X4 - b1.X1 - dy = b1.Y4 - b1.Y1 - - d2 = math.Fabs(((b1.X2-b1.X4)*dy - (b1.Y2-b1.Y4)*dx)) - d3 = math.Fabs(((b1.X3-b1.X4)*dy - (b1.Y3-b1.Y4)*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.X4, b1.Y4) - } - // 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.X1 + 3*curve.X2 - 3*curve.X3 + curve.X4) - bx := (3*curve.X1 - 6*curve.X2 + 3*curve.X3) - cx := (-3*curve.X1 + 3*curve.X2) - ay := (-curve.Y1 + 3*curve.Y2 - 3*curve.Y3 + curve.Y4) - by := (3*curve.Y1 - 6*curve.Y2 + 3*curve.Y3) - cy := (-3*curve.Y1 + 3*curve.Y2) - 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 -} +// 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.X1 + 3*mu*mum1*mum1*c.X2 + 3*mu*mu*mum1*c.X3 + mu3*c.X4 + y = mum13*c.Y1 + 3*mu*mum1*mum1*c.Y2 + 3*mu*mu*mum1*c.Y3 + mu3*c.Y4 + return +} + +func (c *CubicCurveFloat64) SubdivideAt(c1, c2 *CubicCurveFloat64, t float64) (x23, y23 float64) { + inv_t := (1 - t) + c1.X1, c1.Y1 = c.X1, c.Y1 + c2.X4, c2.Y4 = c.X4, c.Y4 + + c1.X2 = inv_t*c.X1 + t*c.X2 + c1.Y2 = inv_t*c.Y1 + t*c.Y2 + + x23 = inv_t*c.X2 + t*c.X3 + y23 = inv_t*c.Y2 + t*c.Y3 + + c2.X3 = inv_t*c.X3 + t*c.X4 + c2.Y3 = inv_t*c.Y3 + t*c.Y4 + + c1.X3 = inv_t*c1.X2 + t*x23 + c1.Y3 = inv_t*c1.Y2 + t*y23 + + c2.X2 = inv_t*x23 + t*c2.X3 + c2.Y2 = inv_t*y23 + t*c2.Y3 + + c1.X4 = inv_t*c1.X3 + t*c2.X2 + c1.Y4 = inv_t*c1.Y3 + t*c2.Y2 + + c2.X1, c2.Y1 = c1.X4, c1.Y4 + return +} + +func (c *CubicCurveFloat64) EstimateDistance() float64 { + dx1 := c.X2 - c.X1 + dy1 := c.Y2 - c.Y1 + dx2 := c.X3 - c.X2 + dy2 := c.Y3 - c.Y2 + dx3 := c.X4 - c.X3 + dy3 := c.Y4 - c.Y3 + 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.X4, c.Y4) +} + +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.X4 - c.X1 + dy := c.Y4 - c.Y1 + + d2 := math.Fabs(((c.X2-c.X4)*dy - (c.Y2-c.Y4)*dx)) + d3 := math.Fabs(((c.X3-c.X4)*dy - (c.Y3-c.Y4)*dx)) + + if (d2+d3)*(d2+d3) < flattening_threshold*(dx*dx+dy*dy) { + t.LineTo(c.X4, c.Y4) + 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.X4, c.Y4) +} + +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.X4 - c.X1 + dy := c.Y4 - c.Y1 + + d2 := math.Fabs(((c.X2-c.X4)*dy - (c.Y2-c.Y4)*dx)) + d3 := math.Fabs(((c.X3-c.X4)*dy - (c.Y3-c.Y4)*dx)) + switch { + case d2 <= CurveCollinearityEpsilon && d3 <= CurveCollinearityEpsilon: + // All collinear OR p1==p4 + //---------------------- + k := dx*dx + dy*dy + if k == 0 { + d2 = squareDistance(c.X1, c.Y1, c.X2, c.Y2) + d3 = squareDistance(c.X4, c.Y4, c.X3, c.Y3) + } else { + k = 1 / k + da1 := c.X2 - c.X1 + da2 := c.Y2 - c.Y1 + d2 = k * (da1*dx + da2*dy) + da1 = c.X3 - c.X1 + da2 = c.Y3 - c.Y1 + 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.X2, c.Y2, c.X1, c.Y1) + } else if d2 >= 1 { + d2 = squareDistance(c.X2, c.Y2, c.X4, c.Y4) + } else { + d2 = squareDistance(c.X2, c.Y2, c.X1+d2*dx, c.Y1+d2*dy) + } + + if d3 <= 0 { + d3 = squareDistance(c.X3, c.Y3, c.X1, c.Y1) + } else if d3 >= 1 { + d3 = squareDistance(c.X3, c.Y3, c.X4, c.Y4) + } else { + d3 = squareDistance(c.X3, c.Y3, c.X1+d3*dx, c.Y1+d3*dy) + } + } + if d2 > d3 { + if d2 < distanceToleranceSquare { + t.LineTo(c.X2, c.Y2) + return + } + } else { + if d3 < distanceToleranceSquare { + t.LineTo(c.X3, c.Y3) + 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.Fabs(math.Atan2(c.Y4-c.Y3, c.X4-c.X3) - math.Atan2(c.Y3-c.Y2, c.X3-c.X2)) + if da1 >= math.Pi { + da1 = 2*math.Pi - da1 + } + + if da1 < angleTolerance { + t.LineTo(c.X2, c.Y2) + t.LineTo(c.X3, c.Y3) + return + } + + if cuspLimit != 0.0 { + if da1 > cuspLimit { + t.LineTo(c.X3, c.Y3) + 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.Fabs(math.Atan2(c.Y3-c.Y2, c.X3-c.X2) - math.Atan2(c.Y2-c.Y1, c.X2-c.X1)) + if da1 >= math.Pi { + da1 = 2*math.Pi - da1 + } + + if da1 < angleTolerance { + t.LineTo(c.X2, c.Y2) + t.LineTo(c.X3, c.Y3) + return + } + + if cuspLimit != 0.0 { + if da1 > cuspLimit { + t.LineTo(c.X2, c.Y2) + 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.Y3-c.Y2, c.X3-c.X2) + da1 := math.Fabs(k - math.Atan2(c.Y2-c.Y1, c.X2-c.X1)) + da2 := math.Fabs(math.Atan2(c.Y4-c.Y3, c.X4-c.X3) - 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.X2, c.Y2) + return + } + + if da2 > cuspLimit { + t.LineTo(c.X3, c.Y3) + 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.X4 - c.X1 + dy = c.Y4 - c.Y1 + + d2 = math.Fabs(((c.X2-c.X4)*dy - (c.Y2-c.Y4)*dx)) + d3 = math.Fabs(((c.X3-c.X4)*dy - (c.Y3-c.Y4)*dx)) + switch { + case i == len(curves)-1: + t.LineTo(c.X4, c.Y4) + i-- + continue + case d2 <= CurveCollinearityEpsilon && d3 <= CurveCollinearityEpsilon: + // All collinear OR p1==p4 + //---------------------- + k = dx*dx + dy*dy + if k == 0 { + d2 = squareDistance(c.X1, c.Y1, c.X2, c.Y2) + d3 = squareDistance(c.X4, c.Y4, c.X3, c.Y3) + } else { + k = 1 / k + da1 := c.X2 - c.X1 + da2 := c.Y2 - c.Y1 + d2 = k * (da1*dx + da2*dy) + da1 = c.X3 - c.X1 + da2 = c.Y3 - c.Y1 + 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.X2, c.Y2, c.X1, c.Y1) + } else if d2 >= 1 { + d2 = squareDistance(c.X2, c.Y2, c.X4, c.Y4) + } else { + d2 = squareDistance(c.X2, c.Y2, c.X1+d2*dx, c.Y1+d2*dy) + } + + if d3 <= 0 { + d3 = squareDistance(c.X3, c.Y3, c.X1, c.Y1) + } else if d3 >= 1 { + d3 = squareDistance(c.X3, c.Y3, c.X4, c.Y4) + } else { + d3 = squareDistance(c.X3, c.Y3, c.X1+d3*dx, c.Y1+d3*dy) + } + } + if d2 > d3 { + if d2 < distanceToleranceSquare { + t.LineTo(c.X2, c.Y2) + i-- + continue + } + } else { + if d3 < distanceToleranceSquare { + t.LineTo(c.X3, c.Y3) + 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.Fabs(math.Atan2(c.Y4-c.Y3, c.X4-c.X3) - math.Atan2(c.Y3-c.Y2, c.X3-c.X2)) + if da1 >= math.Pi { + da1 = 2*math.Pi - da1 + } + + if da1 < angleTolerance { + t.LineTo(c.X2, c.Y2) + t.LineTo(c.X3, c.Y3) + i-- + continue + } + + if cuspLimit != 0.0 { + if da1 > cuspLimit { + t.LineTo(c.X3, c.Y3) + 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.Fabs(math.Atan2(c.Y3-c.Y2, c.X3-c.X2) - math.Atan2(c.Y2-c.Y1, c.X2-c.X1)) + if da1 >= math.Pi { + da1 = 2*math.Pi - da1 + } + + if da1 < angleTolerance { + t.LineTo(c.X2, c.Y2) + t.LineTo(c.X3, c.Y3) + i-- + continue + } + + if cuspLimit != 0.0 { + if da1 > cuspLimit { + t.LineTo(c.X2, c.Y2) + 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.Y3-c.Y2, c.X3-c.X2) + da1 := math.Fabs(k - math.Atan2(c.Y2-c.Y1, c.X2-c.X1)) + da2 := math.Fabs(math.Atan2(c.Y4-c.Y3, c.X4-c.X3) - 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.X2, c.Y2) + i-- + continue + } + + if da2 > cuspLimit { + t.LineTo(c.X3, c.Y3) + i-- + continue + } + } + } + } + + // Continue subdivision + //---------------------- + curves[i+1], curves[i] = c1, c2 + i++ + } + t.LineTo(curve.X4, curve.Y4) +} + + +/********************** 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.X2 - c.X1 + dy := c.Y2 - c.Y1 + l2 := dx*dx + dy*dy + if l2 == 0 { + return 0 + } + l := math.Sqrt(l2) + r := (c.Y2 - c.Y1) / l + s := (c.X1 - c.X2) / l + u := (c.X2*c.Y1 - c.X1*c.Y2) / l + return math.Fabs(r*c.X3 + s*c.Y3 + u) +} + +// Find the number of inflection point +func (c *CubicCurveFloat64) numberOfInflectionPoints() int { + dx21 := (c.X2 - c.X1) + dy21 := (c.Y2 - c.Y1) + dx32 := (c.X3 - c.X2) + dy32 := (c.Y3 - c.Y2) + dx43 := (c.X4 - c.X3) + dy43 := (c.Y4 - c.Y3) + 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.X4 - c.X1 + dy = c.Y4 - c.Y1 + + d2 = math.Fabs(((c.X2-c.X4)*dy - (c.Y2-c.Y4)*dx)) + d3 = math.Fabs(((c.X3-c.X4)*dy - (c.Y3-c.Y4)*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.X4, c.Y4) + 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.X4 - b1.X1 + dy = b1.Y4 - b1.Y1 + + d2 = math.Fabs(((b1.X2-b1.X4)*dy - (b1.Y2-b1.Y4)*dx)) + d3 = math.Fabs(((b1.X3-b1.X4)*dy - (b1.Y3-b1.Y4)*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.X4, b1.Y4) + } + // 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.X1 + 3*curve.X2 - 3*curve.X3 + curve.X4) + bx := (3*curve.X1 - 6*curve.X2 + 3*curve.X3) + cx := (-3*curve.X1 + 3*curve.X2) + ay := (-curve.Y1 + 3*curve.Y2 - 3*curve.Y3 + curve.Y4) + by := (3*curve.Y1 - 6*curve.Y2 + 3*curve.Y3) + cy := (-3*curve.Y1 + 3*curve.Y2) + 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/draw2d/curve/curve_test.go b/draw2d/curve/curve_test.go index 0c49d79..9439a8e 100644 --- a/draw2d/curve/curve_test.go +++ b/draw2d/curve/curve_test.go @@ -1,262 +1,262 @@ -package curve - -import ( - "testing" - "log" - "fmt" - "os" - "bufio" - "image" - "image/png" - "exp/draw" - "draw2d.googlecode.com/hg/draw2d/raster" -) - - -var ( - flattening_threshold float64 = 0.25 - testsCubicFloat64 = []CubicCurveFloat64{ - CubicCurveFloat64{100, 100, 200, 100, 100, 200, 200, 200}, - CubicCurveFloat64{100, 100, 300, 200, 200, 200, 300, 100}, - CubicCurveFloat64{100, 100, 0, 300, 200, 0, 300, 300}, - CubicCurveFloat64{150, 290, 10, 10, 290, 10, 150, 290}, - CubicCurveFloat64{10, 290, 10, 10, 290, 10, 290, 290}, - CubicCurveFloat64{100, 290, 290, 10, 10, 10, 200, 290}, - } - testsQuadFloat64 = []QuadCurveFloat64{ - QuadCurveFloat64{100, 100, 200, 100, 200, 200}, - QuadCurveFloat64{100, 100, 290, 200, 290, 100}, - QuadCurveFloat64{100, 100, 0, 290, 200, 290}, - QuadCurveFloat64{150, 290, 10, 10, 290, 290}, - QuadCurveFloat64{10, 290, 10, 10, 290, 290}, - QuadCurveFloat64{100, 290, 290, 10, 120, 290}, - } -) - -type Path struct { - points []float64 -} - -func (p *Path) LineTo(x, y float64) { - if len(p.points)+2 > cap(p.points) { - points := make([]float64, len(p.points)+2, len(p.points)+32) - copy(points, p.points) - p.points = points - } else { - p.points = p.points[0 : len(p.points)+2] - } - p.points[len(p.points)-2] = x - p.points[len(p.points)-1] = y -} - -func init() { - f, err := os.Create("_test.html") - if err != nil { - log.Println(err) - os.Exit(1) - } - defer f.Close() - 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))) - } - for i := 0; i < len(testsQuadFloat64); i++ { - f.Write([]byte(fmt.Sprintf("
\n
\n", i))) - } - f.Write([]byte("")) - -} - -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 image.Color, s ...float64) image.Image { - 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) - img.Set(x, y-1, c) - img.Set(x+1, y, c) - img.Set(x+1, y+1, c) - img.Set(x+1, y-1, c) - img.Set(x-1, y, c) - 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.X1, curve.Y1) - curve.SegmentRec(&p, flattening_threshold) - img := image.NewNRGBA(300, 300) - raster.PolylineBresenham(img, image.NRGBAColor{0xff, 0, 0, 0xff}, curve.X1, curve.Y1, curve.X2, curve.Y2, curve.X3, curve.Y3, curve.X4, curve.Y4) - raster.PolylineBresenham(img, image.Black, p.points...) - //drawPoints(img, image.NRGBAColor{0, 0, 0, 0xff}, curve.X1, curve.Y1, curve.X2, curve.Y2, curve.X3, curve.Y3, curve.X4, curve.Y4) - drawPoints(img, image.NRGBAColor{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() -} - -func TestCubicCurve(t *testing.T) { - for i, curve := range testsCubicFloat64 { - var p Path - p.LineTo(curve.X1, curve.Y1) - curve.Segment(&p, flattening_threshold) - img := image.NewNRGBA(300, 300) - raster.PolylineBresenham(img, image.NRGBAColor{0xff, 0, 0, 0xff}, curve.X1, curve.Y1, curve.X2, curve.Y2, curve.X3, curve.Y3, curve.X4, curve.Y4) - raster.PolylineBresenham(img, image.Black, p.points...) - //drawPoints(img, image.NRGBAColor{0, 0, 0, 0xff}, curve.X1, curve.Y1, curve.X2, curve.Y2, curve.X3, curve.Y3, curve.X4, curve.Y4) - drawPoints(img, image.NRGBAColor{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.X1, curve.Y1) - curve.AdaptiveSegmentRec(&p, 1, 0, 0) - img := image.NewNRGBA(300, 300) - raster.PolylineBresenham(img, image.NRGBAColor{0xff, 0, 0, 0xff}, curve.X1, curve.Y1, curve.X2, curve.Y2, curve.X3, curve.Y3, curve.X4, curve.Y4) - raster.PolylineBresenham(img, image.Black, p.points...) - //drawPoints(img, image.NRGBAColor{0, 0, 0, 0xff}, curve.X1, curve.Y1, curve.X2, curve.Y2, curve.X3, curve.Y3, curve.X4, curve.Y4) - drawPoints(img, image.NRGBAColor{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.X1, curve.Y1) - curve.AdaptiveSegment(&p, 1, 0, 0) - img := image.NewNRGBA(300, 300) - raster.PolylineBresenham(img, image.NRGBAColor{0xff, 0, 0, 0xff}, curve.X1, curve.Y1, curve.X2, curve.Y2, curve.X3, curve.Y3, curve.X4, curve.Y4) - raster.PolylineBresenham(img, image.Black, p.points...) - //drawPoints(img, image.NRGBAColor{0, 0, 0, 0xff}, curve.X1, curve.Y1, curve.X2, curve.Y2, curve.X3, curve.Y3, curve.X4, curve.Y4) - drawPoints(img, image.NRGBAColor{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.X1, curve.Y1) - curve.ParabolicSegment(&p, flattening_threshold) - img := image.NewNRGBA(300, 300) - raster.PolylineBresenham(img, image.NRGBAColor{0xff, 0, 0, 0xff}, curve.X1, curve.Y1, curve.X2, curve.Y2, curve.X3, curve.Y3, curve.X4, curve.Y4) - raster.PolylineBresenham(img, image.Black, p.points...) - //drawPoints(img, image.NRGBAColor{0, 0, 0, 0xff}, curve.X1, curve.Y1, curve.X2, curve.Y2, curve.X3, curve.Y3, curve.X4, curve.Y4) - drawPoints(img, image.NRGBAColor{0, 0, 0, 0xff}, p.points...) - savepng(fmt.Sprintf("_testParabolic%d.png", i), img) - log.Printf("Num of points: %d\n", len(p.points)) - } - fmt.Println() -} - - -func TestQuadCurve(t *testing.T) { - for i, curve := range testsQuadFloat64 { - var p Path - p.LineTo(curve.X1, curve.Y1) - curve.Segment(&p, flattening_threshold) - img := image.NewNRGBA(300, 300) - raster.PolylineBresenham(img, image.NRGBAColor{0xff, 0, 0, 0xff}, curve.X1, curve.Y1, curve.X2, curve.Y2, curve.X3, curve.Y3) - raster.PolylineBresenham(img, image.Black, p.points...) - drawPoints(img, image.NRGBAColor{0, 0, 0, 0xff}, p.points...) - savepng(fmt.Sprintf("_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.X1, curve.Y1) - 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.X1, curve.Y1) - 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.X1, curve.Y1) - 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.X1, curve.Y1) - 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.X1, curve.Y1) - 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.X1, curve.Y1) - curve.Segment(&p, flattening_threshold) - } - } -} +package curve + +import ( + "testing" + "log" + "fmt" + "os" + "bufio" + "image" + "image/png" + "exp/draw" + "draw2d.googlecode.com/hg/draw2d/raster" +) + + +var ( + flattening_threshold float64 = 0.25 + testsCubicFloat64 = []CubicCurveFloat64{ + CubicCurveFloat64{100, 100, 200, 100, 100, 200, 200, 200}, + CubicCurveFloat64{100, 100, 300, 200, 200, 200, 300, 100}, + CubicCurveFloat64{100, 100, 0, 300, 200, 0, 300, 300}, + CubicCurveFloat64{150, 290, 10, 10, 290, 10, 150, 290}, + CubicCurveFloat64{10, 290, 10, 10, 290, 10, 290, 290}, + CubicCurveFloat64{100, 290, 290, 10, 10, 10, 200, 290}, + } + testsQuadFloat64 = []QuadCurveFloat64{ + QuadCurveFloat64{100, 100, 200, 100, 200, 200}, + QuadCurveFloat64{100, 100, 290, 200, 290, 100}, + QuadCurveFloat64{100, 100, 0, 290, 200, 290}, + QuadCurveFloat64{150, 290, 10, 10, 290, 290}, + QuadCurveFloat64{10, 290, 10, 10, 290, 290}, + QuadCurveFloat64{100, 290, 290, 10, 120, 290}, + } +) + +type Path struct { + points []float64 +} + +func (p *Path) LineTo(x, y float64) { + if len(p.points)+2 > cap(p.points) { + points := make([]float64, len(p.points)+2, len(p.points)+32) + copy(points, p.points) + p.points = points + } else { + p.points = p.points[0 : len(p.points)+2] + } + p.points[len(p.points)-2] = x + p.points[len(p.points)-1] = y +} + +func init() { + f, err := os.Create("_test.html") + if err != nil { + log.Println(err) + os.Exit(1) + } + defer f.Close() + 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))) + } + for i := 0; i < len(testsQuadFloat64); i++ { + f.Write([]byte(fmt.Sprintf("
\n
\n", i))) + } + f.Write([]byte("")) + +} + +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 image.Color, s ...float64) image.Image { + /*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) + img.Set(x, y-1, c) + img.Set(x+1, y, c) + img.Set(x+1, y+1, c) + img.Set(x+1, y-1, c) + img.Set(x-1, y, c) + 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.X1, curve.Y1) + curve.SegmentRec(&p, flattening_threshold) + img := image.NewNRGBA(300, 300) + raster.PolylineBresenham(img, image.NRGBAColor{0xff, 0, 0, 0xff}, curve.X1, curve.Y1, curve.X2, curve.Y2, curve.X3, curve.Y3, curve.X4, curve.Y4) + raster.PolylineBresenham(img, image.Black, p.points...) + //drawPoints(img, image.NRGBAColor{0, 0, 0, 0xff}, curve.X1, curve.Y1, curve.X2, curve.Y2, curve.X3, curve.Y3, curve.X4, curve.Y4) + drawPoints(img, image.NRGBAColor{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() +} + +func TestCubicCurve(t *testing.T) { + for i, curve := range testsCubicFloat64 { + var p Path + p.LineTo(curve.X1, curve.Y1) + curve.Segment(&p, flattening_threshold) + img := image.NewNRGBA(300, 300) + raster.PolylineBresenham(img, image.NRGBAColor{0xff, 0, 0, 0xff}, curve.X1, curve.Y1, curve.X2, curve.Y2, curve.X3, curve.Y3, curve.X4, curve.Y4) + raster.PolylineBresenham(img, image.Black, p.points...) + //drawPoints(img, image.NRGBAColor{0, 0, 0, 0xff}, curve.X1, curve.Y1, curve.X2, curve.Y2, curve.X3, curve.Y3, curve.X4, curve.Y4) + drawPoints(img, image.NRGBAColor{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.X1, curve.Y1) + curve.AdaptiveSegmentRec(&p, 1, 0, 0) + img := image.NewNRGBA(300, 300) + raster.PolylineBresenham(img, image.NRGBAColor{0xff, 0, 0, 0xff}, curve.X1, curve.Y1, curve.X2, curve.Y2, curve.X3, curve.Y3, curve.X4, curve.Y4) + raster.PolylineBresenham(img, image.Black, p.points...) + //drawPoints(img, image.NRGBAColor{0, 0, 0, 0xff}, curve.X1, curve.Y1, curve.X2, curve.Y2, curve.X3, curve.Y3, curve.X4, curve.Y4) + drawPoints(img, image.NRGBAColor{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.X1, curve.Y1) + curve.AdaptiveSegment(&p, 1, 0, 0) + img := image.NewNRGBA(300, 300) + raster.PolylineBresenham(img, image.NRGBAColor{0xff, 0, 0, 0xff}, curve.X1, curve.Y1, curve.X2, curve.Y2, curve.X3, curve.Y3, curve.X4, curve.Y4) + raster.PolylineBresenham(img, image.Black, p.points...) + //drawPoints(img, image.NRGBAColor{0, 0, 0, 0xff}, curve.X1, curve.Y1, curve.X2, curve.Y2, curve.X3, curve.Y3, curve.X4, curve.Y4) + drawPoints(img, image.NRGBAColor{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.X1, curve.Y1) + curve.ParabolicSegment(&p, flattening_threshold) + img := image.NewNRGBA(300, 300) + raster.PolylineBresenham(img, image.NRGBAColor{0xff, 0, 0, 0xff}, curve.X1, curve.Y1, curve.X2, curve.Y2, curve.X3, curve.Y3, curve.X4, curve.Y4) + raster.PolylineBresenham(img, image.Black, p.points...) + //drawPoints(img, image.NRGBAColor{0, 0, 0, 0xff}, curve.X1, curve.Y1, curve.X2, curve.Y2, curve.X3, curve.Y3, curve.X4, curve.Y4) + drawPoints(img, image.NRGBAColor{0, 0, 0, 0xff}, p.points...) + savepng(fmt.Sprintf("_testParabolic%d.png", i), img) + log.Printf("Num of points: %d\n", len(p.points)) + } + fmt.Println() +} + + +func TestQuadCurve(t *testing.T) { + for i, curve := range testsQuadFloat64 { + var p Path + p.LineTo(curve.X1, curve.Y1) + curve.Segment(&p, flattening_threshold) + img := image.NewNRGBA(300, 300) + raster.PolylineBresenham(img, image.NRGBAColor{0xff, 0, 0, 0xff}, curve.X1, curve.Y1, curve.X2, curve.Y2, curve.X3, curve.Y3) + raster.PolylineBresenham(img, image.Black, p.points...) + drawPoints(img, image.NRGBAColor{0, 0, 0, 0xff}, p.points...) + savepng(fmt.Sprintf("_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.X1, curve.Y1) + 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.X1, curve.Y1) + 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.X1, curve.Y1) + 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.X1, curve.Y1) + 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.X1, curve.Y1) + 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.X1, curve.Y1) + curve.Segment(&p, flattening_threshold) + } + } +} diff --git a/draw2d/curve/quad_float64.go b/draw2d/curve/quad_float64.go index 4f776b7..29f651e 100644 --- a/draw2d/curve/quad_float64.go +++ b/draw2d/curve/quad_float64.go @@ -1,53 +1,61 @@ -// Copyright 2010 The draw2d Authors. All rights reserved. -// created: 17/05/2011 by Laurent Le Goff -package curve - -import ( - "math" -) - -type QuadCurveFloat64 struct { - X1, Y1, X2, Y2, X3, Y3 float64 -} - - -func (c *QuadCurveFloat64) Subdivide(c1, c2 *QuadCurveFloat64) { - // Calculate all the mid-points of the line segments - //---------------------- - c1.X1, c1.Y1 = c.X1, c.Y1 - c2.X3, c2.Y3 = c.X3, c.Y3 - c1.X2 = (c.X1 + c.X2) / 2 - c1.Y2 = (c.Y1 + c.Y2) / 2 - c2.X2 = (c.X2 + c.X3) / 2 - c2.Y2 = (c.Y2 + c.Y3) / 2 - c1.X3 = (c1.X2 + c2.X2) / 2 - c1.Y3 = (c1.Y2 + c2.Y2) / 2 - c2.X1, c2.Y1 = c1.X3, c1.Y3 - return -} - - -func (curve *QuadCurveFloat64) Segment(t LineTracer, flattening_threshold float64) { - var curves [CurveRecursionLimit]QuadCurveFloat64 - curves[0] = *curve - i := 0 - // current curve - var c *QuadCurveFloat64 - var dx, dy, d float64 - for i >= 0 { - c = &curves[i] - dx = c.X3 - c.X1 - dy = c.Y3 - c.Y1 - - d = math.Fabs(((c.X2-c.X3)*dy - (c.Y2-c.Y3)*dx)) - - if (d*d) < flattening_threshold*(dx*dx+dy*dy) || i == len(curves)-1 { - t.LineTo(c.X3, c.Y3) - i-- - } else { - // second half of bezier go lower onto the stack - c.Subdivide(&curves[i+1], &curves[i]) - i++ - } - } -} +// Copyright 2010 The draw2d Authors. All rights reserved. +// created: 17/05/2011 by Laurent Le Goff +package curve + +import ( + "math" +) + +type QuadCurveFloat64 struct { + X1, Y1, X2, Y2, X3, Y3 float64 +} + + +func (c *QuadCurveFloat64) Subdivide(c1, c2 *QuadCurveFloat64) { + // Calculate all the mid-points of the line segments + //---------------------- + c1.X1, c1.Y1 = c.X1, c.Y1 + c2.X3, c2.Y3 = c.X3, c.Y3 + c1.X2 = (c.X1 + c.X2) / 2 + c1.Y2 = (c.Y1 + c.Y2) / 2 + c2.X2 = (c.X2 + c.X3) / 2 + c2.Y2 = (c.Y2 + c.Y3) / 2 + c1.X3 = (c1.X2 + c2.X2) / 2 + c1.Y3 = (c1.Y2 + c2.Y2) / 2 + c2.X1, c2.Y1 = c1.X3, c1.Y3 + return +} + + +func (curve *QuadCurveFloat64) Segment(t LineTracer, flattening_threshold float64) { + var curves [CurveRecursionLimit]QuadCurveFloat64 + curves[0] = *curve + i := 0 + // current curve + var c *QuadCurveFloat64 + var dx, dy, d float64 + var lx, ly float64 + distance_threshold := flattening_threshold * 5 + lx, ly = curve.X1, curve.Y1 + + for i >= 0 { + c = &curves[i] + dx = c.X3 - c.X1 + dy = c.Y3 - c.Y1 + + d = math.Fabs(((c.X2-c.X3)*dy - (c.Y2-c.Y3)*dx)) + + if (d*d) < flattening_threshold*(dx*dx+dy*dy) || i == len(curves)-1 { + if !(math.Fabs(lx - c.X3) <= distance_threshold && math.Fabs(ly - c.Y3)<= distance_threshold ) { + t.LineTo(c.X3, c.Y3) + lx, ly = c.X3, c.Y3 + } + + i-- + } else { + // second half of bezier go lower onto the stack + c.Subdivide(&curves[i+1], &curves[i]) + i++ + } + } +}