How to generate the best curved fit for an unknown set of 2D points in C++

Question

I am trying to get the best fit for an unknown set of 2D points. The points are centered points of rivers, and they don't come in a certain order. I've tried to use polynomial regression but I don't know what is the best polynomial order for different sets of data.

I've also tried cubic spline, but I don't want a line through all the points I have, I want an approximation of the best fit line through the points.

I would like to get something like this even for lines that have more curves. This example is computed with polynomial regression, and it works fine.

Is there a way to do some smooth or regression algorithm that can get the best fit line even for a set of points like the following?

    PolynomialRegression<double> pol;
    static int polynomOrder = <whateverPolynomOrderFitsBetter>;
    double error = 0.005f;
    std::vector<double> coeffs;
    pol.fitIt(x, y, polynomOrder, coeffs);
    // get fitted values
    for(std::size_t i = 0; i < points.size(); i++)
    {
        int order = polynomOrder;
        long double yFitted = 0;
        while(order >= 0)
        {
            yFitted += (coeffs[order] * pow(points[i].x, order) + error);
            order --;
        }
        points[i].y = yFitted;
    }

In my implementation with an 35 polynomial order this is all I can get , and also changing the polynomial order with higher values turns into Nan values for the coefficients.

I'm not sure if this is the best approach I can have.