Resample curve into even length segments using C++

Question

What would be the best way of resampling a curve into even length segments using C++? What I have is a set of points that represents a 2d curve. In my example below I have a point struct with x and y components and a vector of points with test positions. Each pair of points represents a segment on the curve. Example resample curves are the images below. The red circles are the original positions the green circles are the target positions after the resample.

struct Point
{
    float x, y;
};


std::vector<Point> Points;
int numPoints = 5;

float positions[] = {
0.0350462,  -0.0589667,
0.0688311,  0.240896,
0.067369,   0.557199,
-0.024258,  0.715255,
0.0533231,  0.948694,
};

// Add points
int offset = 2;
for (int i =0; i < numPoints; i++)
{
    offset = i * 2;
    Point pt;
    pt.x = positions[offset];
    pt.y = positions[offset+1];
    Points.push_back(pt);

}

Answer 1

See if this works for you. Resampled points are equidistant from each other on the linear interpolation of the source vector's points.

#include <iostream>
#include <iomanip>
#include <vector>
#include <cmath>

struct Point {
    double x, y;
};

// Distance gives the Euclidean distance between two Points
double Distance(const Point& a, const Point& b) {
    const double dx = b.x - a.x;
    const double dy = b.y - a.y;
    const double lsq = dx*dx + dy*dy;
    return std::sqrt(lsq);
}

// LinearCurveLength calculates the total length of the linear
// interpolation through a vector of Points.  It is the sum of
// the Euclidean distances between all consecutive points in
// the vector.
double LinearCurveLength(std::vector<Point> const &points) {
    auto start = points.begin();
    if(start == points.end()) return 0;
    auto finish = start + 1;
    double sum = 0;
    while(finish != points.end()) {
        sum += Distance(*start, *finish);
        start = finish++;
    }
    return sum;
}

// Gives a vector of Points which are sampled as equally-spaced segments
// taken along the linear interpolation between points in the source.
// In general, consecutive points in the result will not be equidistant,
// because of a corner-cutting effect.
std::vector<Point> UniformLinearInterpolation(std::vector<Point> const &source, std::size_t target_count) {
    std::vector<Point> result;
    if(source.size() < 2 || target_count < 2) {
        // degenerate source vector or target_count value
        // for simplicity, this returns an empty result
        // but special cases may be handled when appropriate for the application
        return result;
    }
    // total_length is the total length along a linear interpolation
    // of the source points.
    const double total_length = LinearCurveLength(source);

    // segment_length is the length between result points, taken as
    // distance traveled between these points on a linear interpolation
    // of the source points.  The actual Euclidean distance between
    // points in the result vector can vary, and is always less than
    // or equal to segment_length.
    const double segment_length = total_length / (target_count - 1);

    // start and finish are the current source segment's endpoints
    auto start = source.begin();
    auto finish = start + 1;

    // src_segment_offset is the distance along a linear interpolation
    // of the source curve from its first point to the start of the current
    // source segment.
    double src_segment_offset = 0;

    // src_segment_length is the length of a line connecting the current
    // source segment's start and finish points.
    double src_segment_length = Distance(*start, *finish);

    // The first point in the result is the same as the first point
    // in the source.
    result.push_back(*start);

    for(std::size_t i=1; i<target_count-1; ++i) {
        // next_offset is the distance along a linear interpolation
        // of the source curve from its beginning to the location
        // of the i'th point in the result.
        // segment_length is multiplied by i here because iteratively
        // adding segment_length could accumulate error.
        const double next_offset = segment_length * i;

        // Check if next_offset lies inside the current source segment.
        // If not, move to the next source segment and update the
        // source segment offset and length variables.
        while(src_segment_offset + src_segment_length < next_offset) {
            src_segment_offset += src_segment_length;
            start = finish++;
            src_segment_length = Distance(*start, *finish);
        }
        // part_offset is the distance into the current source segment
        // associated with the i'th point's offset.
        const double part_offset = next_offset - src_segment_offset;

        // part_ratio is part_offset's normalized distance into the 
        // source segment. Its value is between 0 and 1,
        // where 0 locates the next point at "start" and 1
        // locates it at "finish".  In-between values represent a
        // weighted location between these two extremes.
        const double part_ratio = part_offset / src_segment_length;

        // Use part_ratio to calculate the next point's components
        // as weighted averages of components of the current
        // source segment's points.
        result.push_back({
            start->x + part_ratio * (finish->x - start->x),
            start->y + part_ratio * (finish->y - start->y)
        });
    }

    // The first and last points of the result are exactly
    // the same as the first and last points from the input,
    // so the iterated calculation above skips calculating
    // the last point in the result, which is instead copied
    // directly from the source vector here.
    result.push_back(source.back());
    return result;
}

int main() {
    std::vector<Point> points = {
        { 0.0350462,  -0.0589667},
        { 0.0688311,   0.240896 },
        { 0.067369,    0.557199 },
        {-0.024258,    0.715255 },
        { 0.0533231,   0.948694 }
    };

    std::cout << "Source Points:\n";
    for(const auto& point : points) {
        std::cout << std::setw(14) << point.x << "    " << std::setw(14) << point.y << '\n';
    }
    std::cout << '\n';
    auto interpolated = UniformLinearInterpolation(points, 7);
    std::cout << "Interpolated Points:\n";
    for(const auto& point : interpolated) {
        std::cout << std::setw(14) << point.x << "    " << std::setw(14) << point.y << '\n';
    }
    std::cout << '\n';
    std::cout << "Source linear interpolated length:           " << LinearCurveLength(points) << '\n';
    std::cout << "Interpolation's linear interpolated length:  " << LinearCurveLength(interpolated) << '\n';
}

Answer 2

For green points equidistant along the polyline:

The first run: walk through point list, calculate length of every segment and cumulative length up to current point. Pseudocode:

cumlen[0] = 0;
for (int i=1; i < numPoints; i++) {
   len = Sqrt((Point[i].x - Point[i-1].x)^2 + (Point[i].y - Point [i-1].y)^2)
   cumlen[i] = cumlen[i-1] + len; 
}

Now find length of every new piece

plen = cumlen[numpoints-1] / numpieces;

Now the second run - walk through point list and insert new points in appropriate segments.

i = 0;
for (ip=0; ip<numpieces; ip++) {
   curr = plen * ip;
   while cumlen[i+1] < curr
     i++;
   P[ip].x = Points[i].x + (curr - cumlen[i]) * (Points[i+1].x - Points[i].x) / 
                            (cumlen[i+1] - cumlen[i]);
    ..the same for y
}

Examples of real output for numpieces > numPoints and vice versa