How to equidistant resample a line (or curve)?

Question

I have a line l_1 given with a point series p_1,...,p_n. I now want a new line l_2 having k points: q_1,...,q_k. But for all i \in {1,...,k-1}: abs( q_i - q_i+1 ) = const, meaning the segments of l_2 are equidistant or uniform.

• k >= 2
• and p_1 and p_n should be in l_2.
• abs( p_i - p_i+1 ) not const

One solution is to approximate a line with a spline, and then subsample this again, to have uniform length segments then. Can I do any better? Is there any C++ code for that?

Ah, I've missed a specific detail: those q_i should be in l_1, meaning either they are on the line segments of l_1 or they are sample points of l_1.

Answer 1

Using a parametric function

You may define a piecewise parametric function:

 f[t_] := Piecewise[
      When x[i] <= t <= x[i + 1]

         f[t]= (y[i+1]-y[i]) (t - x[i]) / (x[i+1]-x[i]) + y[i], 

      For {i, 1 ... N};

Then select your points q, ideally spaced less than the minimum p[i+1]-p[i]

Finally sample f[q] at equal t intervals.

Sample result:

Here you can see the effect of reducing the interval size from the biggest to the smallest in the original sample:

You may evaluate the goodness of the approximation adding up the areas (integrating) between the original and re-sampled curves:

If you Plot the integrals for different interval sizes, you may decide what a good sampling is:

Just for the record, the code in Mathematica is:

a = 0;
p = Table[{   a = a + RandomReal[], RandomReal[]}, {10}];
f[t_, h_] := Piecewise[Table[{(h[[i + 1, 2]] - h[[i, 2]]) (t - h[[i, 1]]) /
                              (h[[i + 1, 1]] - h[[i, 1]]) + h[[i, 2]],
                       h[[i, 1]] <= t <= h[[i + 1, 1]]}, 
                       {i, 1, Length[h] - 1}]];

minSeg[h_] := Min[Table[Norm[h[[i, 1]] - h[[i + 1, 1]]], {i, Length[h] - 1}]];

newSegSize[h_] := (h[[Length@h, 1]] - h[[1, 1]])/
                  Ceiling[(h[[Length@h, 1]] - h[[1, 1]])/minSeg[h]]

qTable = Table[{t, f[t, p]}, {t, p[[1, 1]], p[[Length@p, 1]], newSegSize[p]}];

Edit: Answering your comment

Commented pgm code:

a = 0; (* Accumulator to ensure an increasing X Value*)

p = Table[{a = a + RandomReal[], 
    RandomReal[]}, {10}]; (*Generates 10 {x,y} Rnd points with \
                            increasing x Value*)

f[t_, h_] :=  (* Def. a PWise funct:
                Example of resulting function:
                     f[t,{{1,2},{2,2},{3,4}}]
                Returns teh following function definition:

                    Value          for Range
                     2             1<=t<=2
                 2+2*(-2+t)        2<=t<=3
                     0             True
              *)
  Piecewise[
   Table[{(h[[i + 1, 2]] - 
           h[[i, 2]]) (t - h[[i, 1]])/(h[[i + 1, 1]] - h[[i, 1]]) + h[[i, 2]],
           h[[i, 1]] <= t <= h[[i + 1, 1]]},
           {i, 1, Length[h] - 1}]];

  minSeg[h_] := (* Just lookup the min input point separation*)
               Min[Table[Norm[h[[i, 1]] - h[[i + 1, 1]]], {i, Length[h] - 1}]];

  newSegSize[h_] := (* Determine the new segment size for having
                       the full interval length as a multiple of the
                       segment size *)
                   (h[[Length@h, 1]] - h[[1, 1]])/
                    Ceiling[(h[[Length@h, 1]] - h[[1, 1]])/minSeg[h]]

   qTable =     (*Generates a table of points using the PW function *)
         Table[
               {t, f[t, p]},
               {t, p[[1, 1]], p[[Length@p, 1]],newSegSize[p]}];

   ListLinePlot[{qTable, p}, PlotStyle -> {Red, Blue}] (*Plot*)

Answer 2

It depends on your line points - what are they? If they define a smooth line, then resampling a cubic spline is a good bet.

Essentially, if you are making the points equidistant, you need to define what you want to see between the points - is smoothness more important than staying true to the original line? Is there a speed constraint?

As far as I can see, you are likely to end up with an iterative process here, because if your original points define a smooth line, it is not simple to calculate even the length of that line, let alone divide that up into equal pieces and determine the coordinates of these points.

If you use cubic splines, for each spline you should be able to calculate its length via the formula on Wikipedia's Arc Length article. However, it requires you to do integration - when you do numerical integration it is known as 'quadrature'. For a cubic (to calculate the length of a line segment between two original points), this should end up as a weighted sum of the coefficients of the cubic spline - particularly if you use Gaussian quadrature.

However, you could probably get a reasonable answer using piecewise cubic polynomials (generate a cubic polynomial from 2 points and the 2 points either side of them) and an iterative algorithm that improves guess values of xi to give equidistant points. But that's assuming you want speed rather than accuracy.