缺失数据的 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(ModSplin