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 & 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 <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:
..the set of 256 low-frequency points obtained from one wavelet step...
...through to the 128 points obtained by two steps:
.
