Non least squares formulation to fit catmull rom spline I have a set of unordered points that I'm getting from an image attached. I'd like to simply fit a parametrized curve such as a catmull-rom curve to with n control points (n = 4 to 10, and can be assumed to be known). 

I've explored ways of doing this and it seems there doesn't seem to be some library function that simply does this. I explore b-spline fitting (as I read that it's possible to convert b-spline to catmull-rom afterwards) with scipy that simply seems to unable to fit this properly.

import scipy.interpolate as interpolate
tck, u = interpolate.splprep([x, y], s=6)

N = 100
xmin, xmax = 0, 1.01
unew = np.linspace(xmin, xmax, N)
spline = interpolate.splev(unew, tck)

plt.plot(x, y, 'bo', label='Original points')
plt.plot(spline[0], spline[1], 'r', label='BSpline')
plt.show()


I'm now thinking of doing a robust regression approach now (https://scipy-cookbook.readthedocs.io/items/robust_regression.html). But, I'm lost on how to convert my set of [x, y] points to a parametric formulation for catmull-rom spline to write the loss function needed for least squares fit. 

Any help is appreciated with a relevant code sample.

Result:
.

Original:
Disclaimer: I've done some work in the past with wavelet decomposition of image data but it occurred to me that it may be applicable to your problem.  Admittedly, I don't actually know if the following will work well in practice but, given how simple the approach is, and, as no other answer has been put forward, it seems to warrant investigation

Let's assume you have a sequence of 2D (or nD) vertices, $\small[P_0,P_1,P_2,...P_{N-1}]$, and you want to reduce those to a sequence of just $M$ which form your Catmull-Rom control points. For the moment let's assume $\small N=2^kM$, for some small value of $\small k$. We will apply k wavelet decompositions to reduce N to M points.

For our purposes, a single wavelet decomposition step takes a sequence of 2R points, $\small P=[P_0,P_1,...P_{2R-1}]$, and replaces it with two sequences of R points each: a low frequency set, $\small L=[L_0,...L_{R-1}]$, and a high frequency set, $\small H=[H_0,...H_{R-1}]$. I suspect in your use case you might only be interested in the low frequency sequence(s).

There are multiple different wavelets that can be used - the simplest being the Haar wavelet - but I think a linear (see Sweldens &amp; Schröder for an easy introduction)  or perhaps a 'Catmull-Rom-based' cubic wavelet might be suitable (though I don't know if this appears in any literature).

For each wavelet decomposition step we perform two passes; a 'lifting' operation to create all the high frequency terms, i.e. either...


  Linear:  $$H_i := P_{2i+1} - \frac{(P_{2i} +P_{2i+2})}{2}$$
  CR-cubic:$$H_i := P_{2i+1} - \frac{9(P_{2i} +P_{2i+2})-(P_{2i-2} +P_{2i+4})}{16}$$


...followed by an 'update'  which filters some of the high frequency terms back into the 'even'/low frequency terms as a correction. (We'll just assume, perhaps sub-optimally, this is the same for both linear and cubic)


  $$L_i := P_{2i} + \frac{(H_{i-1} + H_{i})}{4}$$


Notes: Actually, in practice, with clever indexing, these calculations can be done 'in place' on the original array.
This has also glossed over what happens when exceeding the bounds of the arrays, e.g. at $P_{2N}$ and $L_{-1}$, but one approach is to simply repeat values and/or pad with zeros.

To do the next wavelet decomposition, you repeat the process with the low frequency data, e.g. "P:=L" and repeat.

I imagine that the final set of low frequency points, L, will form a reasonable set of control points for a piecewise Catmull-Rom curve.

More notes: In this I've assumed that all the original points were 'equally important' and probably 'reasonably evenly spread'. It could be that you may want to look at the magnitudes of some of the {H_i} terms and if some are relatively large (compared to their neighbours) perhaps take action, e.g. say, re-insert them into the lower frequency set?

As requested in the comments, I've put together a little bit of C code to do the "cubic" decomposition. For simplicity, this just uses separate arrays for each of the lower frequency results rather than, say, doing the calculation "in place" in the original array.

#include &lt;stdio.h&gt;
#include &lt;stdlib.h&gt;
#include &lt;assert.h&gt;

#define IS_EVEN(X) (((X)&amp;1)==0)
#define MAX(X,Y) ((X)&gt;=(Y)?(X):(Y))
#define MIN(X,Y) ((X)&lt;=(Y)?(X):(Y))

typedef unsigned int UINT;


#define NUM_DIMS (2)

typedef struct
{
    float pos[NUM_DIMS];
}vec_t;

/*
// Lift
//
// Applies a wavelet "lift" step, i.e. computes high frequency terms, by
// predicting using Catmull-Rom spline through even points, and computing
// differences from odd
*/
static void Lift(const vec_t  PointsIn[], 
                 const UINT NumPtsIn,
                 vec_t HFOut[])
{
    //this *example* assumes an even number of input values
    assert(IS_EVEN(NumPtsIn) &amp;&amp;  NumPtsIn&gt;=2);

    /*
    // step through odd values doing a "lifting" step, i.e. predicting 
    // a position from neighbourhood of even points and computing the difference
    */
    for(UINT i = 1; i &lt; NumPtsIn; i+= 2)
    {
        /*
        // get 2 neighbours on either side of this odd point. IF we step out of range, then
        // clamp to nearest valid point
        //
        // If you want more efficient code, move these decisions out of the loop
        */
        const vec_t * pv_m3 = &amp;PointsIn[MAX(i, 3) - 3];
        const vec_t * pv_m1 = &amp;PointsIn[i-1];
        const vec_t * pv_p1 = &amp;PointsIn[MIN(i+1, NumPtsIn - 1)];
        const vec_t * pv_p3 = &amp;PointsIn[MIN(i+3, NumPtsIn - 1)];
        const vec_t * pv    = &amp;PointsIn[i];

        vec_t *pResult = &amp;HFOut[i/2];

        for(UINT j = 0; j &lt; NUM_DIMS; j++)
        {
            //Predict using Catmull-Rom interpolation - 
            //i.e. mid point of CR spline through the even points {p[i-3], p[i-1], p[i+1], p[i+3]}
            const float PredictedPos = (9.0/16.0) * (pv_m1-&gt;pos[j] + pv_p1-&gt;pos[j]) - (1.0/16.0) * (pv_m3-&gt;pos[j] + pv_p3-&gt;pos[j]);

            //Compute the difference from actual pos
            pResult-&gt;pos[j] =pv-&gt;pos[j] - PredictedPos;
        }
    }
}

/*
// Update
//
// Applies a wavelet "Update" step, i.e. computes the low frequency terms, 
// by filtering in the 'errors' in the prediction (i.e the HF terms) into the
// 'even' vals
*/
static void Update(const vec_t  PointsIn[], 
                   const UINT NumPtsIn,
                   const vec_t HFIn[],
                   vec_t LFOut[])
{

    static const vec_t TheZeroVec = {0};

    //this *example* assumes an even number of input values
    assert(IS_EVEN(NumPtsIn));

    //init as if for the "i = -1" iteration
    const vec_t * pHFToRight = &amp;TheZeroVec;

    /*
    // step through even values doing an "update" step,
    */
    for(UINT i = 0; i &lt; NumPtsIn; i+= 2)
    {
        /*
        // move one spot over in HF data
        */
        const vec_t *pHFToLeft = pHFToRight;
        pHFToRight = &amp;HFIn[i/2];

        const vec_t * pv    = &amp;PointsIn[i];

        vec_t *pResult = &amp;LFOut[i/2];

        for(UINT j = 0; j &lt; NUM_DIMS; j++)
        {
            //Add in 1/2 of the average of the adjacent HF/Odd terms
            pResult-&gt;pos[j] = pv-&gt;pos[j] + (pHFToLeft-&gt;pos[j] + pHFToRight-&gt;pos[j]) * 0.25;
        }
    }
}

//Some test data ...
#define NUM_ORIG_POINTS (512u)
static const vec_t PointsIn[NUM_ORIG_POINTS]=
{
#include "fjord.h"
};


#define NUM_ADDITIONAL_PASSES (4)
extern int main(int argc, const char *argv[])
{

    vec_t HighFreq[NUM_ORIG_POINTS/2];
    vec_t LowFreqs[NUM_ADDITIONAL_PASSES+1][NUM_ORIG_POINTS/2];

    /*
    // Do initial wavelet pass
    */
    Lift(PointsIn, NUM_ORIG_POINTS, HighFreq);
    Update(PointsIn, NUM_ORIG_POINTS, HighFreq, LowFreqs[0]);

    /*
    // Do a few more wavelet passes
    */
    for(UINT Pass=0; Pass &lt; NUM_ADDITIONAL_PASSES; Pass++)
    {
        const UINT CurrentNumSrcPoints = NUM_ORIG_POINTS &gt;&gt; (Pass+1);

        Lift(LowFreqs[Pass],   CurrentNumSrcPoints, HighFreq);
        Update(LowFreqs[Pass], CurrentNumSrcPoints, HighFreq, LowFreqs[Pass+1]);
    }   
}




The above has been run on the outline on a section of a map (a section of a fjord apparently) which I've plotted using Excel's scatter plot tools (which I hope will give a 'reasonable' approximation of Catmull-Rom splines without me having to write (or find) my own rendering code)

This shows an original set of 512 points:


..the set of 256 low-frequency points obtained from one wavelet step...


...through to the 128 points obtained by two steps: