I have split contours in a image and I want to connect them somehow. I was considering some curve fitting or curve tracing - but I do not know, how can I implement it.
I have exactly 1 px wide split contours:
I know, where are ends of all contours:
As a human I know, where are the curves and what contours are part of each unique curve:
But I do not know, how can I connect these split contours programmatically.
I tried using first and second derivative to predict, but it was not good enough:
Anybody knows how to fit some curve/spline/Bezier/.. to it and recognize the curves as a human?
Thanks.
CodePudding user response：
So you have sets of curves represented by 2d points. You seem to have extrapolated with 1st and 2nd derivative reasonably accurately. In order to determine collisions, I think a solid strategy would be to consider "distance" in 4 dimensional space: (x,y,vx,vy) so both position and current direction. Now iterate through all curves, stepping them one step forward along their respective extrapolation equations (on both ends) and then compare all sets of active endpoints for two properties: cartesian distance, and angular separation. You could make it into a single distance function, but I feel like keeping them separate will be easier since angles will be nonlinear (that being said if it is easier to combine them into a single weighted value, u do u). Any time two end points are sufficiently close and acceptably similar in directional vector (keep in mind that u should test angleBetween and 180-angleBetween since connective curves will likely be extrapolating in opposite directions (you may actually only want the ones opposing each-other)), you can connect those two curves and delete the colliding endpoints, creating a new curve with the remaining two endpoints of the combination. Also reaching the edge of the image should kill an endpoint.
An additional consideration depending on your resources is that you could add intentional variation into the extrapolations (ie create duplicate curves that have slightly adjusted first and second derivatives). So each curve would almost extrapolate in a fan pattern. This will require a fair bit more memory and you will need to handle priority of better collisions overwriting worse collisions, but would definitely yield a stronger solution (especially in the presence of sparse data).
One final note is that you can do these two strategies sequentially, so first start with just stepping and looking for collisions, and pulling out any curves that find good collisions. Then, take the remaining curves and start over (remove the original extrapolation), but apply a shift to the first and or second derivative (u could do it exhaustively or with some random variation depending on your preference). Keep in mind that in this situation, you will only need to apply that perturbation for the first step of the extrapolation and then let its effects carry forward through your regular extrapolation.
Edit: added code to demonstrate the core concept. First thing of note, sorry for how trash the code is. Second thing, there are a ton of "magic variables" in this code (which is kinda inherent to this process). My suggestion is to run it multiple times with different settings iteratively and take result with highest overall continuity. (I will probably have a go at trying something like this after I get some sleep). One last note is that the extrapolations are very sensitive to initial conditions and as I described in the methodology before, adding a little bit of intentional variation to those initial settings (direction and curvature) may yield better results which again can be collected from a large numerical set based on best image continuity. Sorry again for the damage to your eyes.
#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<Point> resampleContour(vector<Point> originalContour, int voxelWidth_px = 3)
{
//quick and dirty voxel filter downsampling to 5px spacing
vector<Point> outputCurve = vector<Point>();
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;
int aveX = std::round((double)sumX / (double)count);
int aveY = std::round((double)sumY / (double)count);
outputCurve.push_back(Point(aveX, aveY));
}
return outputCurve;
}
Point getNN(vector<Point> 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<Point> cloud, Point 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<Point> orderCurve(vector<Point> cloud)
{
if (cloud.size() < 3) { return cloud; }
vector<Point> outCloud = vector<Point>();
int currentIndex = cloud.size() / 2;
double lastDist = -1;
vector<Point> tempCloud = cloud;
while (true)
{
if (tempCloud.size() == 1)
{
outCloud.push_back(tempCloud[0]);
break;
}
int testIndex = 0;
double testDist;
getNN(tempCloud, currentIndex, testIndex, testDist);
if (lastDist == -1)
{
lastDist = testDist;
tempCloud.erase(tempCloud.begin() currentIndex);
if (testIndex > currentIndex) { testIndex--; }
currentIndex = testIndex;
continue;
}
else
{
lastDist = testDist;
tempCloud.erase(tempCloud.begin() currentIndex);
if (testIndex > currentIndex) { testIndex--; }
currentIndex = testIndex;
}
}
tempCloud = cloud;
currentIndex = getNN(cloud, outCloud[0]);
while (tempCloud.size() > 1)
{
int testIndex = 0;
double testDist;
Point tempPoint = getNN(tempCloud, currentIndex, testIndex, testDist);
outCloud.push_back(tempPoint);
tempCloud.erase(tempCloud.begin() currentIndex);
if (testIndex > currentIndex) { testIndex--; }
currentIndex = testIndex;
}
return outCloud;
}
Vec2d normalize(Vec2d inVec)
{
double length = sqrt(pow(inVec(0), 2) pow(inVec(1), 2));
inVec = inVec / length;
return inVec;
}
class PointNorm
{
public:
PointNorm() {}
//point at p1 with direction p1-p2
PointNorm(Point p1, Point 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 = abs((atan2(pn2.dir(0), pn2.dir(1)) - atan2(pn1.dir(0),pn1.dir(1))));
}
};
class TravelingCurve
{
public:
TravelingCurve(vector<Point> orderedCurve, bool _comboFlag = false)
{
totalCurve = orderedCurve;
int end = orderedCurve.size() - 1;
head = PointNorm(orderedCurve[0], orderedCurve[2]);//slight offset gives better estimation of direction if last point is noisy
tail = PointNorm(orderedCurve[end], orderedCurve[end-2]);
comboFlag = _comboFlag;
}
void stepExtrapolate(int stepSize_px = 3, int derivativeLookback = 10)
{
if (derivativeLookback > totalCurve.size() - 1) { derivativeLookback = totalCurve.size() - 1; }
head.X = head.dir(0) * stepSize_px;
head.Y = head.dir(1) * stepSize_px;
totalCurve.insert(totalCurve.begin(), 1,Point(round(head.X),round(head.Y)));
{
double dirChangeAngle = computeSecondDerivative(0, derivativeLookback);
head.dir = normalize(rotateByAngle(head.dir, dirChangeAngle));
}
tail.X = tail.dir(0) * stepSize_px;
tail.Y = tail.dir(1) * stepSize_px;
totalCurve.push_back(Point(round(tail.X), round(tail.Y)));
{
double dirChangeAngle = computeSecondDerivative(totalCurve.size()-1, totalCurve.size() - (1 derivativeLookback));
tail.dir = normalize(rotateByAngle(tail.dir, dirChangeAngle));
}
}
vector<Point> getTotalCurve()
{
return totalCurve;
}
vector<PointNorm> getEndpoints()
{
vector<PointNorm> endPoints = vector<PointNorm>();
endPoints.push_back(head);
endPoints.push_back(tail);
return endPoints;
}
static TravelingCurve combineCurves(TravelingCurve c1, TravelingCurve c2)
{
//this could def be done cleaner...
vector<Point> cloud1 = c1.getTotalCurve();
vector<Point> cloud2 = c2.getTotalCurve();
double dist1 = sqDist(Point(c1.head.X, c1.head.Y), Point(c2.head.X, c2.head.Y));
double dist2 = sqDist(Point(c1.head.X, c1.head.Y), Point(c2.tail.X, c2.tail.Y));
double dist3 = sqDist(Point(c1.tail.X, c1.tail.Y), Point(c2.head.X, c2.head.Y));
double dist4 = sqDist(Point(c1.tail.X, c1.tail.Y), Point(c2.tail.X, c2.tail.Y));
double minDist = min(dist1, min(dist2, min(dist3, dist4)));
if (minDist == dist1)
{
reverse(cloud1.begin(), cloud1.end());
cloud1.insert(cloud1.end(), cloud2.begin(), cloud2.end());
return TravelingCurve(cloud1, true);
}
else if (minDist == dist2)
{
cloud2.insert(cloud2.end(), cloud1.begin(), cloud1.end());
return TravelingCurve(cloud2, true);
}
else if (minDist == dist3)
{
cloud1.insert(cloud1.end(), cloud2.begin(), cloud2.end());
return TravelingCurve(cloud1, true);
}
else if (minDist == dist4)
{
reverse(cloud2.begin(), cloud2.end());
cloud1.insert(cloud1.end(), cloud2.begin(), cloud2.end());
return TravelingCurve(cloud1,true);
}
}
bool isComboCurve() { return comboFlag; }
private:
vector<Point> totalCurve;
PointNorm head;
PointNorm tail;
bool comboFlag;
double computeSecondDerivative(int startIndex, int endIndex)
{
int increment = 1;
if (endIndex < startIndex) { increment = -1; }
double aveAngle = 0;
int count = 0;
for (int i = startIndex; i 2 * increment != endIndex; i = increment)
{
count ;
Point p1 = totalCurve[i];
Point p2 = totalCurve[i increment];
Point 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);
v1 = normalize(v1);
v2 = normalize(v2);
double tempAngle = (atan2(v1(1), v1(0)) - atan2(v2(1), v2(0)));
if (tempAngle > M_PI) { tempAngle -= 2 * M_PI; }
else if (tempAngle <= -M_PI) { tempAngle = 2 * M_PI; }
aveAngle = tempAngle;
}
aveAngle = aveAngle / (double)count;
return aveAngle;
}
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)
{
bool debugFlag = true;
std::string path = "C:/Local Software/voyDICOM/resources/images/SparseCurves.png";
Mat originalImg = cv::imread(path, cv::IMREAD_GRAYSCALE);
if (debugFlag) { 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<Point>> curves = vector<vector<Point>>();
for (int i = 0; i < contours.size(); i )
{
if (contours[i].size() >= 30)
{
curves.push_back(contours[i]);
}
}
Mat initialCurves = debugImg.clone();
for (int i = 0; i < curves.size(); i )
{
curves[i] = resampleContour(curves[i], 5);
curves[i] = orderCurve(curves[i]);
vector<vector<Point>> singleCurve = vector<vector<Point>>();
singleCurve.push_back(curves[i]);
cv::polylines(initialCurves, singleCurve, false, colorWheel[colorWheelIndex], 2);
circle(initialCurves, curves[i][0], 5, Scalar(0, 255, 0), 1);
circle(initialCurves, curves[i][curves[i].size() - 1], 5, Scalar(0, 0, 255), 1);
rotateColorWheel();
}
imshow("initialCurves", initialCurves);
vector<TravelingCurve> travCurves = vector<TravelingCurve>();
for (int i = 0; i < curves.size(); i )
{
travCurves.push_back(TravelingCurve(curves[i]));
}
int iterations = 0;
while(iterations<20)
{
iterations ;
//check for any curve collisions
bool collisionFound = true;
while (collisionFound)
{
collisionFound = false;
vector<PointNorm> allEndPoints = vector<PointNorm>();
for (int i = 0; i < travCurves.size(); i )
{
vector<PointNorm> tempEndPoints = travCurves[i].getEndpoints();
allEndPoints.insert(allEndPoints.end(), tempEndPoints.begin(), tempEndPoints.end());
}
for (int i = 0; i < allEndPoints.size(); i )
{
for (int ii = 0; ii < allEndPoints.size(); ii )
{
if (i / 2 == ii / 2) { continue; }
double dist, angle;
PointNorm::dif(allEndPoints[i], allEndPoints[ii], dist, angle);
if (dist < 6 && abs(M_PI - angle) < M_PI * 0.2)
{
std::cout << "collision" << std::endl;
int index1 = i / 2;
int index2 = ii / 2;
TravelingCurve newCurve = TravelingCurve::combineCurves(travCurves[index1], travCurves[index2]);
travCurves.erase(travCurves.begin() max(index1, index2));
travCurves.erase(travCurves.begin() min(index1, index2));
travCurves.push_back(newCurve);
collisionFound = true;
}
if (collisionFound) { break; }
}
if (collisionFound) { break; }
}
}
//extrapolate all curve endpoints
Mat growingCurves = debugImg.clone();
for (int i = 0; i < travCurves.size(); i )
{
travCurves[i].stepExtrapolate(3);
vector<vector<Point>> singleCurve = vector<vector<Point>>();
singleCurve.push_back(travCurves[i].getTotalCurve());
cv::polylines(growingCurves, singleCurve, false, colorWheel[colorWheelIndex], 2);
rotateColorWheel();
}
imshow("debug curves", growingCurves);
waitKey(0);
}
Mat finalCurves = debugImg.clone();
for (int i = 0; i < travCurves.size(); i )
{
if (!travCurves[i].isComboCurve()) { continue; }
vector<vector<Point>> singleCurve = vector<vector<Point>>();
singleCurve.push_back(travCurves[i].getTotalCurve());
cv::polylines(finalCurves, singleCurve, false, colorWheel[colorWheelIndex], 2);
rotateColorWheel();
}
imshow("final curves", finalCurves);
imshow("initial Curves", initialCurves);
imshow("original Image", originalImg);
waitKey(0);
}
If you have any questions about the code, let me know (this problem is pretty fun).