Как подогнать круг к набору точек с ограниченным радиусом?

Почему бы не использовать HoughCircles()? Это позволяет вам указать оба параметра. Пример взят из связанных документов, изменен для добавления параметров минимального и максимального радиуса:

#include <opencv2/imgproc.hpp>
#include <opencv2/highgui.hpp>
#include <math.h>

using namespace cv;

int main(int argc, char** argv)
{
Mat img, gray;
img = imread('ring.png')
cvtColor(img, gray, COLOR_BGR2GRAY);
// smooth it, otherwise a lot of false circles may be detected
medianBlur(gray, gray, 7);
vector<Vec3f> circles;
int minRadius = 90
int maxRadius = 100
HoughCircles(gray, circles, HOUGH_GRADIENT,
1, 100, 20, 20, minRadius, maxRadius);
for( size_t i = 0; i < circles.size(); i++ )
{
Point center(cvRound(circles[i][0]), cvRound(circles[i][1]));
int radius = cvRound(circles[i][2]);
// draw the circle center
circle( img, center, 3, Scalar(0,255,0), -1, 8, 0 );
// draw the circle outline
circle( img, center, radius, Scalar(0,0,255), 3, 8, 0 );
}
namedWindow( "circles", 1 );
imshow( "circles", img );
return 0;
}

Пример выполнения (обратите внимание, я сделал это с Python, но с теми же параметрами):

Изображение украдено из аспирант в КМУ.

Используйте уравнение координат из

https://de.wikipedia.org/wiki/Kreis

     G: a*x + b*y + c + y^2 + x^2 = 0
with
a: -2*xm
b: -2*ym
c: ym^2+xm^2-r^2

(Это точное уравнение не упоминается в английской Википедии, но упоминается базовая форма.)

Чтобы получить линейную регрессию, дифференцируйте G ^ 2 для a / b / c

Тогда вы получите в виде матрицы

           A        B         C      R        R
G1:  +2*a*x^2 +2*b*x*y  +2*c*x +2*x^3    +2*x*y^2  = 0
G2:  +2*a*x*y +2*b*y^2  +2*c*y +2*y^3    +2*x^2*y  = 0
G3:  +2*a*x   +2*b*y    +2*c   +2*y^2    +2*x^2    = 0

Это производные уравнения вы можете решить с вашими точками для a / b / c. Из a / b / c вы можете заменить, чтобы получить xm / ym / r. Чтобы отклонить определенный радиус, просто проверьте его пределы.

Пример:

РЕДАКТИРОВАТЬ: Источник

#include <opencv2/opencv.hpp>

#define _USE_MATH_DEFINES
#include <math.h>

#include <vector>
#include <random>

void fit_circle(const std::vector<cv::Point2d> &pnts, cv::Point2d &centre,
double &radius)
{
/*
A        B         C      R        R
G1:  +2*a*x^2 +2*b*x*y  +2*c*x +2*x^3    +2*x*y^2  = 0
G2:  +2*a*x*y +2*b*y^2  +2*c*y +2*y^3    +2*x^2*y  = 0
G3:  +2*a*x   +2*b*y    +2*c   +2*y^2    +2*x^2    = 0
*/

static const int rows = 3;
static const int cols = 3;
cv::Mat LHS(rows, cols, CV_64F, 0.0);
cv::Mat RHS(rows, 1, CV_64F, 0.0);
cv::Mat solution(rows, 1, CV_64F, 0.0);

if (pnts.size() < 3)
{
throw std::runtime_error("To less points");
}

for (int i = 0; i < static_cast<int>(pnts.size()); i++)
{
double x1 = pnts[i].x;
double x2 = std::pow(pnts[i].x, 2);
double x3 = std::pow(pnts[i].x, 3);
double y1 = pnts[i].y;
double y2 = std::pow(pnts[i].y, 2);
double y3 = std::pow(pnts[i].y, 3);

// col 0 = A / col 1 = B / col 2 = C
// Row 0 = G1
LHS.at<double>(0, 0) += 2 * x2;
LHS.at<double>(0, 1) += 2 * x1 * y1;
LHS.at<double>(0, 2) += 2 * x1;

RHS.at<double>(0, 0) -= 2 * x3 + 2 * x1 * y2;

// Row 1 = G2
LHS.at<double>(1, 0) += 2 * x1 * y1;
LHS.at<double>(1, 1) += 2 * y2;
LHS.at<double>(1, 2) += 2 * y1;

RHS.at<double>(1, 0) -= 2 * y3 + 2 * x2 * y1;

// Row 2 = G3
LHS.at<double>(2, 0) += 2 * x1;
LHS.at<double>(2, 1) += 2 * y1;
LHS.at<double>(2, 2) += 2;

RHS.at<double>(2, 0) -= 2 * y2 + 2 * x2;
}

cv::solve(LHS, RHS, solution);

std::vector<double> abc{solution.at<double>(0, 0),
solution.at<double>(1, 0),
solution.at<double>(2, 0)};

centre.x = abc[0] / -2.0;
centre.y = abc[1] / -2.0;
radius = std::sqrt(
std::abs(std::pow(centre.x, 2) + std::pow(centre.y, 2) - abc[2]));
}

int main(int argc, const char *argv[])
{
try
{
std::random_device rd;
std::mt19937 gen(rd());
std::uniform_real_distribution<> dis(-0.5, 0.5);

// Generate reandom points
double radius_in = 25;
double xm_in = 10;
double ym_in = 20;

std::vector<cv::Point2d> pnts;
{
for (double ang = 0; ang <= 90; ang += 10)
{
cv::Point2d p(
cos(ang / 180.0 * M_PI) * radius_in + xm_in + dis(gen),
sin(ang / 180.0 * M_PI) * radius_in + ym_in + dis(gen));
pnts.push_back(p);
}
}

cv::Point2d c_out;
double radius_out = 0;
fit_circle(pnts, c_out, radius_out);

double dxm = c_out.x - xm_in;
double dym = c_out.y - ym_in;
double dradius = radius_out - radius_in;

std::cout << "Deltas: " << dxm << " " << dym << " " << dradius
<< std::endl;
return EXIT_SUCCESS;
}
catch (const std::exception &ex)
{
std::cerr << ex.what() << std::endl;
}

return EXIT_FAILURE;
}

Это выводы
Дельты: 0,180365 0,190016 -0,231563

РЕДАКТИРОВАТЬ: Добавить источник для создания статистики и статистики об алгоритме

#include <opencv2/opencv.hpp>

#define _USE_MATH_DEFINES
#include <math.h>

#include <vector>
#include <random>

void fit_circle(const std::vector<cv::Point2d> &pnts, cv::Point2d &centre,
double &radius)
{
/*
A        B         C      R        R
G1:  +2*a*x^2 +2*b*x*y  +2*c*x +2*x^3    +2*x*y^2  = 0
G2:  +2*a*x*y +2*b*y^2  +2*c*y +2*y^3    +2*x^2*y  = 0
G3:  +2*a*x   +2*b*y    +2*c   +2*y^2    +2*x^2    = 0
*/

static const int rows = 3;
static const int cols = 3;
cv::Mat LHS(rows, cols, CV_64F, 0.0);
cv::Mat RHS(rows, 1, CV_64F, 0.0);
cv::Mat solution(rows, 1, CV_64F, 0.0);

if (pnts.size() < 3)
{
throw std::runtime_error("To less points");
}

for (int i = 0; i < static_cast<int>(pnts.size()); i++)
{
double x1 = pnts[i].x;
double x2 = std::pow(pnts[i].x, 2);
double x3 = std::pow(pnts[i].x, 3);
double y1 = pnts[i].y;
double y2 = std::pow(pnts[i].y, 2);
double y3 = std::pow(pnts[i].y, 3);

// col 0 = A / col 1 = B / col 2 = C
// Row 0 = G1
LHS.at<double>(0, 0) += 2 * x2;
LHS.at<double>(0, 1) += 2 * x1 * y1;
LHS.at<double>(0, 2) += 2 * x1;

RHS.at<double>(0, 0) -= 2 * x3 + 2 * x1 * y2;

// Row 1 = G2
LHS.at<double>(1, 0) += 2 * x1 * y1;
LHS.at<double>(1, 1) += 2 * y2;
LHS.at<double>(1, 2) += 2 * y1;

RHS.at<double>(1, 0) -= 2 * y3 + 2 * x2 * y1;

// Row 2 = G3
LHS.at<double>(2, 0) += 2 * x1;
LHS.at<double>(2, 1) += 2 * y1;
LHS.at<double>(2, 2) += 2;

RHS.at<double>(2, 0) -= 2 * y2 + 2 * x2;
}

cv::solve(LHS, RHS, solution);

std::vector<double> abc{solution.at<double>(0, 0),
solution.at<double>(1, 0),
solution.at<double>(2, 0)};

centre.x = abc[0] / -2.0;
centre.y = abc[1] / -2.0;
radius = std::sqrt(
std::abs(std::pow(centre.x, 2) + std::pow(centre.y, 2) - abc[2]));
}

int main(int argc, const char *argv[])
{
int count = 100;
if (argc == 2)
{
count = std::atoi(argv[1]);
}

try
{
std::random_device rd;
std::mt19937 gen(rd());
std::uniform_real_distribution<> jitter(-0.5, 0.5);
std::uniform_real_distribution<> pos(-100, 100);

std::uniform_real_distribution<> w_start(0, 360);
std::uniform_real_distribution<> w_length(10, 180);
std::uniform_real_distribution<> rad(10, 180);

std::cout << "start\tlen\tdelta xm\tdelta ym\tdelta radius" << std::endl;

for (int i = 0; i < count; ++i)
{
// Generate reandom points
double radius_in = rad(gen);
double xm_in = pos(gen);
double ym_in = pos(gen);

double start = w_start(gen);
double len = w_length(gen);

std::vector<cv::Point2d> pnts;
{
for (double ang = start; ang <= start + len; ang += 1)
{
cv::Point2d p(
cos(ang / 180.0 * M_PI) * radius_in + xm_in + jitter(gen),
sin(ang / 180.0 * M_PI) * radius_in + ym_in + jitter(gen));
pnts.push_back(p);
}
}

cv::Point2d c_out;
double radius_out = 0;
fit_circle(pnts, c_out, radius_out);

double dxm = c_out.x - xm_in;
double dym = c_out.y - ym_in;
double dradius = radius_out - radius_in;

std::cout << start << "\t" << len << "\t" << dxm << "\t" << dym << "\t" << dradius
<< std::endl;
}

return EXIT_SUCCESS;
}
catch (const std::exception &ex)
{
std::cerr << ex.what() << std::endl;
}

return EXIT_FAILURE;
}

100000 кругов

Участок сформирован:

set multiplot layout 3,1 rowsfirst
set yrange [-50:50]
set grid
plot 'stat.txt' using 2:3 title 'segment-length vs delta x' with points ps 0.1
set yrange [-50:50]
set grid
plot 'stat.txt' using 2:4 title 'segment-length vs delta y' with points ps 0.1
set yrange [-50:50]
set grid
plot 'stat.txt' using 2:5 title 'segment-length vs delta radius' with points ps 0.1