线性回归：最小二乘、泰尔森估计、RANSAC

最小二乘法（LSM）

x=[ 0,1,2,3…]
y=[ 2,4,6,8…]

矩阵法公式：

for (int i = 0; i < num; i++)
{
    
	sum_x += x[i], 
	sum_xx += x[i] * x[i], 
	sum_y += y[i], 
	sum_xy += x[i] * y[i];
}
	a0 = (sum_x * sum_xx - sum_xy * sum_x) / (sum_xx * area  - sum_x  * sum_x);
	a1 = (sum_y * sum_x  - sum_xy * num) / (sum_x  * sum_x - sum_xx * num); //

泰尔森估算(Theil–Sen estimator)

两两计算斜率，然后从小到大排序，取中值

//按x从大到小排序,或者[x(t),y(t)],按t排序
sort( [x,y] );

for (int i = 0; i < num; i++) {
    
	x1 =x[i];y1 =y[i];
	for (int j = i + 1; j < num; j++) {
    
		x2 = x[j];y2 =y[j];
		
		//分离式
		x_delta_histogram[x1 - x2]++; // 加入直方图，方便排序
		y_delta_histogram[y1 - y2]++; // 
		//合并式
		//deltaxy_histogram[(y1 - y2)/(x1 - x2)]++;
	}
}

mediank= getmid(x_delta_histogram) /getmid(x_delta_histogram)
//mediank= getmid(deltaxy_histogram)

RANSAC算法

RANSAC算法的输入是一组观测数据（往往含有较大的噪声或无效点），它是一种重采样技术（resampling technique），通过估计模型参数所需的最小的样本点数，来得到备选模型集合，然后在不断的对集合进行扩充，其算法步骤为：

1. 随机的选择估计模型参数所需的最少的样本点 ( 直线的话就是两个点，圆是三个点)。
2. 估计出模型的参数（计算直线或圆方程）。
3. 找出在误差内，有多少点适合当前这个模型，并将这些点标记为模型内点
4. 如果内点的数目占总样本点的比例达到了事先设定的阈值 ，那么基于这些内点重新估计模型的参数(最小二乘)，并以此为最终模型， 终止程序。
5. 否则重复执行1到4步。

RANSAC算法是从输入样本集合的内点的随机子集中学习模型。

RANSAC算法是一个非确定性算法（non-deterministic algorithm），这个算法只能得以一定的概率得到一个还不错的结果，在基本模型已定的情况下，结果的好坏程度主要取决于算法最大的迭代次数。

省略重新估计模型版本：
（1）随机从观测点中选择两个点，得到通过该点的直线；
（2）用（1）中的直线去测试其他观测点，由点到直线的距离确定观测点是否为局内点或者局外点；
（3）如果局内点足够多，并且局内点多于原有“最佳”直线的局内点，那么将这次迭代的直线设为“最佳”直线；
（4）重复（1）~（3）步直到找到最佳直线。

直线估计：

maxnum=0;
bestline;
//最大迭代50次
for(int c=0;c< 50;c++)
{
    
	//随机索引
	idx1= get_randon(); 
	idx2= get_randon();
	innum=2;//内点数量
	//计算直线参数
	line= calculate(point[idx1] , point[idx2]);
	inpoints=[];
	for(int i=0;i< num;i++)
	{
    
	 	if((i!=idx1)&&(i!=idx2))
	 	{
    
 			//计算距离
 			disdance = pointdis(line,point[i]);
 			if(disdance < threshold)
 			{
    
 				innum++；
 				inpoints=[inpoints,point[i]];//将该点加入内点队列
 			}
		}
	}
	//找到n次迭代的最佳直线
	if(innum > maxnum)
	{
    
		bestline= line;
	}
	if(innum > numthreshold)
	{
    
	 	//重新计算直线参数
		bestline= calculate_line(inpoints[]);   
	    break;
	}
}

return bestline;

其中两点计算直线方程如下：

//（如果直线为y=ax+b形式，则C为Nan；如果直线为x=c形式，则A和B为Nan）

calculate(PointF p1, PointF p2)
{
    
     if (p1.X == p2.X)
     {
    
        A = NaN;
        B = NaN;
        C = p1.X;
     }
     else
     {
    
        A = 1d * (p1.Y - p2.Y) / (p1.X - p2.X);
        B = p1.Y - A * p1.X;
        C = NaN;
     }
}

计算点到直线的距离：

double pointdis(LineF line,PointF p)
{
    
    double d = 0d;
    if (line.C = NaN)
    {
    
         //y=ax+b相当于ax-y+b=0
         d = abs(1d * (line.A * p.X - p.Y + line.B) / Sqrt(line.A * line.A + 1*1));
     }
     else
     {
    
         d = abs(C - p.X);
     }
     return d;
}

圆估计：

maxnum=0;
bestcircle;
//最大迭代50次
for(int c=0;c< 50;c++)
{
    
	//随机索引
	idx1= get_randon(); 
	idx2= get_randon();
	idx3= get_randon();
	innum=2;//内点数量
	//计算圆参数
	circle= calculate(point[idx1] , point[idx2], point[idx2]);
	inpoints=[];
	for(int i=0;i< num;i++)
	{
    
	 	if((i!=idx1)&&(i!=idx2)&&(i!=idx3))
	 	{
    
 			//计算距离
 			disdance = pointdis(circle,point[i]);
 			if(disdance < threshold)
 			{
    
 				innum++；
 				inpoints=[inpoints,point[i]];//将该点加入内点队列
 			}
		}
	}
	//找到n次迭代的最佳圆
	if(innum > maxnum)
	{
    
		bestcircle= circle;
	}
	if(innum > numthreshold)
	{
    
	 	//重新计算圆参数
		bestcircle= calculate_line(inpoints[]);   
	    break;
	}
}

return bestline;

其中点到圆的距离为：

double pointdis(PointF p)
{
    
    return abs(Cr - Sqrt( Pow(p.X - Cx, 2) + Pow(p.Y - Cy, 2)));
}

三点构造圆的方程：

void calculate(PointF p1,PointF p2,PointF p3)
{
    
    float xMove = p1.X;
    float yMove = p1.Y;
    
    p1.X = 0;
    p1.Y = 0;
    
    p2.X = p2.X - xMove;
    p2.Y = p2.Y - yMove;
    
    p3.X = p3.X - xMove;
    p3.Y = p3.Y - yMove;
    
    float x1 = p2.X, y1 = p2.Y, x2 = p3.X, y2 = p3.Y;
    float m = 2.0 * (x1 * y2 - y1 * x2);
    if (m == 0)
       throw new ArgumentException("参数错误，提供的三点不能构成圆。");
    float  x0 = (x1 * x1 * y2 - x2 * x2 * y1 + y1 * y2 * (y1 - y2)) / m;
    float y0 = (x1 * x2 * (x2 - x1) - y1 * y1 * x2 + x1 * y2 * y2) / m;
    
    Cr = Sqrt(x0 * x0 + y0 * y0);
    Cx = x0 + xMove;
    Cy = y0 + yMove;

}