From 9606b186e02bcad8d1278a627df4c87be4d9a9bf Mon Sep 17 00:00:00 2001 From: Laurent Le Goff Date: Wed, 18 May 2011 18:31:32 +0200 Subject: [PATCH] derecursify casteljau --- draw2d/curve/curve_float64.go | 156 +++++++++++++++++++++++++--------- draw2d/curve/curve_test.go | 77 ++++++++++++++++- 2 files changed, 190 insertions(+), 43 deletions(-) diff --git a/draw2d/curve/curve_float64.go b/draw2d/curve/curve_float64.go index c34c3ae..f08a9f5 100644 --- a/draw2d/curve/curve_float64.go +++ b/draw2d/curve/curve_float64.go @@ -7,56 +7,102 @@ import ( ) var ( - m_distance_tolerance float64 = 0.25 + flattening_threshold float64 = 0.25 ) type CubicCurveFloat64 struct { x1, y1, x2, y2, x3, y3, x4, y4 float64 - segments []float64 } -func NewCubicCurveFloat64(x1, y1, x2, y2, x3, y3, x4, y4 float64) (*CubicCurveFloat64){ - return &CubicCurveFloat64{x1, y1, x2, y2, x3, y3, x4, y4, make([]float64, 0, 2)} +func NewCubicCurveFloat64(x1, y1, x2, y2, x3, y3, x4, y4 float64) *CubicCurveFloat64 { + return &CubicCurveFloat64{x1, y1, x2, y2, x3, y3, x4, y4} } -func (c *CubicCurveFloat64) addPoint(x, y float64) { - c.segments = append(c.segments, x, y) +//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) Subdivide() (c1, c2 *CubicCurveFloat64) { +func (c *CubicCurveFloat64) SubdivideAt(c1, c2 *CubicCurveFloat64, t 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 +} + +func (c *CubicCurveFloat64) Subdivide(c1, c2 *CubicCurveFloat64) { // Calculate all the mid-points of the line segments //---------------------- - x12 := (c.x1 + c.x2) / 2 - y12 := (c.y1 + c.y2) / 2 + 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 - x34 := (c.x3 + c.x4) / 2 - y34 := (c.y3 + c.y4) / 2 - x123 := (x12 + x23) / 2 - y123 := (y12 + y23) / 2 - x234 := (x23 + x34) / 2 - y234 := (y23 + y34) / 2 - x1234 := (x123 + x234) / 2 - y1234 := (y123 + y234) / 2 - c1 = &CubicCurveFloat64{c.x1, c.y1, x12, y12, x123, y123, x1234, y1234, c.segments} - c2 = &CubicCurveFloat64{x1234, y1234, x234, y234, x34, y34, c.x4, c.y4, c.segments} - return + 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 +} + +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 Casteljau subdivision // and computing minimal distance tolerance -func (c *CubicCurveFloat64) SegmentCasteljau() []float64{ +func (c *CubicCurveFloat64) SegmentCasteljauRec(segments []float64) []float64 { // reinit segments - c.segments = make([]float64, 0, 2) - c.addPoint(c.x1, c.y1) - c.segmentCasteljauRec() - c.addPoint(c.x4, c.y4) - return c.segments + segments = segments[0 : len(segments)+2] + segments[len(segments)-2] = c.x1 + segments[len(segments)-1] = c.y1 + segments = c.segmentCasteljauRec(segments) + segments = segments[0 : len(segments)+2] + segments[len(segments)-2] = c.x4 + segments[len(segments)-1] = c.y4 + return segments } -func (c *CubicCurveFloat64) segmentCasteljauRec() { - - c1, c2 := c.Subdivide() +func (c *CubicCurveFloat64) segmentCasteljauRec(segments []float64) []float64 { + var c1, c2 CubicCurveFloat64 + c.Subdivide(&c1, &c2) // Try to approximate the full cubic curve by a single straight line //------------------ @@ -66,15 +112,47 @@ func (c *CubicCurveFloat64) segmentCasteljauRec() { 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) < m_distance_tolerance*(dx*dx+dy*dy) { - c.addPoint(c2.x1, c2.y1) - return - } else { - // Continue subdivision - //---------------------- - c1.segmentCasteljauRec() - c2.segments = c1.segments - c2.segmentCasteljauRec() - c.segments = c2.segments + if (d2+d3)*(d2+d3) < flattening_threshold*(dx*dx+dy*dy) { + segments = segments[0 : len(segments)+2] + segments[len(segments)-2] = c2.x1 + segments[len(segments)-1] = c2.y1 + return segments } + // Continue subdivision + //---------------------- + segments = c1.segmentCasteljauRec(segments) + segments = c2.segmentCasteljauRec(segments) + return segments } + +func (curve *CubicCurveFloat64) SegmentCasteljau(segments []float64) ([]float64) { + var curves [32]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 { + segments = segments[0 : len(segments)+2] + segments[len(segments)-2] = c.x1 + segments[len(segments)-1] = c.y1 + i--; + } else { + // second half of bezier go lower onto the stack + c.Subdivide(&curves[i+1], &curves[i]) + i++; + } + } + segments = segments[0 : len(segments)+2] + segments[len(segments)-2] = curve.x1 + segments[len(segments)-1] = curve.y1 + return segments +} \ No newline at end of file diff --git a/draw2d/curve/curve_test.go b/draw2d/curve/curve_test.go index 787f51b..cf621f5 100644 --- a/draw2d/curve/curve_test.go +++ b/draw2d/curve/curve_test.go @@ -8,12 +8,81 @@ import ( var ( cf64Test1 = NewCubicCurveFloat64(100, 100, 200, 100, 100, 200, 200, 200) + cf64Test2 = NewCubicCurveFloat64(100, 100, 300, 200, 200, 200, 200, 100) ) +func TestCubicCurveCasteljauRecTest1(t *testing.T) { + var s []float64 + d := cf64Test1.EstimateDistance() + log.Printf("Distance estimation: %f\n", d) + numSegments := int(d * 0.25) + log.Printf("Max segments estimation: %d\n", numSegments) + s = make([]float64, 0, numSegments) + s = cf64Test1.SegmentCasteljauRec(s) + log.Printf("points: %v\n", s) + log.Printf("Num of points: %d\n", len(s)) +} + +func TestCubicCurveCasteljauTest1(t *testing.T) { + var s []float64 + d := cf64Test1.EstimateDistance() + log.Printf("Distance estimation: %f\n", d) + numSegments := int(d * 0.25) + log.Printf("Max segments estimation: %d\n", numSegments) + s = make([]float64, 0, numSegments) + s = cf64Test1.SegmentCasteljau(s) + log.Printf("points: %v\n", s) + log.Printf("Num of points: %d\n", len(s)) +} + + +func BenchmarkCubicCurveCasteljauRecTest1(b *testing.B) { + var s []float64 + d := cf64Test1.EstimateDistance() + log.Printf("Distance estimation: %f\n", d) + numSegments := int(d * 0.25) + log.Printf("Max segments estimation: %d\n", numSegments) + for i := 0; i < b.N; i++ { + s = make([]float64, 0, numSegments) + s = cf64Test1.SegmentCasteljauRec(s) + } + log.Printf("Num of points: %d\n", len(s)) +} + +func BenchmarkCubicCurveCasteljauRecTest2(b *testing.B) { + var s []float64 + d := cf64Test1.EstimateDistance() + log.Printf("Distance estimation: %f\n", d) + numSegments := int(d * 0.25) + log.Printf("Max segments estimation: %d\n", numSegments) + for i := 0; i < b.N; i++ { + s = make([]float64, 0, numSegments) + s = cf64Test2.SegmentCasteljauRec(s) + } + log.Printf("Num of points: %d\n", len(s)) +} func BenchmarkCubicCurveCasteljauTest1(b *testing.B) { var s []float64 + d := cf64Test1.EstimateDistance() + log.Printf("Distance estimation: %f\n", d) + numSegments := int(d * 0.25) + log.Printf("Max segments estimation: %d\n", numSegments) for i := 0; i < b.N; i++ { - s = cf64Test1.SegmentCasteljau() - } - log.Printf("Num of points: %d\n", len(s)); -} \ No newline at end of file + s = make([]float64, 0, numSegments) + s = cf64Test1.SegmentCasteljau(s) + } + log.Printf("Num of points: %d\n", len(s)) +} + +func BenchmarkCubicCurveCasteljauTest2(b *testing.B) { + var s []float64 + d := cf64Test1.EstimateDistance() + log.Printf("Distance estimation: %f\n", d) + numSegments := int(d * 0.25) + log.Printf("Max segments estimation: %d\n", numSegments) + for i := 0; i < b.N; i++ { + s = make([]float64, 0, numSegments) + s = cf64Test2.SegmentCasteljau(s) + } + log.Printf("Num of points: %d\n", len(s)) +}