基于opencv求直线方程

在图像处理中，得到目标边缘的点，希望拟合精确的直线方程，一下代码给出一些二维点，即可输出直线方程，精确到亚像素级：

 

//求直线方程
struct PtOnLine
{
    double d;
    cv::Point2d pt;
};
bool cmpACS(PtOnLine &a, PtOnLine &b)
{
    return a.d > b.d;
}
void RemovePtX(vector<Point> &points, double smallpart, double bigpart)//对X坐标从大到小排列,smallpart为尾部的点集，bigpart为首部的点集
{
    int removeNum1 = bigpart * points.size();
    int removeNum2 = smallpart * points.size();
    vector<PtOnLine> vecline(points.size());
    for (int i = 0; i < points.size(); i++)
    {
        PtOnLine ptOnCircle;
        vecline[i].d = points[i].x;
        vecline[i].pt = points[i];
    }
    sort(vecline.begin(), vecline.end(), cmpACS);
    points.clear();
    for (int i = removeNum1; i < vecline.size() - removeNum2; i++)
    {
        points.push_back(vecline[i].pt);
    }
}
void RemovePtY(vector<Point> &points, double smallpart, double bigpart)//对Y坐标从大到小排列，smallpart为尾部的点集，bigpart为首部的点集
{
    int removeNum1 = bigpart * points.size();
    int removeNum2 = smallpart * points.size();
    vector<PtOnLine> vecline(points.size());
    for (int i = 0; i < points.size(); i++)
    {
        PtOnLine ptOnCircle;
        vecline[i].d = points[i].y;
        vecline[i].pt = points[i];
    }
    sort(vecline.begin(), vecline.end(), cmpACS);
    points.clear();
    for (int i = removeNum1; i < vecline.size() - removeNum2; i++)
    {
        points.push_back(vecline[i].pt);
    }
}
bool LineFit(vector<Point>rec, double &a1, double &b1, double &c1)
{
    int size = rec.size();
    double Xall = 0;
    double count = 0;
    for (int i = 0; i < size; i++)
    {
        Xall = Xall + rec[i].x;
    }
    for (int i = 0; i < size; i++)
    {
        if (rec[i].x == Xall / (double)size)
        {
            count++;
        }
    }
    if (count == size)
    {
        a1 = 1;
        b1 = 0;
        c1 = -Xall / size;
        return true;
    }
    double a = 0, b = 0, c = 0;
    if (size < 2)
    {
        a1 = 0;
        b1 = 0;
        c1 = 0;
        return false;
    }

    double x_mean = 0;
    double y_mean = 0;
    for (int i = 0; i < size; i++)
    {
        x_mean += rec[i].x;
        y_mean += rec[i].y;
    }
    x_mean /= size;
    y_mean /= size; //至此，计算出了 x y 的均值
    double Dxx = 0, Dxy = 0, Dyy = 0;
    for (int i = 0; i < size; i++)
    {
        Dxx += (rec[i].x - x_mean) * (rec[i].x - x_mean);
        Dxy += (rec[i].x - x_mean) * (rec[i].y - y_mean);
        Dyy += (rec[i].y - y_mean) * (rec[i].y - y_mean);
    }
    double lambda = ((Dxx + Dyy) - sqrt((Dxx - Dyy) * (Dxx - Dyy) + 4 * Dxy * Dxy)) / 2.0;
    double den = sqrt(Dxy * Dxy + (lambda - Dxx) * (lambda - Dxx));
    a = Dxy / den;
    b = (lambda - Dxx) / den;
    c = -a * x_mean - b * y_mean;
    for (int i = 0; i<size-1; i++)
    {
        if (fabs(a*rec[i].x + b*rec[i].y + c) / sqrt(a*a + b*b)>1.5)
        {
            rec[i].x = rec[i + 1].x;
            rec[i].y = rec[i + 1].y;
            size--;
        }
    }
    vector<Point>recNew(size);
    for (int i = 0; i < size; i++)
    {
        recNew[i].x = rec[i].x;
        recNew[i].y = rec[i].y;
    }
    if (size < 2)
    {
        a1 = 0;
        b1 = 0;
        c1 = 0;
        return false;
    }
    double x_mean1 = 0;
    double y_mean1 = 0;
    for (int i = 0; i < size; i++)
    {
        x_mean1 += recNew[i].x;
        y_mean1 += recNew[i].y;
    }
    x_mean1 /= size;
    y_mean1 /= size; //至此，计算出了 x y 的均值
    double Dxx1 = 0, Dxy1 = 0, Dyy1 = 0;
    for (int i = 0; i < size; i++)
    {
        Dxx1 += (recNew[i].x - x_mean1) * (recNew[i].x - x_mean1);
        Dxy1 += (recNew[i].x - x_mean1) * (recNew[i].y - y_mean1);
        Dyy1 += (recNew[i].y - y_mean1) * (recNew[i].y - y_mean1);
    }
    double lambda1 = ((Dxx1 + Dyy1) - sqrt((Dxx1 - Dyy1) * (Dxx1 - Dyy1) + 4 * Dxy1 * Dxy1)) / 2.0;
    double den1 = sqrt(Dxy * Dxy + (lambda - Dxx) * (lambda - Dxx));
    a1 = Dxy1 / den1;
    b1 = (lambda1 - Dxx1) / den1;
    c1 = -a1 * x_mean1 - b1 * y_mean1;
    return true;
}