Linear regression: least squares, Tellson estimation, RANSAC

Least square method （LSM）

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

Matrix method formula ：

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); //

Tyson estimates (Theil–Sen estimator)

Calculate the slope in pairs , And then sort it from small to large , Take the median

// Press x Sort from large to small , perhaps [x(t),y(t)], Press t Sort 
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];
		
		// Separate type 
		x_delta_histogram[x1 - x2]++; //  Add histogram , Easy to sort 
		y_delta_histogram[y1 - y2]++; // 
		// Combined 
		//deltaxy_histogram[(y1 - y2)/(x1 - x2)]++;
	}
}

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

RANSAC Algorithm

RANSAC The input of the algorithm is a set of observation data （ It often contains large noise or invalid points ）, It is a resampling technique （resampling technique）, By estimating the minimum number of sample points required for model parameters , To get a set of alternative models , And then continue to expand the set , The algorithm steps are ：

1. Randomly select the minimum sample points required to estimate the model parameters ( A straight line is two points , A circle is three points ).
2. Estimate the parameters of the model （ Calculate a line or circle equation ）.
3. Find out if it is within the error , How many points fit the current model , And mark these points as interior points of the model
4. If the ratio of the number of interior points to the total sample points reaches the preset threshold , Then re estimate the parameters of the model based on these interior points ( least square ), And take this as the final model , To terminate the program .
5. Otherwise repeat 1 To 4 Step .

RANSAC The algorithm learns the model from a random subset of the interior points of the input sample set .

RANSAC The algorithm is a non deterministic algorithm （non-deterministic algorithm）, This algorithm can only get a good result with a certain probability , When the basic model has been determined , The result depends on the maximum number of iterations .

Omit re estimating the model version ：
（1） Randomly select two points from the observation points , Get a straight line through this point ;
（2） use （1） To test other observation points , The distance from the point to the line determines whether the observation point is an internal point or an external point ;
（3） If there are enough local points , And the local points are more than the original “ The best ” The local interior point of a line , Then set the iteration line as “ The best ” A straight line ;
（4） repeat （1）~（3） Step until the best line is found .

Line estimation ：

maxnum=0;
bestline;
// Maximum iteration 50 Time 
for(int c=0;c< 50;c++)
{
    
	// Random index 
	idx1= get_randon(); 
	idx2= get_randon();
	innum=2;// Number of interior points 
	// Calculate line parameters 
	line= calculate(point[idx1] , point[idx2]);
	inpoints=[];
	for(int i=0;i< num;i++)
	{
    
	 	if((i!=idx1)&&(i!=idx2))
	 	{
    
 			// Calculated distance 
 			disdance = pointdis(line,point[i]);
 			if(disdance < threshold)
 			{
    
 				innum++;
 				inpoints=[inpoints,point[i]];// Add the point to the inner point queue 
 			}
		}
	}
	// find n The best line for the second iteration 
	if(innum > maxnum)
	{
    
		bestline= line;
	}
	if(innum > numthreshold)
	{
    
	 	// Recalculate line parameters 
		bestline= calculate_line(inpoints[]);   
	    break;
	}
}

return bestline;

Two of them calculate the linear equation as follows ：

//（ If the line is y=ax+b form , be C by Nan; If the line is x=c form , be A and B by 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;
     }
}

Calculate the distance from the point to the line ：

double pointdis(LineF line,PointF p)
{
    
    double d = 0d;
    if (line.C = NaN)
    {
    
         //y=ax+b amount to 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;
}

Circle estimation ：

maxnum=0;
bestcircle;
// Maximum iteration 50 Time 
for(int c=0;c< 50;c++)
{
    
	// Random index 
	idx1= get_randon(); 
	idx2= get_randon();
	idx3= get_randon();
	innum=2;// Number of interior points 
	// Calculate circle parameters 
	circle= calculate(point[idx1] , point[idx2], point[idx2]);
	inpoints=[];
	for(int i=0;i< num;i++)
	{
    
	 	if((i!=idx1)&&(i!=idx2)&&(i!=idx3))
	 	{
    
 			// Calculated distance 
 			disdance = pointdis(circle,point[i]);
 			if(disdance < threshold)
 			{
    
 				innum++;
 				inpoints=[inpoints,point[i]];// Add the point to the inner point queue 
 			}
		}
	}
	// find n The best circle of iterations 
	if(innum > maxnum)
	{
    
		bestcircle= circle;
	}
	if(innum > numthreshold)
	{
    
	 	// Recalculate circle parameters 
		bestcircle= calculate_line(inpoints[]);   
	    break;
	}
}

return bestline;

The distance from the point to the circle is ：

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

The equation for constructing a circle at three points ：

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(" Parameter error , The three points provided cannot form a circle .");
    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;

}