Clean curve package

This commit is contained in:
Laurent Le Goff 2015-04-22 19:07:03 +02:00
parent f8fb5a2052
commit 216d3f60dd
6 changed files with 65 additions and 864 deletions

View file

@ -1,12 +1,14 @@
// Copyright 2010 The draw2d Authors. All rights reserved. // Copyright 2010 The draw2d Authors. All rights reserved.
// created: 21/11/2010 by Laurent Le Goff // created: 21/11/2010 by Laurent Le Goff
package curve package curve
import ( import (
"math" "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 end := start + angle
clockWise := true clockWise := true
if angle < 0 { if angle < 0 {

View file

@ -1,5 +1,7 @@
// Copyright 2010 The draw2d Authors. All rights reserved. // Copyright 2010 The draw2d Authors. All rights reserved.
// created: 17/05/2011 by Laurent Le Goff // created: 17/05/2011 by Laurent Le Goff
// Package curve implements Bezier Curve Subdivision using De Casteljau's algorithm
package curve package curve
import ( import (
@ -10,38 +12,48 @@ const (
CurveRecursionLimit = 32 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 CubicCurveFloat64 [8]float64
type LineTracer interface { // Subdivide a Bezier cubic curve in 2 equivalents Bezier cubic curves.
LineTo(x, y float64) // 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
func (c *CubicCurveFloat64) Subdivide(c1, c2 *CubicCurveFloat64) (x23, y23 float64) {
// Calculate all the mid-points of the line segments
//----------------------
c1[0], c1[1] = c[0], c[1] 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] c2[6], c2[7] = c[6], c[7]
// Subdivide segment using midpoints
c1[2] = (c[0] + c[2]) / 2 c1[2] = (c[0] + c[2]) / 2
c1[3] = (c[1] + c[3]) / 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[4] = (c[4] + c[6]) / 2
c2[5] = (c[5] + c[7]) / 2 c2[5] = (c[5] + c[7]) / 2
c1[4] = (c1[2] + x23) / 2
c1[5] = (c1[3] + y23) / 2 c1[4] = (c1[2] + midX) / 2
c2[2] = (x23 + c2[4]) / 2 c1[5] = (c1[3] + midY) / 2
c2[3] = (y23 + c2[5]) / 2
c2[2] = (midX + c2[4]) / 2
c2[3] = (midY + c2[5]) / 2
c1[6] = (c1[4] + c2[2]) / 2 c1[6] = (c1[4] + c2[2]) / 2
c1[7] = (c1[5] + c2[3]) / 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] 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 var curves [CurveRecursionLimit]CubicCurveFloat64
curves[0] = *curve curves[0] = *curve
i := 0 i := 0
// current curve // current curve
var c *CubicCurveFloat64 var c *CubicCurveFloat64
@ -52,9 +64,10 @@ func (curve *CubicCurveFloat64) Segment(t LineTracer, flattening_threshold float
dx = c[6] - c[0] dx = c[6] - c[0]
dy = c[7] - c[1] dy = c[7] - c[1]
d2 = math.Abs(((c[2]-c[6])*dy - (c[3]-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)) 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 { if (d2+d3)*(d2+d3) < flattening_threshold*(dx*dx+dy*dy) || i == len(curves)-1 {
t.LineTo(c[6], c[7]) t.LineTo(c[6], c[7])
i-- i--

View file

@ -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
}

View file

@ -1,16 +1,15 @@
package curve package curve
import ( import (
"bufio"
"fmt" "fmt"
"image" "image"
"image/color" "image/color"
"image/draw" "image/draw"
"image/png"
"log" "log"
"os" "os"
"testing" "testing"
"github.com/llgcode/draw2d"
"github.com/llgcode/draw2d/raster" "github.com/llgcode/draw2d/raster"
) )
@ -51,7 +50,8 @@ func (p *Path) LineTo(x, y float64) {
} }
func init() { func init() {
f, err := os.Create("_test.html") os.Mkdir("test_results", 0666)
f, err := os.Create("test_results/_test.html")
if err != nil { if err != nil {
log.Println(err) log.Println(err)
os.Exit(1) os.Exit(1)
@ -60,7 +60,7 @@ func init() {
log.Printf("Create html viewer") log.Printf("Create html viewer")
f.Write([]byte("<html><body>")) f.Write([]byte("<html><body>"))
for i := 0; i < len(testsCubicFloat64); i++ { for i := 0; i < len(testsCubicFloat64); i++ {
f.Write([]byte(fmt.Sprintf("<div><img src='_testRec%d.png'/>\n<img src='_test%d.png'/>\n<img src='_testAdaptiveRec%d.png'/>\n<img src='_testAdaptive%d.png'/>\n<img src='_testParabolic%d.png'/>\n</div>\n", i, i, i, i, i))) f.Write([]byte(fmt.Sprintf("<div><img src='_test%d.png'/></div>\n", i)))
} }
for i := 0; i < len(testsQuadFloat64); i++ { for i := 0; i < len(testsQuadFloat64); i++ {
f.Write([]byte(fmt.Sprintf("<div><img src='_testQuad%d.png'/>\n</div>\n", i))) f.Write([]byte(fmt.Sprintf("<div><img src='_testQuad%d.png'/>\n</div>\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 { 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) x, y := int(s[i]+0.5), int(s[i+1]+0.5)
img.Set(x, y, c) img.Set(x, y, c)
img.Set(x, y+1, 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)
img.Set(x-1, y-1, c) img.Set(x-1, y-1, c)
}*/ }
return img 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()
}
func TestCubicCurve(t *testing.T) { func TestCubicCurve(t *testing.T) {
for i, curve := range testsCubicFloat64 { for i, curve := range testsCubicFloat64 {
var p Path var p Path
p.LineTo(curve[0], curve[1]) 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)) img := image.NewNRGBA(image.Rect(0, 0, 300, 300))
raster.PolylineBresenham(img, color.NRGBA{0xff, 0, 0, 0xff}, curve[:]...) raster.PolylineBresenham(img, color.NRGBA{0xff, 0, 0, 0xff}, curve[:]...)
raster.PolylineBresenham(img, image.Black, p.points...) raster.PolylineBresenham(img, image.Black, p.points...)
//drawPoints(img, image.NRGBAColor{0, 0, 0, 0xff}, curve[:]...) //drawPoints(img, image.NRGBAColor{0, 0, 0, 0xff}, curve[:]...)
drawPoints(img, color.NRGBA{0, 0, 0, 0xff}, p.points...) drawPoints(img, color.NRGBA{0, 0, 0, 0xff}, p.points...)
savepng(fmt.Sprintf("_test%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()
}
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)
log.Printf("Num of points: %d\n", len(p.points)) log.Printf("Num of points: %d\n", len(p.points))
} }
fmt.Println() fmt.Println()
@ -190,74 +106,24 @@ func TestQuadCurve(t *testing.T) {
for i, curve := range testsQuadFloat64 { for i, curve := range testsQuadFloat64 {
var p Path var p Path
p.LineTo(curve[0], curve[1]) 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)) img := image.NewNRGBA(image.Rect(0, 0, 300, 300))
raster.PolylineBresenham(img, color.NRGBA{0xff, 0, 0, 0xff}, curve[:]...) raster.PolylineBresenham(img, color.NRGBA{0xff, 0, 0, 0xff}, curve[:]...)
raster.PolylineBresenham(img, image.Black, p.points...) raster.PolylineBresenham(img, image.Black, p.points...)
//drawPoints(img, image.NRGBAColor{0, 0, 0, 0xff}, curve[:]...) //drawPoints(img, image.NRGBAColor{0, 0, 0, 0xff}, curve[:]...)
drawPoints(img, color.NRGBA{0, 0, 0, 0xff}, p.points...) 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)) log.Printf("Num of points: %d\n", len(p.points))
} }
fmt.Println() 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) { func BenchmarkCubicCurve(b *testing.B) {
for i := 0; i < b.N; i++ { for i := 0; i < b.N; i++ {
for _, curve := range testsCubicFloat64 { for _, curve := range testsCubicFloat64 {
p := Path{make([]float64, 0, 32)} p := Path{make([]float64, 0, 32)}
p.LineTo(curve[0], curve[1]) p.LineTo(curve[0], curve[1])
curve.Segment(&p, flattening_threshold) curve.Trace(&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)
} }
} }
} }

View file

@ -1,19 +1,25 @@
// Copyright 2010 The draw2d Authors. All rights reserved. // Copyright 2010 The draw2d Authors. All rights reserved.
// created: 17/05/2011 by Laurent Le Goff // created: 17/05/2011 by Laurent Le Goff
package curve package curve
import ( import (
"math" "math"
) )
//X1, Y1, X2, Y2, X3, Y3 float64 //x1, y1, cpx1, cpy2, x2, y2 float64
type QuadCurveFloat64 [6]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) { 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] 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] c2[4], c2[5] = c[4], c[5]
// Subdivide segment using midpoints
c1[2] = (c[0] + c[2]) / 2 c1[2] = (c[0] + c[2]) / 2
c1[3] = (c[1] + c[3]) / 2 c1[3] = (c[1] + c[3]) / 2
c2[2] = (c[2] + c[4]) / 2 c2[2] = (c[2] + c[4]) / 2
@ -24,7 +30,10 @@ func (c *QuadCurveFloat64) Subdivide(c1, c2 *QuadCurveFloat64) {
return 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 var curves [CurveRecursionLimit]QuadCurveFloat64
curves[0] = *curve curves[0] = *curve
i := 0 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)) 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 { if (d*d) < flattening_threshold*(dx*dx+dy*dy) || i == len(curves)-1 {
t.LineTo(c[4], c[5]) t.LineTo(c[4], c[5])
i-- i--

6
curve/tracer.go Normal file
View file

@ -0,0 +1,6 @@
package curve
// LineTracer is an interface that help segmenting curve into small lines
type LineTracer interface {
LineTo(x, y float64)
}