圆拟合Taubin fit 方法

参考网站：

https://people.cas.uab.edu/~mosya/cl/

C++版本源码：

#pragma once
//头文件引用声明
#include <stdio.h>
#include <tchar.h>
#include <iostream>
#include <cmath>
#include <limits>
#include <iomanip>
#include <cstdlib>using namespace std;typedef double reals;typedef long long integers;//常用数学常量声明
const reals One = 1.0, Two = 2.0, Three = 3.0, Four = 4.0, Five = 5.0, Six = 6.0, Ten = 10.0;
const reals Pi = 3.141592653589793238462643383L;template<typename T>
inline T SQR(T t) { return t*t; };//data数据类型定义
class Data
{
public:int n;reals *X;      //space is allocated in the constructorsreals *Y;       //space is allocated in the constructorsreals meanX, meanY;// constructorsData();Data(int N);Data(int N, reals X[], reals Y[]);// routinesvoid means(void);void center(void);void scale(void);void print(void);// destructors~Data();
};/************************************************************************
BODY OF THE MEMBER ROUTINES
************************************************************************/
// Default constructor
Data::Data()
{n = 0;X = new reals[n];Y = new reals[n];for (int i = 0; i<n; i++){X[i] = 0.;Y[i] = 0.;}
}// Constructor with assignment of the field N
Data::Data(int N)
{n = N;X = new reals[n];Y = new reals[n];for (int i = 0; i<n; i++){X[i] = 0.;Y[i] = 0.;}
}// Constructor with assignment of each field
Data::Data(int N, reals dataX[], reals dataY[])
{n = N;X = new reals[n];Y = new reals[n];for (int i = 0; i<n; i++){X[i] = dataX[i];Y[i] = dataY[i];}
}// Routine that computes the x- and y- sample means (the coordinates of the centeroid)void Data::means(void)
{meanX = 0.; meanY = 0.;for (int i = 0; i<n; i++){meanX += X[i];meanY += Y[i];}meanX /= n;meanY /= n;
}// Routine that centers the data set (shifts the coordinates to the centeroid)void Data::center(void)
{reals sX = 0., sY = 0.;int i;for (i = 0; i<n; i++){sX += X[i];sY += Y[i];}sX /= n;sY /= n;for (i = 0; i<n; i++){X[i] -= sX;Y[i] -= sY;}meanX = 0.;meanY = 0.;
}// Routine that scales the coordinates (makes them of order one)void Data::scale(void)
{reals sXX = 0., sYY = 0., scaling;int i;for (i = 0; i<n; i++){sXX += X[i] * X[i];sYY += Y[i] * Y[i];}scaling = sqrt((sXX + sYY) / n / Two);for (i = 0; i<n; i++){X[i] /= scaling;Y[i] /= scaling;}
}// Printing routinevoid Data::print(void)
{cout << endl << "The data set has " << n << " points with coordinates :" << endl;for (int i = 0; i<n - 1; i++) cout << setprecision(7) << "(" << X[i] << "," << Y[i] << "), ";cout << "(" << X[n - 1] << "," << Y[n - 1] << ")\n";
}// Destructor
Data::~Data()
{delete[] X;delete[] Y;
}//Circle.h
class Circle
{
public:// The fields of a Circlereals a, b, r, s, g, Gx, Gy;int i, j;// constructorsCircle();Circle(reals aa, reals bb, reals rr);// routinesvoid print(void);// no destructor we didn't allocate memory by hand.
};/************************************************************************
BODY OF THE MEMBER ROUTINES
************************************************************************/
// Default constructorCircle::Circle()
{a = 0.; b = 0.; r = 1.; s = 0.; i = 0; j = 0;
}// Constructor with assignment of the circle parameters onlyCircle::Circle(reals aa, reals bb, reals rr)
{a = aa; b = bb; r = rr;
}// Printing routinevoid Circle::print(void)
{cout << endl;cout << setprecision(10) << "center (" << a << "," << b << ")  radius "<< r << "  sigma " << s << "  gradient " << g << "  iter " << i << "  inner " << j << endl;
}//函数声明
//****************** Sigma ************************************
//
//   estimate of Sigma = square root of RSS divided by N
//   gives the root-mean-square error of the geometric circle fitreals Sigma(Data& data, Circle& circle)
{reals sum = 0., dx, dy;for (int i = 0; i<data.n; i++){dx = data.X[i] - circle.a;dy = data.Y[i] - circle.b;sum += SQR(sqrt(dx*dx + dy*dy) - circle.r);}return sqrt(sum / data.n);
}Circle CircleFitByTaubin(Data& data)
/*
参数声明：
Input:  data     - the class of data (contains the given points):data.n   - the number of data points
data.X[] - the array of X-coordinates
data.Y[] - the array of Y-coordinatesOutput:
circle - parameters of the fitting circle:circle.a - the X-coordinate of the center of the fitting circle
circle.b - the Y-coordinate of the center of the fitting circle
circle.r - the radius of the fitting circle
circle.s - the root mean square error (the estimate of sigma)
circle.j - the total number of iterations*/
{int i, iter, IterMAX = 99;reals Xi, Yi, Zi;reals Mz, Mxy, Mxx, Myy, Mxz, Myz, Mzz, Cov_xy, Var_z;reals A0, A1, A2, A22, A3, A33;reals Dy, xnew, x, ynew, y;reals DET, Xcenter, Ycenter;Circle circle;data.means();   // Compute x- and y- sample means (via a function in the class "data")//     computing momentsMxx = Myy = Mxy = Mxz = Myz = Mzz = 0.;for (i = 0; i<data.n; i++){Xi = data.X[i] - data.meanX;   //  centered x-coordinatesYi = data.Y[i] - data.meanY;   //  centered y-coordinatesZi = Xi*Xi + Yi*Yi;Mxy += Xi*Yi;Mxx += Xi*Xi;Myy += Yi*Yi;Mxz += Xi*Zi;Myz += Yi*Zi;Mzz += Zi*Zi;}Mxx /= data.n;Myy /= data.n;Mxy /= data.n;Mxz /= data.n;Myz /= data.n;Mzz /= data.n;//      computing coefficients of the characteristic polynomialMz = Mxx + Myy;Cov_xy = Mxx*Myy - Mxy*Mxy;Var_z = Mzz - Mz*Mz;A3 = Four*Mz;A2 = -Three*Mz*Mz - Mzz;A1 = Var_z*Mz + Four*Cov_xy*Mz - Mxz*Mxz - Myz*Myz;A0 = Mxz*(Mxz*Myy - Myz*Mxy) + Myz*(Myz*Mxx - Mxz*Mxy) - Var_z*Cov_xy;A22 = A2 + A2;A33 = A3 + A3 + A3;for (x = 0., y = A0, iter = 0; iter<IterMAX; iter++)  // usually, 4-6 iterations are enough{Dy = A1 + x*(A22 + A33*x);xnew = x - y / Dy;if ((xnew == x) || (!_finite(xnew)/*!isinf(xnew)*/)) break;ynew = A0 + xnew*(A1 + xnew*(A2 + xnew*A3));if (abs(ynew) >= abs(y))  break;x = xnew;  y = ynew;}//       computing paramters of the fitting circleDET = x*x - x*Mz + Cov_xy;Xcenter = (Mxz*(Myy - x) - Myz*Mxy) / DET / Two;Ycenter = (Myz*(Mxx - x) - Mxz*Mxy) / DET / Two;//       assembling the outputcircle.a = Xcenter + data.meanX;circle.b = Ycenter + data.meanY;circle.r = sqrt(Xcenter*Xcenter + Ycenter*Ycenter + Mz);circle.s = Sigma(data, circle);circle.i = 0;circle.j = iter;  //  return the number of iterations, tooreturn circle;}

Matlab：

function Par = TaubinNTN(XY)%--------------------------------------------------------------------------
%
%     Algebraic circle fit by Taubin, based on Newton's method
%      G. Taubin, "Estimation Of Planar Curves, Surfaces And Nonplanar
%                  Space Curves Defined By Implicit Equations, With
%                  Applications To Edge And Range Image Segmentation",
%      IEEE Trans. PAMI, Vol. 13, pages 1115-1138, (1991)
%
%     Input:  XY(n,2) is the array of coordinates of n points x(i)=XY(i,1), y(i)=XY(i,2)
%
%     Output: Par = [a b R] is the fitting circle:
%                           center (a,b) and radius R
%
%     Note: this is a version optimized for speed, not for stability
%           it does not use built-in matrix functions,
%           so it can be easily programmed in any language
%
%--------------------------------------------------------------------------n = size(XY,1);      % number of data pointscentroid = mean(XY);   % the centroid of the data set%     computing moments (note: all moments will be normed, i.e. divided by n)Mxx = 0; Myy = 0; Mxy = 0; Mxz = 0; Myz = 0; Mzz = 0;for i=1:nXi = XY(i,1) - centroid(1);  %  centering dataYi = XY(i,2) - centroid(2);  %  centering dataZi = Xi*Xi + Yi*Yi;Mxy = Mxy + Xi*Yi;Mxx = Mxx + Xi*Xi;Myy = Myy + Yi*Yi;Mxz = Mxz + Xi*Zi;Myz = Myz + Yi*Zi;Mzz = Mzz + Zi*Zi;
endMxx = Mxx/n;
Myy = Myy/n;
Mxy = Mxy/n;
Mxz = Mxz/n;
Myz = Myz/n;
Mzz = Mzz/n;%    computing the coefficients of the characteristic polynomialMz = Mxx + Myy;
Cov_xy = Mxx*Myy - Mxy*Mxy;
A3 = 4*Mz;
A2 = -3*Mz*Mz - Mzz;
A1 = Mzz*Mz + 4*Cov_xy*Mz - Mxz*Mxz - Myz*Myz - Mz*Mz*Mz;
A0 = Mxz*Mxz*Myy + Myz*Myz*Mxx - Mzz*Cov_xy - 2*Mxz*Myz*Mxy + Mz*Mz*Cov_xy;
A22 = A2 + A2;
A33 = A3 + A3 + A3;xnew = 0;
ynew = 1e+20;
epsilon = 1e-12;
IterMax = 20;% Newton's method starting at x=0for iter=1:IterMaxyold = ynew;ynew = A0 + xnew*(A1 + xnew*(A2 + xnew*A3));if abs(ynew) > abs(yold)disp('Newton-Taubin goes wrong direction: |ynew| > |yold|');xnew = 0;break;endDy = A1 + xnew*(A22 + xnew*A33);xold = xnew;xnew = xold - ynew/Dy;if (abs((xnew-xold)/xnew) < epsilon), break, endif (iter >= IterMax)disp('Newton-Taubin will not converge');xnew = 0;endif (xnew<0.)fprintf(1,'Newton-Taubin negative root:  x=%f\n',xnew);xnew = 0;end
end%  computing the circle parametersDET = xnew*xnew - xnew*Mz + Cov_xy;
Center = [Mxz*(Myy-xnew)-Myz*Mxy , Myz*(Mxx-xnew)-Mxz*Mxy]/DET/2;Par = [Center+centroid , sqrt(Center*Center'+Mz)];end    %    TaubinNTN

圆拟合Taubin fit 方法相关推荐

圆拟合算法（距离之和最小）
上一篇博客介绍了最小二乘法拟合圆的方法.这种方法对误差符合正态分布的数据点很有效.但是在机器视觉应用中经常会碰到一些干扰点.这些干扰点多数时候是偏向某一个方向的.这时要是用最小二乘法拟合,拟合出的圆会 ...
最小二乘法-圆拟合（不啰嗦）
原理:原理部分网上大部分可以搜得到,以一句很简单的话就是是通过最小化误差的平方和找到一组数据的最佳函数匹配(自行百度). 作用:如果现在有一张图片,需要你拟合图片中的圆. 需要拟合的圆图片: 方法:最 ...
防止模型过拟合的必备方法！
↑↑↑关注后"星标"Datawhale 每日干货 & 每月组队学习,不错过 Datawhale干货作者:Mahitha,来源:机器之心正如巴菲特所言:「近似的正确好过精 ...
python求交点坐标_Python求两个圆的交点坐标或三个圆的交点坐标方法
计算两个圆的交点代码如下: # -*- coding: utf-8 -*- import math import numpy as np def insec(p1,r1,p2,r2): x = p1 ...
直线和圆交点 halcon_初中数学三角形、四边形、圆辅助线的添加方法
今天,小编为大家整理了初中数学三角形.四边形.圆的辅助线添加方法,速来看!! 1三角形中常见辅助线的添加与角平分线有关的 (1)可向两边作垂线: (2)可作平行线,构造等腰三角形: (3)在角的两边 ...
欠拟合、过拟合及其解决方法
欠拟合.过拟合及其解决方法参考文章: (1)欠拟合.过拟合及其解决方法 (2)https://www.cnblogs.com/alan666/p/8311809.html 备忘一下.
机器学习中防止过拟合的处理方法
原文地址:一只鸟的天空,http://blog.csdn.net/heyongluoyao8/article/details/49429629 防止过拟合的处理方法过拟合我们都知道,在进行数据 ...
防止过拟合的处理方法
原文地址:一只鸟的天空,http://blog.csdn.net/heyongluoyao8/article/details/49429629 防止过拟合的处理方法过拟合我们都知道,在进行数据 ...
解决神经网络过拟合问题—Dropout方法、python实现
解决神经网络过拟合问题-Dropout方法一.what is Dropout?如何实现? 二.使用和不使用Dropout的训练结果对比一.what is Dropout?如何实现? 如果网络模型复 ...

圆拟合Taubin fit 方法

参考网站：

C++版本源码：

Matlab：

圆拟合Taubin fit 方法相关推荐

最新文章

热门文章