缺少数据的 2D 图像中的曲线拟合

Question

我在图像中有分割轮廓，我想将它们连接起来，假设它们是参数化的曲线。我正在考虑一些曲线拟合或曲线 跟踪 - 但我不知道如何实现它。

我正好有 1 px 宽的分割轮廓：

import itertools
import cv2

# Get all pixel directions
directions = list(itertools.product([-1, 0, 1], repeat=2))
directions.remove((0, 0))

# Load image
img = cv2.imread("image.png", cv2.IMREAD_COLOR)
gray = cv2.cvtColor(img, cv2.COLOR_RGB2GRAY)

cv2.imshow("Input", img)

# Find contours
_, thresh = cv2.threshold(gray, 25, 255, cv2.THRESH_BINARY)
contours, _ = cv2.findContours(thresh, cv2.RETR_LIST, cv2.CHAIN_APPROX_NONE)

我知道，所有轮廓的终点在哪里：

for cnt in contours:   
    ends = [] 

    # Get ends of contour
    # - first pixel is not always the first pixel of open contour
    # - last pixel is not mostly the last pixel in contour

    for pix in cnt:
        pixel = tuple(pix[0])
        pixel_x, pixel_y = pixel

        total_neighbours = 0

        # Count pixel neighbours
        for plus_x, plus_y in directions:
            current_x = pixel_x + plus_x
            current_y = pixel_y + plus_y

            if current_y < 0 or current_y >= gray.shape[0]:
                continue

            if current_x < 0 or current_x >= gray.shape[1]:
                continue

            if gray[current_y][current_x]:
                total_neighbours += 1

        # If it is end pixel
        if total_neighbours == 1:
            ends.append(pixel)
            cv2.circle(img, pixel, 3, (0, 255, 255), 1)

作为一个人类，我知道曲线在哪里，每条独特曲线的轮廓是什么：

为了更好的人类想象力，有原始图像：

但我不知道，如何以编程方式连接这些分割轮廓。我尝试使用一阶和二阶导数进行预测，但还不够好：

    # Make contour "predictions" by first derivative
    # - go from end to second end and calculate first derivative
    # - "predict" on second end contour connection

    for end in ends:
        pixel_x, pixel_y = end

        def predict(pixel_x, pixel_y, was):
            for plus_x, plus_y in directions:
                current_x = pixel_x + plus_x
                current_y = pixel_y + plus_y

                if current_y < 0 or current_y >= gray.shape[0]:
                    continue

                if current_x < 0 or current_x >= gray.shape[1]:
                    continue

                if (current_x, current_y) not in ends:
                    if (current_x, current_y) not in was:
                        if gray[current_y][current_x]:
                            was.append((current_x, current_y))
                            predict(current_x, current_y, was)

                else:
                    derivatives_x = []
                    derivatives_y = []

                    # Calculate derivative
                    for pix in range(len(was) - 3):
                        x1, y1 = was[pix]
                        x2, y2 = was[pix + 3]

                        derivatives_x.append(x1 - x2)
                        derivatives_y.append(y1 - y2)

                    if not derivatives_x:
                        continue

                    # Get last N derivatives and average them
                    avg_x = derivatives_x[-20:]
                    avg_y = derivatives_y[-20:]

                    avg_x = sum(avg_x) / len(avg_x)
                    avg_y = sum(avg_y) / len(avg_y)

                    # Predict N pixels
                    for i in range(7):
                        pos = round(x2 - avg_x * i), round(y2 - avg_y * i)
                        cv2.circle(img, pos, 1, (0, 0, 255), cv2.FILLED)

        predict(pixel_x, pixel_y, [])

cv2.imshow("First derivative", img)
cv2.waitKey(0)

有人知道如何拟合一些曲线/样条线/贝塞尔曲线/.. 并将曲线识别为人类吗？

还有更多的二值化源图片：

谢谢。

Answer 1

众所周知，外推法容易受到初始条件的影响（在这种情况下是进行切割的样条状态。话虽如此，我确实相信一阶和二阶导数足以合理解决这些数据集。

我认为您的外推代码中有两个主要错误：（1）整数像素数据的一阶和二阶导数计算将产生与整数空间的离散性质相关的额外噪声（样条是连续的）。这可以通过将整数精度曲线（使用体素网格滤波器或类似的东西）下采样成双精度曲线来解决。 (2) 您计算二阶导数的方法是沿样条线取两个后续方向向量之间的笛卡尔差。一个更好的策略是取而代之的是后续方向向量之间的角度变化。除了更准确之外，这还允许样条线在自身上回绕并持续传播无损曲线。

我确定碰撞的方法是通过将每条曲线与其他曲线进行外推。在外推的每一步之后，评估了 3 个适应度估计以确定局部最小值（以及随后的曲线组合的最佳适应度）。 (1) 每条曲线端点之间的笛卡尔距离。 (2) 每条曲线端点之间的角距离（180 度偏移被认为是最佳的）。 (3) 累积外推距离。总适应度是这 3 个值的加权和。每个测试候选的外推继续进行，直到总适应度下降。在测试每条曲线的每个组合后，如果最佳适应度低于可配置的阈值，则组合这两条曲线并重复该过程。

每次确定样条合并的理想候选者时，都会对两者进行外推，直到它们在中间相遇。然后使用正弦加权平均这些外插点，以改善样条曲线的曲率/连续性。然后重建三组点（样条 1、桥、样条 2）并重构样条。

代码（c++11，opencv）：

#include <stdio.h>
#include <opencv2/opencv.hpp>
#include <Windows.h>
#include <string>

#define M_PI 3.14159

using namespace std;
using namespace cv;

double sqDist(Point p1, Point p2)
{
    return std::pow((double)p1.x - p2.x,2) + std::pow((double)p1.y - p2.y,2);
}

vector<Point2d> resampleContour(vector<Point> originalContour, int voxelWidth_px = 3)
{
    //quick and dirty voxel filter downsampling to 5px spacing
    vector<Point2d> outputCurve = vector<Point2d>();
    while(originalContour.size()>0)
    {
        int binX = std::floor(originalContour[0].x / voxelWidth_px);
        int binY = std::floor(originalContour[0].y / voxelWidth_px);
        int sumX = originalContour[0].x;
        int sumY = originalContour[0].y;
        int count = 1;
        vector<Point> remainingPoints = vector<Point>();
        for (int i = 1; i < originalContour.size(); i++)
        {
            int tempBinX = std::floor(originalContour[i].x / voxelWidth_px);
            int tempBinY = std::floor(originalContour[i].y / voxelWidth_px);
            if (tempBinX == binX && tempBinY == binY)
            {
                sumX += originalContour[i].x;
                sumY += originalContour[i].y;
                count++;
            }
            else
            {
                remainingPoints.push_back(originalContour[i]);
            }
        }
        originalContour = remainingPoints;
        double aveX = (double)sumX / (double)count;
        double aveY = (double)sumY / (double)count;
        outputCurve.push_back(Point2d(aveX, aveY));
    }

    return outputCurve;
}

Point2d getNN(vector<Point2d> cloud, int targetIndex, int &nnIndex, double &dist)
{
    nnIndex = -1;
    double shortestDist = 0;
    for (int i = 0; i < cloud.size(); i++)
    {
        if (i == targetIndex) { continue; }
        double tempDist = sqDist(cloud[targetIndex], cloud[i]);
        if (nnIndex == -1 || tempDist < shortestDist)
        {
            nnIndex = i;
            shortestDist = tempDist;
        }
    }
    dist = shortestDist;
    return cloud[nnIndex];
}

int getNN(vector<Point2d> cloud, Point2d pt)
{
    int nnIndex = -1;
    double shortestDist = 0;
    for (int i = 0; i < cloud.size(); i++)
    {
        double tempDist = sqDist(pt, cloud[i]);
        if (nnIndex == -1 || tempDist < shortestDist)
        {
            nnIndex = i;
            shortestDist = tempDist;
        }
    }
    return nnIndex;
}

vector<Point2d> getNNs(vector<Point2d> cloud, int targetIndex, int count)
{
    Point2d tPt = cloud[targetIndex];
    sort(cloud.begin(), cloud.end(), [tPt](const Point2d& p1, const Point2d& p2) {
        return sqDist(tPt,p1) < sqDist(tPt, p2);
        });

    vector<Point2d> outCloud = vector<Point2d>();
    for (int i = 1; i <= count && i < cloud.size(); i++) //first point will be same point
    {
        outCloud.push_back(cloud[i]);
    }
    return outCloud;
}

Vec2d normalize(Vec2d inVec)
{
    double length = sqrt(pow(inVec(0), 2) + pow(inVec(1), 2));
    if (length == 0)
    {
        throw;
    }
    inVec = inVec / length;
    return inVec;
}

double angleBetween(Vec2d vec1, Vec2d vec2, bool directionalFlag = false)
{
    vec1 = normalize(vec1);
    vec2 = normalize(vec2);

    double angle = (atan2(vec1(1), vec1(0)) - atan2(vec2(1), vec2(0)));
    if (angle > M_PI) { angle -= 2 * M_PI; }
    else if (angle <= -M_PI) { angle += 2 * M_PI; }

    if (!directionalFlag)
    { angle = abs(angle); }
    return angle;
}

vector<Point> roundPoints(vector<Point2d> cloud)
{
    vector<Point> outCloud = vector<Point>();
    for (int i = 0; i < cloud.size(); i++)
    {
        outCloud.push_back(Point(round(cloud[i].x), round(cloud[i].y)));
    }
    return outCloud;
}

class PointNorm
{
public:
    PointNorm() {}

    //point at p1 with direction p1-p2
    PointNorm(Point2d p1, Point2d p2)
    {
        X = p1.x;
        Y = p1.y;
        dir = Vec2d(p1.x - p2.x, p1.y - p2.y);
        dir = normalize(dir);
    }

    PointNorm(double x,double y,double dx,double dy)
    {
        X = x;
        Y = y;
        dir = Vec2d(dx, dy);
        dir = normalize(dir);
    }
    double X = 0;
    double Y = 0;
    Vec2d dir = Vec2d();

    static void dif(PointNorm pn1, PointNorm pn2, double& distance, double& angle)
    {
        distance = sqrt(pow(pn1.X - pn2.X, 2) + pow(pn1.Y - pn2.Y, 2));
        
        angle = angleBetween(pn1.dir, pn2.dir, false);
    }
};

vector<Point2d> orderCurve(vector<Point2d> cloud)
{
    if (cloud.size() < 3) { return cloud; }

    vector<Point2d> outCloud = vector<Point2d>();

    int currentIndex = cloud.size() / 2;
    double lastDist = -1;
    vector<Point2d> tempCloud = cloud;
    for (int i = 0; i < cloud.size(); i++)
    {
        vector<Point2d> twoNeighbors = getNNs(cloud, i, 5); //technically u can increase this count to increase endpoint confidence
        bool endFlag = true;
        for (int ii = 0; ii < twoNeighbors.size()-1 && endFlag; ii++)
        {
            for (int iii = 0; iii < twoNeighbors.size() - 1 && endFlag; iii++)
            {
                if (ii == iii) { continue; }
                PointNorm pn1 = PointNorm(cloud[i], twoNeighbors[ii]);
                PointNorm pn2 = PointNorm(cloud[i], twoNeighbors[iii]);
                double tempAngle = angleBetween(pn1.dir, pn2.dir);

                if (tempAngle > M_PI / 2)
                { endFlag = false; }
            }
        }

        if (endFlag)
        {
            outCloud.push_back(cloud[i]);
            break;
        }
    }

    tempCloud = cloud;
    currentIndex = getNN(cloud, outCloud[0]);

    while (tempCloud.size() > 1)
    {
        int testIndex = 0;
        double testDist;
        Point2d tempPoint = getNN(tempCloud, currentIndex, testIndex, testDist);
        outCloud.push_back(tempPoint);
        tempCloud.erase(tempCloud.begin() + currentIndex);
        if (testIndex > currentIndex) { testIndex--; }
        currentIndex = testIndex;
    }
    return outCloud;
}

class ModSpline
{
public:

    ModSpline(vector<Point2d> orderedCurve, int dirEstimationOffset, bool _comboCurve = false)
    {
        //totalCurve = orderedCurve;
        totalCurve = vector<Point2d>();
        copy(orderedCurve.begin(), orderedCurve.end(), back_inserter(totalCurve));
        int end = orderedCurve.size() - 1;
        head = PointNorm(orderedCurve[0], orderedCurve[dirEstimationOffset]);//slight offset gives better estimation of direction if last point is noisy
        tail = PointNorm(orderedCurve[end], orderedCurve[end- dirEstimationOffset]);
        comboCurve = _comboCurve;
    }

    void stepExtrapolate_Head(int stepSize_px, int derivativeLookback)
    {
        int tempLookback = derivativeLookback;
        if (tempLookback > totalCurve.size() - 1) { tempLookback = totalCurve.size() - 1; }

        head.X += head.dir(0) * (double)stepSize_px;
        head.Y += head.dir(1) * (double)stepSize_px;
        totalCurve.insert(totalCurve.begin(), 1, Point2d(head.X, head.Y));
        {
            double dirChangeAngle = computeSecondDerivative(0, tempLookback);
            head.dir = normalize(rotateByAngle(head.dir, dirChangeAngle * (double)stepSize_px));
        }
    }
    void stepExtrapolate_Tail(int stepSize_px, int derivativeLookback)
    {
        int tempLookback = derivativeLookback;
        if (tempLookback > totalCurve.size() - 1) { tempLookback = totalCurve.size() - 1; }
        
        tail.X += tail.dir(0) * stepSize_px;
        tail.Y += tail.dir(1) * stepSize_px;
        totalCurve.push_back(Point2d(tail.X, tail.Y));
        {
            double dirChangeAngle = computeSecondDerivative(totalCurve.size() - 1, totalCurve.size() - (1 + tempLookback));
            tail.dir = normalize(rotateByAngle(tail.dir, dirChangeAngle * (double)stepSize_px));
        }
    }

    void stepExtrapolate_Bidirectional(int stepSize_px, int derivativeLookback)
    {
        stepExtrapolate_Head(stepSize_px, derivativeLookback);
        stepExtrapolate_Tail(stepSize_px, derivativeLookback);
    }

    void stepExtrapolate_SingleDirection(int stepSize_px, int derivativeLookback, bool headFlag)
    {
        if (headFlag) { stepExtrapolate_Head(stepSize_px, derivativeLookback); }
        else { stepExtrapolate_Tail(stepSize_px, derivativeLookback); }
    }

    static double estimateSpineDistance(ModSpline spl1, ModSpline spl2,int stepSize_px, int derivativeLookback, double distWeight, double angleWeight, double travelWeight)
    {
        bool dir_1_head, dir_2_head;
        getNearestEndpoints(spl1, spl2, dir_1_head, dir_2_head);

        //todo: run multiple times with adjusted dir and derivatives

        double lastAngle, lastDist, lastComboDist;
        PointNorm::dif(spl1.getEndpointnorm(dir_1_head), spl2.getEndpointnorm(dir_2_head), lastDist, lastAngle);
        lastAngle = abs(M_PI - lastAngle);//looking for opposing vectors;
        lastComboDist = lastAngle * angleWeight + lastDist * distWeight;

        double totalExtrapolationDistance = 0;
        while (true)
        {
            spl1.stepExtrapolate_SingleDirection(stepSize_px, derivativeLookback, dir_1_head);
            spl2.stepExtrapolate_SingleDirection(stepSize_px, derivativeLookback, dir_2_head);
            totalExtrapolationDistance += (stepSize_px + stepSize_px);
            double tempAngle, tempDist, tempComboDist;
            PointNorm::dif(spl1.getEndpointnorm(dir_1_head), spl2.getEndpointnorm(dir_2_head), tempDist, tempAngle);
            tempAngle = abs(M_PI - tempAngle);//looking for opposing vectors;

            tempComboDist = tempAngle * angleWeight + tempDist * distWeight +totalExtrapolationDistance* travelWeight;
            if (tempComboDist > lastComboDist) { break; }
            else
            {
                lastAngle = tempAngle;
                lastDist = tempDist;
                lastComboDist = tempComboDist;
            }
        }
        return lastComboDist;
    }

    static ModSpline mergeSplines(ModSpline spl1, ModSpline spl2, int stepSize_px, int derivativeLookback, int dirEstimationOffset, vector<vector<Point2d>> &debugClouds)
    {
        bool dir_1_head, dir_2_head;
        getNearestEndpoints(spl1, spl2, dir_1_head, dir_2_head);

        double lastAngle, lastDist, lastComboDist;
        int extrapolationCount = 0;
        PointNorm::dif(spl1.getEndpointnorm(dir_1_head), spl2.getEndpointnorm(dir_2_head),lastDist, lastAngle);
        while (true)
        {
            extrapolationCount += 1;
            spl1.stepExtrapolate_SingleDirection(stepSize_px, derivativeLookback, dir_1_head);
            spl2.stepExtrapolate_SingleDirection(stepSize_px, derivativeLookback, dir_2_head);
            double tempAngle, tempDist, tempComboDist;
            PointNorm::dif(spl1.getEndpointnorm(dir_1_head), spl2.getEndpointnorm(dir_2_head), tempDist, tempAngle);

            if (tempDist > lastDist) { break; }
            else { lastDist = tempDist; }
        }

        vector<Point2d> outCloud = vector<Point2d>();

        vector<Point2d> spline1Cloud = spl1.getTotalCurve();
        if (dir_1_head) { reverse(spline1Cloud.begin(), spline1Cloud.end()); }
        vector<Point2d> spline2Cloud = spl2.getTotalCurve();
        if (!dir_2_head) { reverse(spline2Cloud.begin(), spline2Cloud.end()); }

        //spline 1
        {
            vector<Point2d> debug = vector<Point2d>();
            for (int i = 0; i < spline1Cloud.size() - extrapolationCount; i++)
            {
                outCloud.push_back(spline1Cloud[i]);
                debug.push_back(spline1Cloud[i]);
            }
            debugClouds.push_back(debug);
        }
        //fused region
        {
            vector<Point2d> debug = vector<Point2d>();
            double theta = 0;
            double theta_Step = (M_PI / 2.0) / (double)extrapolationCount;
            for (int i = 1; i < extrapolationCount-1; i++)
            {
                int invI = extrapolationCount - i;
                Point2d p1 = spline1Cloud[spline1Cloud.size() - (1 + invI)];
                Point2d p2 = spline2Cloud[i];

                double alpha = sin(theta);
                theta += theta_Step;

                //sinusoidal interpolation
                Point2d fusedPt = Point2d(alpha*p2.x+(1.0-alpha)*p1.x, alpha * p2.y + (1.0 - alpha) * p1.y);

                outCloud.push_back(fusedPt);
                debug.push_back(fusedPt);
            }
            debugClouds.push_back(debug);
        }
        //spline 2
        {
            vector<Point2d> debug = vector<Point2d>();
            for (int i = extrapolationCount; i < spline2Cloud.size(); i++)
            {
                outCloud.push_back(spline2Cloud[i]);
                debug.push_back(spline2Cloud[i]);
            }
            debugClouds.push_back(debug);
        }
        return ModSpline(outCloud,dirEstimationOffset,true);
    }

    vector<Point2d> getTotalCurve()
    {
        return totalCurve;
    }

    int getPtCount() { return totalCurve.size(); }

    vector<PointNorm> getEndpoints()
    {
        vector<PointNorm> endPoints = vector<PointNorm>();
        endPoints.push_back(head);
        endPoints.push_back(tail);
        return endPoints;
    }

    bool isComboCurve() { return comboCurve; }

    Point2d getEndpoint(bool headFlag)
    {
        if (headFlag) { return Point2d(head.X, head.Y); }
        else { return Point2d(tail.X, tail.Y); }
    }

    PointNorm getEndpointnorm(bool headFlag)
    {
        if (headFlag) { return head; }
        else { return tail; }
    }

    static void getNearestEndpoints(ModSpline spl1, ModSpline spl2, bool& headFlag1, bool& headFlag2)
    {
        double dist1 = sqDist(Point2d(spl1.head.X, spl1.head.Y), Point2d(spl2.head.X, spl2.head.Y));
        double dist2 = sqDist(Point2d(spl1.head.X, spl1.head.Y), Point2d(spl2.tail.X, spl2.tail.Y));
        double dist3 = sqDist(Point2d(spl1.tail.X, spl1.tail.Y), Point2d(spl2.head.X, spl2.head.Y));
        double dist4 = sqDist(Point2d(spl1.tail.X, spl1.tail.Y), Point2d(spl2.tail.X, spl2.tail.Y));

        double minDist = min(dist1, min(dist2, min(dist3, dist4)));
        if (minDist == dist1) { headFlag1 = true; headFlag2 = true; }
        else if (minDist == dist2) { headFlag1 = true; headFlag2 = false; }
        else if (minDist == dist3) { headFlag1 = false; headFlag2 = true; }
        else if (minDist == dist4) { headFlag1 = false; headFlag2 = false; }
    }

private:
    vector<Point2d> totalCurve;
    PointNorm head;
    PointNorm tail;
    bool comboCurve = false;

    double computeSecondDerivative(int startIndex, int endIndex)
    {
        int increment = 1;
        if (endIndex < startIndex) { increment = -1; }

        double totAngle = 0;
        double totalDistance = 0;
        int count = 0;
        for (int i = startIndex; i + 2 * increment != endIndex; i += increment)
        {
            count++;

            Point2d p1 = totalCurve[i];
            Point2d p2 = totalCurve[i + increment];
            Point2d p3 = totalCurve[i + 2 * increment];

            Vec2d v1 = Vec2d(p1.x - p2.x, p1.y - p2.y);
            Vec2d v2 = Vec2d(p2.x - p3.x, p2.y - p3.y);

            double tempAngle = angleBetween(v1, v2, true);
            double aveDist = (sqrt(sqDist(p1, p2)) + sqrt(sqDist(p2, p3))) / 2.0;
            totalDistance += aveDist;
            totAngle += tempAngle;
        }

        totAngle = totAngle / totalDistance;
        return totAngle;
    }

    static Vec2d rotateByAngle(Vec2d inVec, double angle)
    {
        Vec2d outVec = Vec2d();
        outVec(0) = inVec(0)*cos(angle) - inVec(1)*sin(angle);
        outVec(1) = inVec(0)*sin(angle) + inVec(1)*cos(angle);
        return outVec;
    }
};


vector<Scalar> colorWheel = { 
    Scalar(255,0,0),
    Scalar(0,255,0),
    Scalar(0,0,255),
    Scalar(255,255,0),
    Scalar(255,0,255),
    Scalar(0,255,255) };
int colorWheelIndex = 0;
void rotateColorWheel()
{
    colorWheelIndex++;
    if (colorWheelIndex >= colorWheel.size()) { colorWheelIndex = 0; }
}


int main(int argc, char** argv)
{

    //control downsampling and extrapolation (several these could be static members of modspline instead)
    const int stepSize_px = 1;//how far each extrapolation step travels
    const int derivativeLookback = 15;//how far back in curve is considered in determining 2nd derivative
    const int dirEstimationOffset = 2;//min 1.  determines how much deviations at ends of curves effect extrapolation (lower equals more impact)
    const int voxelWidth = 5; //voxel dimension for averaging intial pixels into point doubles

    //control spline similarity calculation (each of these contributes to a weighted sum of "spline distance")
    const double distWeighting = 1.0; //how much distance between two splines after optimal extrapolation
    const double angleWeighting = 10.0; //how much angle between two splines after optimal extrapolation
    const double travelWeighting = 0.05; //how far splines have to travel to connect
    const double maxTotalSplineError = 15; //unitless weighted combination of "spline distance"

    bool debugFlag = true;// flag to enable or disable debug visualizers


    std::string path = "C:/Local Software/voyDICOM/resources/images/SparseCurves6.png";

    Mat originalImg = cv::imread(path, cv::IMREAD_GRAYSCALE);
    if (debugFlag) { cv::imshow("orignal", originalImg); }

    Mat debugImg = cv::imread(path, cv::IMREAD_COLOR);

    Mat bwImg;
    threshold(originalImg, bwImg, 150, 255, THRESH_BINARY);

    vector<vector<Point> > contours;
    findContours(bwImg, contours, RETR_LIST, CHAIN_APPROX_NONE);

    vector<vector<Point2d>> dCurves = vector<vector<Point2d>>();

    Mat initialCurves = debugImg.clone();
    for (int i = 0; i < contours.size(); i++)
    {
        vector<Point2d> tempDCurve = resampleContour(contours[i], voxelWidth);
        if (tempDCurve.size() < 7) { continue; }

        tempDCurve = orderCurve(tempDCurve);
        dCurves.push_back(tempDCurve);

        vector<Point> tempCloud = roundPoints(tempDCurve);
        cv::polylines(initialCurves, { tempCloud }, false, colorWheel[colorWheelIndex], 2);
        circle(initialCurves, tempCloud[0], 5, Scalar(0, 255, 0), 1);
        circle(initialCurves, tempCloud[tempCloud.size() - 1], 5, Scalar(0, 0, 255), 1);
        rotateColorWheel();
    }
    if (debugFlag) { cv::imshow("initialCurves", initialCurves); }
        
    

    vector<ModSpline> travCurves = vector<ModSpline>();
    vector<ModSpline> tempCurves = vector<ModSpline>();
    for (int i = 0; i < dCurves.size(); i++)
    {
        travCurves.push_back(ModSpline(dCurves[i], dirEstimationOffset));
        tempCurves.push_back(ModSpline(dCurves[i], dirEstimationOffset));
    }


    if (debugFlag) 
    {
        for (int i = 0; i < 25; i++)
        {
            colorWheelIndex = 0;
            Mat extrapViewer = debugImg.clone();
            for (int ii = 0; ii < tempCurves.size(); ii++)
            {
                //if (!(ii == 7 || ii == 13)) { continue; }
                tempCurves[ii].stepExtrapolate_Bidirectional(stepSize_px,derivativeLookback);
                vector<Point2d> tempCloud = tempCurves[ii].getTotalCurve();
                for (int iii = 0; iii < tempCloud.size(); iii++)
                {
                    cv::circle(extrapViewer, tempCloud[iii], 1, colorWheel[colorWheelIndex], 1);
                }
                rotateColorWheel();
            }
            cv::imshow("extrapolation debug", extrapViewer);
            cv::waitKey(100);
        }
    }
    
    
    int fusionCounter = 1;
    while (true)
    {
        vector<tuple<int, int>> validMerges = vector<tuple<int, int>>();
        for (int i = 0; i < travCurves.size(); i++)
        {
            for (int ii = 0; ii < travCurves.size(); ii++)
            {
                if (i == ii) { continue; }
                double tempComboDist = ModSpline::estimateSpineDistance(
                    travCurves[i], 
                    travCurves[ii], 
                    stepSize_px,
                    derivativeLookback,
                    distWeighting,
                    angleWeighting,
                    travelWeighting);
                if (tempComboDist < maxTotalSplineError)
                {
                    validMerges.push_back(tuple<int,int>(i,ii));
                }
            }
        }

        if (validMerges.size()>0)
        {
            vector<int> splineSizes = vector<int>();
            for (int i = 0; i < travCurves.size(); i++)
            {
                splineSizes.push_back(travCurves[i].getPtCount());
            }

            //sort valid merges by combined spline size (favor larger spline merges before smaller ones)
            sort(validMerges.begin(), validMerges.end(), [splineSizes](const tuple<int, int>& spl1, const tuple<int, int>& spl2) {
                return  splineSizes[get<0>(spl1)] + splineSizes[get<1>(spl1)] > splineSizes[get<0>(spl2)] + splineSizes[get<1>(spl2)];
                });

            int splineIndex1 = get<0>(validMerges[0]);
            int splineIndex2 = get<1>(validMerges[0]);

            vector<vector<Point2d>> debugClouds = vector<vector<Point2d>>();
            ModSpline newCurve = ModSpline::mergeSplines(
                travCurves[splineIndex1], 
                travCurves[splineIndex2], 
                stepSize_px, 
                derivativeLookback,
                dirEstimationOffset,
                debugClouds);

            travCurves.erase(travCurves.begin() + max(splineIndex1, splineIndex2));
            travCurves.erase(travCurves.begin() + min(splineIndex1, splineIndex2));
            travCurves.push_back(newCurve);

            if (debugFlag)
            {
                Mat debugFusions = debugImg.clone();
                cv::polylines(debugFusions, { roundPoints(debugClouds[0]) }, false, Scalar(255, 0, 0), 2);
                cv::polylines(debugFusions, { roundPoints(debugClouds[1]) }, false, Scalar(0, 255, 0), 1);
                cv::polylines(debugFusions, { roundPoints(debugClouds[2]) }, false, Scalar(0, 0, 255), 2);
                cv::imshow("debugFusion"+std::to_string(fusionCounter++), debugFusions);
                cv::waitKey(500);
            }
        }
        else
        {
            break;
        }
    }

    Mat finalCurves = debugImg.clone();
    colorWheelIndex = 0;
    for (int i = 0; i < travCurves.size(); i++)
    {
        if (!travCurves[i].isComboCurve()) { continue; }
        cv::polylines(finalCurves, { roundPoints(travCurves[i].getTotalCurve()) }, false, colorWheel[colorWheelIndex], 2);
        rotateColorWheel();
    }
    cv::imshow("final curves", finalCurves);
    cv::imshow("initialCurves", initialCurves);
    cv::imshow("original", originalImg);
    
    cv::waitKey(0);
}

---

其他数据集（所有 3 个数据集的参数相同）：

最后说明：在主函数顶部使用变量可以完全控制样条碰撞检测的特异性。使其更健壮的下一步将是建立一个大型数据集并根据需要修改参数，以便为每个数据集提供最佳解决方案。然后，分析所有这些单独的设置配置，您应该能够为设置建立中心点和置信范围，以便单个配置适用于最大数量的测试图像。