go-gl / mathgl

A pure Go 3D math library.
BSD 3-Clause "New" or "Revised" License
554 stars 65 forks source link

possible error in BezierCurve2D()? #65

Closed quillaja closed 6 years ago

quillaja commented 6 years ago

I've been using the bezier curve functions in mgl64 lately, and I might have found an error in BezierCurve2D(). I noticed this when I plotted a cubic curve (4 control points) using the BezierCurve2D() and it did not produce the expected curve. I plotted the same curve using CubicBezierCurve2D() and got the expected result. In this case, I used the control points {4, 5}, {5.2, 6.3556}, {6.2, 4.2444}, {6, 3}.

Here's an image showing the 2 curves. Red is from BezierCurve2D() and green is the expected curve from CubicBezierCurve2D(). error_bezier

The code which produced this is below. I used the gosl/plt package to create the plot.

package main

import (
    "github.com/cpmech/gosl/plt"
    "github.com/go-gl/mathgl/mgl64"
)

func main() {

    // create 4 control points (for cubic curve)
    cPts := []mgl64.Vec2{{4, 5}, {5.2, 6.3556}, {6.2, 4.2444}, {6, 3}}
    // slices for points of bezier curves
    bezierCurve := []mgl64.Vec2{}
    cubicCurve := []mgl64.Vec2{}

    // use control points over t=[0,1] to get points for plot
    for t := 0.0; t <= 1.0; t += 0.01 {
        bezierCurve = append(bezierCurve,
            mgl64.BezierCurve2D(t, cPts))
        cubicCurve = append(cubicCurve,
            mgl64.CubicBezierCurve2D(t, cPts[0], cPts[1], cPts[2], cPts[3]))
    }

    // plot using gosl's plt package
    plt.Reset(true, nil)

    // with BezierCurve2D()
    plt.Plot(MapVec2(X, bezierCurve), MapVec2(Y, bezierCurve),
        &plt.A{C: "#FF0000", Closed: false})
    // with CubicBezierCurve2D()
    plt.Plot(MapVec2(X, cubicCurve), MapVec2(Y, cubicCurve),
        &plt.A{C: "#00FF00", Closed: false})

    for _, p := range cPts { // plot control points
        plt.PlotOne(p.X(), p.Y(), &plt.A{C: "#000000", M: "o"})
    }

    plt.Equal()
    // output image to same directory
    plt.Save(".", "error_bezier")
}

// functions that make it easier to convert a slice of Vec2 into 2 slices of
// x and y floats for plotting

// MapVec2 applies f to each element of slice.
func MapVec2(f func(v mgl64.Vec2) float64, slice []mgl64.Vec2) (rval []float64) {
    for _, v := range slice {
        rval = append(rval, f(v))
    }
    return
}

// X gets the x component of v
func X(v mgl64.Vec2) float64 { return v.X() }

// Y gets the y component of v
func Y(v mgl64.Vec2) float64 { return v.Y() }

Another example with the control points {0, 0}, {1, 3}, {3, 3}, {6, 0}: error_bezier

quillaja commented 6 years ago

Hi. I decided to look into this a little more today, and I think I've found the error. The error is in choose() and is due to integer division and the truncation that happens. Using a float internally to accumulate your result then casting to int at the end fixes the issue. Very easy fix. See below.

I haven't tested this with the actual bezier functions, but I noticed that CubicBezierCurve2D() explicitly lists the coefficients which BezierCurve2D() calculates via choose(). These coefficients are 1, 3, 3, 1, as expected. I copied the choose() function into the go playground and messed around with it. The current choose() produces 1, 3, 2, 1. This is clearly not the binomial coefficient, and explains why in the above examples the 3rd control point appears "ignored".

This error will also affect the other arbitrary-length control point bezier functions (eg BezierCurve3D). It would not be noticeable for a 3 point bezier curve where the coefficients are 1, 2, 1, but will be very noticeable as the number of control points increases.

Here is the simple fixed version.

func choose(n, k int) int {
    if k == 0 {
        return 1
    } else if n == 0 {
        return 0
    }
    result := float64(n - (k - 1))
    for i := 2; i <= k; i++ {
        result *= float64(n-(k-i)) / float64(i)
    }

    return int(result)
}
quillaja commented 6 years ago

I ran the above test plots again with the code in PR #71. I had to increase the width of the line for BezierCurve2D (the red one) so it was visible under the expected green one from CubicBezierCurve2D.

error_bezier_1

error_bezier_2