derecursify casteljau

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){
+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)}
+	return &CubicCurveFloat64{x1, y1, x2, y2, x3, y3, x4, y4}
 }
 
-func (c *CubicCurveFloat64) addPoint(x, y float64) {
+//mu ranges from 0 to 1, start to end of curve
-	c.segments = append(c.segments, x, y)
+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
+	c1.x1, c1.y1 = c.x1, c.y1
-	y12 := (c.y1 + c.y2) / 2
+	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
+	c2.x3 = (c.x3 + c.x4) / 2
-	y34 := (c.y3 + c.y4) / 2
+	c2.y3 = (c.y3 + c.y4) / 2
-	x123 := (x12 + x23) / 2
+	c1.x3 = (c1.x2 + x23) / 2
-	y123 := (y12 + y23) / 2
+	c1.y3 = (c1.y2 + y23) / 2
-	x234 := (x23 + x34) / 2
+	c2.x2 = (x23 + c2.x3) / 2
-	y234 := (y23 + y34) / 2
+	c2.y2 = (y23 + c2.y3) / 2
-	x1234 := (x123 + x234) / 2
+	c1.x4 = (c1.x3 + c2.x2) / 2
-	y1234 := (y123 + y234) / 2
+	c1.y4 = (c1.y3 + c2.y2) / 2
-	c1 = &CubicCurveFloat64{c.x1, c.y1, x12, y12, x123, y123, x1234, y1234, c.segments}
+	c2.x1, c2.y1 = c1.x4, c1.y4
-	c2 = &CubicCurveFloat64{x1234, y1234, x234, y234, x34, y34, c.x4, c.y4, c.segments}
+}
-	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 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)
+	segments = segments[0 : len(segments)+2]
-	c.addPoint(c.x1, c.y1)
+	segments[len(segments)-2] = c.x1
-	c.segmentCasteljauRec()
+	segments[len(segments)-1] = c.y1
-	c.addPoint(c.x4, c.y4)
+	segments = c.segmentCasteljauRec(segments)
-	return c.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() {
+func (c *CubicCurveFloat64) segmentCasteljauRec(segments []float64) []float64 {
-
+	var c1, c2 CubicCurveFloat64
-	c1, c2 := c.Subdivide()
+	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) {
+	if (d2+d3)*(d2+d3) < flattening_threshold*(dx*dx+dy*dy) {
-		c.addPoint(c2.x1, c2.y1)
+		segments = segments[0 : len(segments)+2]
-		return
+		segments[len(segments)-2] = c2.x1
-	} else {
+		segments[len(segments)-1] = c2.y1
+		return segments
+	}
 	// Continue subdivision
 	//----------------------
-		c1.segmentCasteljauRec()
+	segments = c1.segmentCasteljauRec(segments)
-		c2.segments = c1.segments 
+	segments = c2.segmentCasteljauRec(segments)
-		c2.segmentCasteljauRec()
+	return segments
-		c.segments = c2.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
+}

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()
+		s = make([]float64, 0, numSegments)
+		s = cf64Test1.SegmentCasteljau(s)
 	}
-	log.Printf("Num of points: %d\n", len(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))
 }