<p><em><strong>Disclaimer</strong>: I've done some work in the past with wavelet decomposition of <strong>image</strong> 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</em></p>

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

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

<p>There are multiple different wavelets that can be used - the simplest being the Haar wavelet - but I think a linear (see <a href="https://9p.io/who/wim/papers/athome/" rel="nofollow noreferrer">Sweldens &amp; Schröder</a> 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).</p>

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

<blockquote>
  <p><em>Linear</em>:  <span class="math-container">$$H_i := P_{2i+1} - \frac{(P_{2i} +P_{2i+2})}{2}$$</span>
  <em>CR-cubic</em>:<span class="math-container">$$H_i := P_{2i+1} - \frac{9(P_{2i} +P_{2i+2})-(P_{2i-2} +P_{2i+4})}{16}$$</span></p>
</blockquote>

<p>...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)</p>

<blockquote>
  <p><span class="math-container">$$L_i := P_{2i} + \frac{(H_{i-1} + H_{i})}{4}$$</span></p>
</blockquote>

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

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

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

<p><strong><em>More notes</strong>: 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?</em></p>

<p>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.</p>

<pre><code>#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]);
    }   
}

</code></pre>

<p><br/>
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 (<em>which I hope will give a 'reasonable' approximation of Catmull-Rom splines without me having to write (or find) my own rendering code</em>)</p>

<p>This shows an original set of 512 points:
<a href="https://i.sstatic.net/t456F.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/t456F.png" alt="enter image description here"></a></p>

<p>..the set of 256 low-frequency points obtained from one wavelet step...
<a href="https://i.sstatic.net/ZktUd.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/ZktUd.png" alt="One Wavelet Step"></a></p>

<p>...through to the 128 points obtained by two steps:<a href="https://i.sstatic.net/OfdKt.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/OfdKt.png" alt="enter image description here"></a></p>

CC BY-SA 4.0
2019-10-04T14:24:55.960Z
It's hard to grasp the method from what you've said so that I can simply try it out.
dev_nut
CC BY-SA 4.0
2019-10-04T15:10:46.670Z
OK.  If I get a chance next week, I'll try making some diagrams.
Simon F
CC BY-SA 4.0
2019-10-04T15:29:17.993Z
Some code would be the best help. Maybe, using scipy (python) or eigen (c++), just so you don't have to write a 1000 lines.
dev_nut
CC BY-SA 4.0
2019-10-08T13:19:03.743Z
As requested, have added some C code and some example results. Also fixed a mistake in the 'update'  equation.
Simon F
CC BY-SA 4.0
2019-10-08T15:38:28.780Z
This looks really cool. Thanks! Essentially a $(log_2 n) - 1$ iterations need to be performed to get 4 points from this.
dev_nut
CC BY-SA 4.0
2019-10-08T16:21:56.537Z
Given that the initial end points aren't really being handled correctly in the above, I'm not sure just 4 points would be terribly useful.  The plot with just 4 (and also 8 points) looks like the outline of a nose! :-/.   I might need to fix the example code.
Simon F
CC BY-SA 4.0
2019-10-04T09:32:23.427Z
2019-11-08T12:07:13.197Z
2019-11-08T12:07:13.197Z
Simon F
<p>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). </p>

<p>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.</p>

<pre><code>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()
</code></pre>

<p>I'm now thinking of doing a robust regression approach now <a href="https://scipy-cookbook.readthedocs.io/items/robust_regression.html" rel="nofollow noreferrer">(https://scipy-cookbook.readthedocs.io/items/robust_regression.html</a>). 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. </p>

<p>Any help is appreciated with a relevant code sample.</p>

<p>Result:
<a href="https://i.sstatic.net/9Qc1y.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/9Qc1y.png" alt="as"></a>.</p>

<p>Original:</p>

<p><a href="https://i.sstatic.net/j3hxf.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/j3hxf.png" alt="enter image description here"></a></p>

<curve>
Non least squares formulation to fit catmull rom spline
dev_nut
# Non least squares formulation to fit catmull rom spline

(Asked by: username_0 on 2019-10-02)

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](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:
[![as](https://i.sstatic.net/9Qc1y.png)](https://i.sstatic.net/9Qc1y.png).

Original:

[![enter image description here](https://i.sstatic.net/j3hxf.png)](https://i.sstatic.net/j3hxf.png)

**Tags:** curve


---

## 1 answers

### Answer by username_1 on 2019-10-04. Score: 1

***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 **n**D) 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 **2**R 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 & Schröder](https://9p.io/who/wim/papers/athome/) 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 <stdio.h>
#include <stdlib.h>
#include <assert.h>

#define IS_EVEN(X) (((X)&1)==0)
#define MAX(X,Y) ((X)>=(Y)?(X):(Y))
#define MIN(X,Y) ((X)<=(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) &&  NumPtsIn>=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 < 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 = &PointsIn[MAX(i, 3) - 3];
        const vec_t * pv_m1 = &PointsIn[i-1];
        const vec_t * pv_p1 = &PointsIn[MIN(i+1, NumPtsIn - 1)];
        const vec_t * pv_p3 = &PointsIn[MIN(i+3, NumPtsIn - 1)];
        const vec_t * pv    = &PointsIn[i];

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

        for(UINT j = 0; j < 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->pos[j] + pv_p1->pos[j]) - (1.0/16.0) * (pv_m3->pos[j] + pv_p3->pos[j]);

            //Compute the difference from actual pos
            pResult->pos[j] =pv->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 = &TheZeroVec;

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

        const vec_t * pv    = &PointsIn[i];

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

        for(UINT j = 0; j < NUM_DIMS; j++)
        {
            //Add in 1/2 of the average of the adjacent HF/Odd terms
            pResult->pos[j] = pv->pos[j] + (pHFToLeft->pos[j] + pHFToRight->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 < NUM_ADDITIONAL_PASSES; Pass++)
    {
        const UINT CurrentNumSrcPoints = NUM_ORIG_POINTS >> (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:
[![enter image description here](https://i.sstatic.net/t456F.png)](https://i.sstatic.net/t456F.png)

..the set of 256 low-frequency points obtained from one wavelet step...
[![One Wavelet Step](https://i.sstatic.net/ZktUd.png)](https://i.sstatic.net/ZktUd.png)

...through to the 128 points obtained by two steps:[![enter image description here](https://i.sstatic.net/OfdKt.png)](https://i.sstatic.net/OfdKt.png)

> **Comments:**
> **username_0**: It's hard to grasp the method from what you've said so that I can simply try it out.  
> **username_1**: OK.  If I get a chance next week, I'll try making some diagrams.  
> **username_0**: Some code would be the best help. Maybe, using scipy (python) or eigen (c++), just so you don't have to write a 1000 lines.  
> **username_1**: As requested, have added some C code and some example results. Also fixed a mistake in the 'update'  equation.  
> **username_0**: This looks really cool. Thanks! Essentially a $(log_2 n) - 1$ iterations need to be performed to get 4 points from this.  
> **username_1**: Given that the initial end points aren't really being handled correctly in the above, I'm not sure just 4 points would be terribly useful.  The plot with just 4 (and also 8 points) looks like the outline of a nose! :-/.   I might need to fix the example code.  

---

