圆拟合 --- 点到圆距离绝对值之和最小

2024-01-26 19:12

本文主要是介绍圆拟合 --- 点到圆距离绝对值之和最小,希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!

https://www.researchgate.net/publication/220146330_Robust_Fitting_of_Circle_Arcs 根据这篇文章介绍的算法进行了实现,经过验证结果与论文中结果一致

AGE.h

#pragma once
#include <vector>struct Point2d {double x;double z;
};struct Circle2d
{Point2d center;double radius;
};// https://www.researchgate.net/publication/220146330_Robust_Fitting_of_Circle_Arcs
bool fitCircleAGE(const std::vector<Point2d>& points, const double epsilon, double lambda,int iterations, Circle2d& circle);

AGE.cpp

#include <cmath>
#include <algorithm>
#include <numeric>#include "AGE.h"bool calcCircleCoefficients(const std::vector<Point2d>& points, Circle2d& circle)
{if (points.size() != 3)return false;const Point2d p1 = points[0], p2 = points[1], p3 = points[2];const double a = 2 * (p2.x - p1.x), b = 2 * (p2.z - p1.z);const double c = p2.x * p2.x + p2.z * p2.z - p1.x * p1.x - p1.z * p1.z;const double d = 2 * (p3.x - p2.x), e = 2 * (p3.z - p2.z);const double f = p3.x * p3.x + p3.z * p3.z - p2.x * p2.x - p2.z * p2.z;const double x = (b * f - e * c) / (b * d - e * a);const double z = (d * c - a * f) / (b * d - e * a);const double r = sqrt((x - p1.x) * (x - p1.x) + (z - p1.z) * (z - p1.z));circle = { {x, z}, r };return true;
}double calcDistancePointToPoint(const Point2d& p1, const Point2d& p2)
{return std::sqrt(std::pow(p1.x - p2.x, 2) + std::pow(p1.z - p2.z, 2));
}double calcMedian(const std::vector<double>& values)
{auto vec = values;const int size = vec.size();const int mid = size / 2;if (size & 1) {std::nth_element(vec.begin(), vec.begin() + mid, vec.end());return vec[mid];}std::nth_element(vec.begin(), vec.begin() + mid - 1, vec.end());return (vec[mid - 1] + *std::min_element(vec.begin() + mid, vec.end())) / 2.0;
}// https://www.researchgate.net/publication/220146330_Robust_Fitting_of_Circle_Arcs
bool fitCircleAGE(const std::vector<Point2d>& points, const double epsilon, double lambda,int iterations, Circle2d& circle)
{if (points.size() < 3)return false;// step 1: initail center of circle.std::vector<Point2d> points3 = { points[0], points[1], points[2] };calcCircleCoefficients(points3, circle);lambda = 1;const int pointsNum = points.size();std::vector<double> distances(pointsNum);std::vector<int> I, O, C;I.reserve(pointsNum);O.reserve(pointsNum);C.reserve(pointsNum);std::vector<double> cosTheta(pointsNum), sinTheta(pointsNum);double error = DBL_MAX;while (iterations-- > 0) {// step 2: compute the distances between each point to the circle// center and the median r of these distances.
#pragma omp parallel forfor (int i = 0; i < pointsNum; i++)distances[i] = calcDistancePointToPoint(circle.center, points[i]);circle.radius = calcMedian(distances);// step 3: determine the sets I = {i : distance[i] < r}, O = {i : distance[i] > r}, C = {i :// distance[i] == r}.I.clear();O.clear();C.clear();
#pragma omp parallel{std::vector<int> localI, localO, localC;#pragma omp forfor (int i = 0; i < pointsNum; i++) {if (distances[i] < circle.radius)localI.push_back(i);else if (distances[i] > circle.radius)localO.push_back(i);elselocalC.push_back(i);}
#pragma omp critical{I.insert(I.end(), localI.begin(), localI.end());O.insert(O.end(), localO.begin(), localO.end());C.insert(C.end(), localC.begin(), localC.end());}}// step 4: calculate the values cos(θ_i).
#pragma omp parallel forfor (int i = 0; i < pointsNum; i++) {cosTheta[i] = (points[i].x - circle.center.x) / distances[i];sinTheta[i] = (points[i].z - circle.center.z) / distances[i];}// step 5: compute EaPos, EaNeg, EbPos, EbNeg.double sumCosThetaI = std::accumulate(I.begin(), I.end(), 0.0, [&](double a, int i) { return a + cosTheta[i]; });double sumSinThetaI = std::accumulate(I.begin(), I.end(), 0.0, [&](double a, int i) { return a + sinTheta[i]; });double sumCosThetaO = std::accumulate(O.begin(), O.end(), 0.0, [&](double a, int i) { return a + cosTheta[i]; });double sumSinThetaO = std::accumulate(O.begin(), O.end(), 0.0, [&](double a, int i) { return a + sinTheta[i]; });double sumCosThetaC = std::accumulate(C.begin(), C.end(), 0.0, [&](double a, int i) { return a + abs(cosTheta[i]); });double sumSinThetaC = std::accumulate(C.begin(), C.end(), 0.0, [&](double a, int i) { return a + abs(sinTheta[i]); });double EaPos = sumCosThetaI - sumCosThetaO + sumCosThetaC;double EaNeg = sumCosThetaI - sumCosThetaO - sumCosThetaC;double EbPos = sumSinThetaI - sumSinThetaO + sumSinThetaC;double EbNeg = sumSinThetaI - sumSinThetaO - sumSinThetaC;// step 6: determine the new center.double alpha = EaNeg >= 0 ? 1 : EaPos <= 0 ? 0 : 0.5;double beta = EbNeg >= 0 ? 1 : EbPos <= 0 ? 0 : 0.5;double d1 = -(alpha * EaNeg + (1 - alpha) * EaPos);double d2 = -(beta * EbNeg + (1 - beta) * EbPos);circle.center.x += d1 * lambda;circle.center.z += d2 * lambda;if (std::sqrt(d1 * d1 + d2 * d2) * lambda < epsilon)return true;double newError = std::accumulate(distances.begin(), distances.end(), 0.0,[&](double sum, double b) { return sum + std::abs(b - circle.radius); });if (newError < error) {lambda *= 1.1;error = newError;}else {lambda *= 0.9;}}return true;
}

用了一些加速手段,但好像不明显

这篇关于圆拟合 --- 点到圆距离绝对值之和最小的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



http://www.chinasem.cn/article/647787

相关文章

poj 1258 Agri-Net(最小生成树模板代码)

感觉用这题来当模板更适合。 题意就是给你邻接矩阵求最小生成树啦。~ prim代码:效率很高。172k...0ms。 #include<stdio.h>#include<algorithm>using namespace std;const int MaxN = 101;const int INF = 0x3f3f3f3f;int g[MaxN][MaxN];int n

poj 1287 Networking(prim or kruscal最小生成树)

题意给你点与点间距离,求最小生成树。 注意点是,两点之间可能有不同的路,输入的时候选择最小的,和之前有道最短路WA的题目类似。 prim代码: #include<stdio.h>const int MaxN = 51;const int INF = 0x3f3f3f3f;int g[MaxN][MaxN];int P;int prim(){bool vis[MaxN];

poj 2349 Arctic Network uva 10369(prim or kruscal最小生成树)

题目很麻烦,因为不熟悉最小生成树的算法调试了好久。 感觉网上的题目解释都没说得很清楚,不适合新手。自己写一个。 题意:给你点的坐标,然后两点间可以有两种方式来通信:第一种是卫星通信,第二种是无线电通信。 卫星通信:任何两个有卫星频道的点间都可以直接建立连接,与点间的距离无关; 无线电通信:两个点之间的距离不能超过D,无线电收发器的功率越大,D越大,越昂贵。 计算无线电收发器D

poj 1734 (floyd求最小环并打印路径)

题意: 求图中的一个最小环,并打印路径。 解析: ans 保存最小环长度。 一直wa,最后终于找到原因,inf开太大爆掉了。。。 虽然0x3f3f3f3f用memset好用,但是还是有局限性。 代码: #include <iostream>#include <cstdio>#include <cstdlib>#include <algorithm>#incl

hdu 1102 uva 10397(最小生成树prim)

hdu 1102: 题意: 给一个邻接矩阵,给一些村庄间已经修的路,问最小生成树。 解析: 把已经修的路的权值改为0,套个prim()。 注意prim 最外层循坏为n-1。 代码: #include <iostream>#include <cstdio>#include <cstdlib>#include <algorithm>#include <cstri

poj 2175 最小费用最大流TLE

题意: 一条街上有n个大楼,坐标为xi,yi,bi个人在里面工作。 然后防空洞的坐标为pj,qj,可以容纳cj个人。 从大楼i中的人到防空洞j去避难所需的时间为 abs(xi - pi) + (yi - qi) + 1。 现在设计了一个避难计划,指定从大楼i到防空洞j避难的人数 eij。 判断如果按照原计划进行,所有人避难所用的时间总和是不是最小的。 若是,输出“OPETIMAL",若

poj 2135 有流量限制的最小费用最大流

题意: 农场里有n块地,其中约翰的家在1号地,二n号地有个很大的仓库。 农场有M条道路(双向),道路i连接着ai号地和bi号地,长度为ci。 约翰希望按照从家里出发,经过若干块地后到达仓库,然后再返回家中的顺序带朋友参观。 如果要求往返不能经过同一条路两次,求参观路线总长度的最小值。 解析: 如果只考虑去或者回的情况,问题只不过是无向图中两点之间的最短路问题。 但是现在要去要回

poj 3422 有流量限制的最小费用流 反用求最大 + 拆点

题意: 给一个n*n(50 * 50) 的数字迷宫,从左上点开始走,走到右下点。 每次只能往右移一格,或者往下移一格。 每个格子,第一次到达时可以获得格子对应的数字作为奖励,再次到达则没有奖励。 问走k次这个迷宫,最大能获得多少奖励。 解析: 拆点,拿样例来说明: 3 2 1 2 3 0 2 1 1 4 2 3*3的数字迷宫,走两次最大能获得多少奖励。 将每个点拆成两个

poj 2195 bfs+有流量限制的最小费用流

题意: 给一张n * m(100 * 100)的图,图中” . " 代表空地, “ M ” 代表人, “ H ” 代表家。 现在,要你安排每个人从他所在的地方移动到家里,每移动一格的消耗是1,求最小的消耗。 人可以移动到家的那一格但是不进去。 解析: 先用bfs搞出每个M与每个H的距离。 然后就是网络流的建图过程了,先抽象出源点s和汇点t。 令源点与每个人相连,容量为1,费用为

poj 3068 有流量限制的最小费用网络流

题意: m条有向边连接了n个仓库,每条边都有一定费用。 将两种危险品从0运到n-1,除了起点和终点外,危险品不能放在一起,也不能走相同的路径。 求最小的费用是多少。 解析: 抽象出一个源点s一个汇点t,源点与0相连,费用为0,容量为2。 汇点与n - 1相连,费用为0,容量为2。 每条边之间也相连,费用为每条边的费用,容量为1。 建图完毕之后,求一条流量为2的最小费用流就行了