2018-12-21 15:44:23 yql_617540298 阅读数 682

一、最小二乘法拟合曲线

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

#自定义函数 e指数形式
def func(x, a, b,c):
    return a*np.sqrt(x)*(b*np.square(x)+c)

#定义x、y散点坐标
x = [10,20,30,40,50,60,70,80]
x = np.array(x)
y = [158,455,265,152,263,813,562,126]
y = np.array(y)

#非线性最小二乘法拟合
popt, pcov = curve_fit(func, x, y)
#获取popt里面是拟合系数
print(popt)
a = popt[0]
b = popt[1]
c = popt[2]
yvals = func(x,a,b,c) #拟合y值
print('popt:', popt)
print('系数a:', a)
print('系数b:', b)
print('系数c:', c)
print('系数pcov:', pcov)
print('系数yvals:', yvals)
#绘图
plot1 = plt.plot(x, y, 's',label='original values')
plot2 = plt.plot(x, yvals, 'r',label='polyfit values')
plt.xlabel('x')
plt.ylabel('y')
plt.legend(loc=4) #指定legend的位置右下角
plt.title('curve_fit')
plt.show()

二、高斯分布拟合曲线

import numpy as np
import math
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit


#自定义函数 e指数形式
def func(x, a,u, sig):
    return  a*(np.exp(-(x - u) ** 2 /(2* sig **2))/(math.sqrt(2*math.pi)*sig))*(431+(4750/x))


#定义x、y散点坐标
x = [10,20,30,40,50,60,70,80]
x=np.array(x)
# x = np.array(range(20))
print('x is :\n',x)
y = [158,455,265,152,263,813,562,126]
y = np.array(y)
print('y is :\n',y)

popt, pcov = curve_fit(func, x, y,p0=[3.1,4.2,3.3])
#获取popt里面是拟合系数
a = popt[0]
u = popt[1]
sig = popt[2]


yvals = func(x,a,u,sig) #拟合y值
print(u'系数a:', a)
print(u'系数u:', u)
print(u'系数sig:', sig)

#绘图
plot1 = plt.plot(x, y, 's',label='original values')
plot2 = plt.plot(x, yvals, 'r',label='polyfit values')
plt.xlabel('x')
plt.ylabel('y')
plt.legend(loc=4) #指定legend的位置右下角
plt.title('curve_fit')
plt.show()

三、多项式拟合曲线

import numpy as np
import matplotlib.pyplot as plt

#定义x、y散点坐标
#x = [10,20,30,40,50,60,70,80]
#y = [158,455,265,152,263,813,562,126]
#x = [16,32,48,64,80,96,112,128,144,160,176,192,208,224,240,256,272,288,304]
x = [0,16,32,48,64,80,96,112,128,144,160,176,192,208,224,240,256,272,288,304,320,336,352,368,384,400,416,432,448,464,480,496,]
x = np.array(x)
#y = [506,506,506,506,506,506,506,506,506,505,505,505,505,505,505,505,506,505,504]
y = [340,338,345,348,348,349,50,350,350,350,350,350,350,350,350,350,350,349,348,348,348,347,347,348,348,347,347,347,347,348,347,346]
y = np.array(y)
#用3次多项式拟合
f1 = np.polyfit(x, y, 3)
p1 = np.poly1d(f1)
yvals = p1(x)#拟合y值

#x1 = [32,48,64,80,96,112,128,144,160,176,192,208,224,240,256,272,288,304,320]
x1 = [16,32,48,64,80,96,112,128,144,160,176,192,208,224,240,256,272,288,304,320,336,352,368,384,400,416,432,448,464,480,496,512]
x1 = np.array(x1)
#y1 = [529,528,528,529,529,529,529,529,529,529,528,528,528,528,528,528,528,528,529]
y1 = [370,367,376,378,378,379,379,380,380,380,379,379,379,379,378,378,378,378,377,376,376,375,375,372,372,372,372,375,375,372,372,373]
y1 = np.array(y1)
f2 = np.polyfit(x1,y1,3)
p2 = np.poly1d(f2)
yvals_2 = p2(x1)
#也可使用yvals=np.polyval(f1, x)
#绘图
plot1 = plt.plot(x, y, 's',label='original values')
plot2 = plt.plot(x, yvals, 'r',label='polyfit values')
plt.xlabel('x')
plt.ylabel('y')

plot3 = plt.plot(x1, y1, 's', label = 'original values')
plot4 = plt.plot(x1, yvals_2, 'r',label='',color='green')

#plt.axis('off')
#plt.legend(loc=4) #指定legend的位置右下角
#plt.title('polyfitting')
#plt.show()
plt.savefig("F:/a.jpg")

2016-07-16 12:43:41 zhuason 阅读数 3141

2016/7/16


 

在一次提取发光管的中心线程序中,由于我们只拍到了断续而弯曲的发光管,所以无法使用光带中心线提取的方法进行提取。

在此背景下,我想到了拟合。之前有学过直线拟合的方法,名为最小二乘法。其基本步骤如下:

(1)    设需要拟合的直线为y=a*x+b。

(2)    首先选取进行拟合的点集,选取方法可以为阈值分割,模板匹配等,设最后选出的点集为。

(3)    求该点集到直线的距离平方和,

(4)    对Sum分别求关于x,y的偏导函数。

(5)    根据偏导求出该距离平方和最小时的a,b值即为拟合的曲线的参数。

为了与之后的二次曲线拟合做对比,我写了下直线拟合的函数,先从源图像中根据阈值分割选取拟合点集(这里我认为灰度值超过45即被选取),再根据上述步骤计算a,b的值。源代码如下:

/*输入为三通道图像*/

/*对图像中的亮点进行直线拟合*/

void cvLineFit1D(IplImage* src_getin)

{

    IplImage*src =cvCloneImage(src_getin);

    IplImage*image_threshold = cvCreateImage(cvGetSize(src),8,1);

    cvCvtColor(src,image_threshold,CV_BGR2GRAY);

    cvThreshold(image_threshold,image_threshold,45,255,CV_THRESH_BINARY);

    cvShowImage("cvLineFit1D[Threshold]",image_threshold);

 

    //设拟合的二次曲线方程为y=ax+b;

    //先求出各点到拟合直线上的距离的平方和;

    //求出使得该平方和最小的a,b的值 ;

    long long int k1=0;

    long long int k2=0;

    long long int k3=0;

    long long int k4=0;

    long long int k5=0;

    long long int k6=0;

    for (inti=0;i<image_threshold->height;i++)

    {

        for (intj=0;j<image_threshold->width;j++)

        {

            if(cvGetReal2D(image_threshold,i,j)==255)

            {

                k1+=2*i*i;

                k2+=2*i;

                k3+=2*i*j;

                k4+=2*i;

                k5+=2;

                k6+=2*j;

            }

        }

    }

 

 

    double a,b;

    a = (double)(k3*k5-k6*k2)/(k1*k5-k4*k2);

    b = (double)(k3-k1*a)/k2;

   

    int bottom = image_threshold->height;

    int top = 0;

    for(int i =0;i<image_threshold->height;i++)

        for(int j =0;j<image_threshold->width;j++)

        {

            if(cvGetReal2D(image_threshold,i,j) == 255)

            {

                if(i>top)

                    top= i;

                if(i<bottom)

                    bottom= i;

            }

        }

 

    for (inti=bottom;i<top;i++)

    {

        int temp_y = a*i+b;

        if(temp_y<src->width&&temp_y>0)

        {

            cvCircle(src,cvPoint(temp_y,i),1,cvScalar(255,255,255));

            //cvSet2D(src,i,temp_y,cvScalar(255,255,255));

        }

        cvShowImage("直线拟合",src);

    }

    cvReleaseImage(&src);

}

拟合效果如下:


 

而曲线拟合,则是在直线拟合的基础上,进行相似的偏导计算,只是设所拟合的曲线方程为y=a*x*x+b*x+c

 

代码如下:

 

/*对图像中的亮点进行二次曲线拟合;*/

void cvCurveFit2D(IplImage* src)

{

 

    IplImage*image_threshold = cvCreateImage(cvGetSize(src),8,1);

    cvCvtColor(src,image_threshold,CV_BGR2GRAY);

 

    cvThreshold(image_threshold,image_threshold,50,255,CV_THRESH_BINARY);

    //cvShowImage("cvCurveFit2D[Threshold]",image_threshold);

      

    //设拟合的二次曲线方程为y=ax2+bx+c;

    //先求出各点到拟合曲线上的距离的平方和;

    //求出使得该平方和最小的a,b,c的值 ;

    double a,b,c;

    long long int k1=0;

    long long int k2=0;

    long long int k3=0;

    long long int k4=0;

    long long int k5=0;

    long long int k6=0;

    long long int k7=0;

    long long int k8=0;

    long long int k9=0;

    long long int k10=0;

    long long int k11=0;

    long long int k12=0;

    for (inti=0;i<image_threshold->height;i+=5)

    {

        for (intj=0;j<image_threshold->width;j+=5)

        {

            if(cvGetReal2D(image_threshold,i,j)==255)

            {

                //cvCircle(src,cvPoint(j,i),1,cvScalar(0,0,0));

                //如果该点为亮点,则将其加入拟合点;

                //cout<<"k1="<<k1<<",i="<<i<<",j="<<j<<endl;

               

                k1= k1+(long longint)i*i*i*i; //k1==2x^4

                //cout<<"k1="<<k1<<endl;

                //getchar();

                k2= k2 + (long longint)i*i*i;   //k2==2x^3

                k3= k3 + (long longint)i*i;     //k3==2x^2

                k4= k4 + (long longint)i*i*j;   //k4==-2y*x^2

                k5= k5 + (long longint)i*i*i;

                k6= k6 + (long longint)i*i;

                k7= k7 + (long longint)i;

                k8= k8 + (long longint)i*j;

                k9= k9 + (long longint)i*i;

                k10= k10 + (long longint)i;

                k11= k11 + 1;

                k12= k12 + (long longint)j;

            }

        }

    }

 

    //根据偏导求出a,b,c的值;

 

    long long int t1 = k1/k3-k5/k7;

    long long int t2 = k2/k3-k6/k7;

    long long int t3 = k4/k3-k8/k7;

    long long int t4 = k1/k3-k9/k11;

    long long int t5 = k2/k3-k10/k11;

    long long int t6 = k4/k3-k12/k11;

 

 

    a = (double)(t3*t5-t6*t2)/(t1*t5-t4*t2);

    b = (double)(t3-t1*a)/t2;

    double b1=(double)(t6-t4*a)/t5;

    c = (double)(k4-k1*a-k2*b)/k3;

 

    int bottom = image_threshold->height;

    int top = 0;

    for(int i =0;i<image_threshold->height;i++)

        for(int j =0;j<image_threshold->width;j++)

        {

            if(cvGetReal2D(image_threshold,i,j) == 255)

            {

                if(i>top)

                    top= i;

                if(i<bottom)

                    bottom= i;

            }

        }

 

    for (inti=bottom;i<top;i++)

    {

        int temp_y = a*i*i+b*i+c;

        if(temp_y<src->width&&temp_y>0)

            cvCircle(src,cvPoint(temp_y,i),1,cvScalar(255,255,255));

            //cvSet2D(src,i,temp_y,cvScalar(255,255,255));

        cvShowImage("二次曲线拟合",src);

    }

 }


拟合效果如下:


2017-04-21 11:45:07 guduruyu 阅读数 56259

在数据处理和绘图中,我们通常会遇到直线或曲线的拟合问题,python中scipy模块的子模块optimize中提供了一个专门用于曲线拟合的函数curve_fit()。

下面通过示例来说明一下如何使用curve_fit()进行直线和曲线的拟合与绘制。

代码如下:

import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize

#直线方程函数
def f_1(x, A, B):
    return A*x + B

#二次曲线方程
def f_2(x, A, B, C):
    return A*x*x + B*x + C

#三次曲线方程
def f_3(x, A, B, C, D):
    return A*x*x*x + B*x*x + C*x + D

def plot_test():

    plt.figure()

    #拟合点
    x0 = [1, 2, 3, 4, 5]
    y0 = [1, 3, 8, 18, 36]

    #绘制散点
    plt.scatter(x0[:], y0[:], 25, "red")

    #直线拟合与绘制
    A1, B1 = optimize.curve_fit(f_1, x0, y0)[0]
    x1 = np.arange(0, 6, 0.01)
    y1 = A1*x1 + B1
    plt.plot(x1, y1, "blue")

    #二次曲线拟合与绘制
    A2, B2, C2 = optimize.curve_fit(f_2, x0, y0)[0]
    x2 = np.arange(0, 6, 0.01)
    y2 = A2*x2*x2 + B2*x2 + C2 
    plt.plot(x2, y2, "green")

    #三次曲线拟合与绘制
    A3, B3, C3, D3= optimize.curve_fit(f_3, x0, y0)[0]
    x3 = np.arange(0, 6, 0.01)
    y3 = A3*x3*x3*x3 + B3*x3*x3 + C3*x3 + D3 
    plt.plot(x3, y3, "purple")

    plt.title("test")
    plt.xlabel('x')
    plt.ylabel('y')

    plt.show()

    return


拟合和绘制解果如下:



当然,curve_fit()函数不仅可以用于直线、二次曲线、三次曲线的拟合和绘制,仿照代码中的形式,可以适用于任意形式的曲线的拟合和绘制,只要定义好合适的曲线方程即可。

如高斯曲线拟合,曲线函数形式如下:

def f_gauss(x, A, B, C, sigma):
    return A*np.exp(-(x-B)**2/(2*sigma**2)) + C

2017.04.21

2017-06-05 13:17:56 guduruyu 阅读数 32279

最小二乘法多项式曲线拟合,是常见的曲线拟合方法,有着广泛的应用,这里在借鉴最小二乘多项式曲线拟合原理与实现的原理的基础上,介绍如何在OpenCV来实现基于最小二乘的多项式曲线拟合。


概念

最小二乘法多项式曲线拟合,根据给定的m个点,并不要求这条曲线精确地经过这些点,而是曲线y=f(x)的近似曲线y= φ(x)。

原理

给定数据点pi(xi,yi),其中i=1,2,…,m。求近似曲线y= φ(x)。并且使得近似曲线与y=f(x)的偏差最小。近似曲线在点pi处的偏差δi= φ(xi)-y,i=1,2,...,m。 


常见的曲线拟合方法:

1.使偏差绝对值之和最小

     

2.使偏差绝对值最大的最小

     

3.使偏差平方和最小

     

按偏差平方和最小的原则选取拟合曲线,并且采取二项式方程为拟合曲线的方法,称为最小二乘法。


推导过程:

1. 设拟合多项式为:

          

2.各点到这条曲线的距离之和,即偏差平方和如下:

          

3.为了求得符合条件的a值,对等式右边求ai偏导数,因而我们得到了: 

          

          

                         .......

          

4.将等式左边进行一下化简,然后应该可以得到下面的等式:

          

          

                     .......

          


5.把这些等式表示成矩阵的形式,就可以得到下面的矩阵:

          

6.即X*A=Y。


我们只要解出这个线性方程,即可求得拟合曲线多项式的系数矩阵。而在opencv中,有一个专门用于求解线性方程的函数,即cv::solve(),具体调用形式如下:

int cv::solve(
	cv::InputArray X, // 左边矩阵X, nxn
	cv::InputArray Y, // 右边矩阵Y,nx1
	cv::OutputArray A, // 结果,系数矩阵A,nx1
	int method = cv::DECOMP_LU // 估算方法
);


我们只需要按照上述原理,构造出矩阵X和Y,即可调用该函数,计算出多项式的系数矩阵A。

opencv中支持的估算方法如下图所示:



实现如下:


bool polynomial_curve_fit(std::vector<cv::Point>& key_point, int n, cv::Mat& A)
{
	//Number of key points
	int N = key_point.size();

	//构造矩阵X
	cv::Mat X = cv::Mat::zeros(n + 1, n + 1, CV_64FC1);
	for (int i = 0; i < n + 1; i++)
	{
		for (int j = 0; j < n + 1; j++)
		{
			for (int k = 0; k < N; k++)
			{
				X.at<double>(i, j) = X.at<double>(i, j) +
					std::pow(key_point[k].x, i + j);
			}
		}
	}

	//构造矩阵Y
	cv::Mat Y = cv::Mat::zeros(n + 1, 1, CV_64FC1);
	for (int i = 0; i < n + 1; i++)
	{
		for (int k = 0; k < N; k++)
		{
			Y.at<double>(i, 0) = Y.at<double>(i, 0) +
				std::pow(key_point[k].x, i) * key_point[k].y;
		}
	}

	A = cv::Mat::zeros(n + 1, 1, CV_64FC1);
	//求解矩阵A
	cv::solve(X, Y, A, cv::DECOMP_LU);
	return true;
}


测试代码如下:

int main()
{
	//创建用于绘制的深蓝色背景图像
	cv::Mat image = cv::Mat::zeros(480, 640, CV_8UC3);
	image.setTo(cv::Scalar(100, 0, 0));

	//输入拟合点  
	std::vector<cv::Point> points;
	points.push_back(cv::Point(100., 58.));
	points.push_back(cv::Point(150., 70.));
	points.push_back(cv::Point(200., 90.));
	points.push_back(cv::Point(252., 140.));
	points.push_back(cv::Point(300., 220.));
	points.push_back(cv::Point(350., 400.));

	//将拟合点绘制到空白图上  
	for (int i = 0; i < points.size(); i++)
	{
		cv::circle(image, points[i], 5, cv::Scalar(0, 0, 255), 2, 8, 0);
	}

	//绘制折线
	cv::polylines(image, points, false, cv::Scalar(0, 255, 0), 1, 8, 0);

	cv::Mat A;

	polynomial_curve_fit(points, 3, A);
	std::cout << "A = " << A << std::endl;

	std::vector<cv::Point> points_fitted;

	for (int x = 0; x < 400; x++)
	{
		double y = A.at<double>(0, 0) + A.at<double>(1, 0) * x +
			A.at<double>(2, 0)*std::pow(x, 2) + A.at<double>(3, 0)*std::pow(x, 3);

		points_fitted.push_back(cv::Point(x, y));
	}
	cv::polylines(image, points_fitted, false, cv::Scalar(0, 255, 255), 1, 8, 0);

	cv::imshow("image", image);

	cv::waitKey(0);
	return 0;
}


绘制结果:



2017.06.05




2018-06-12 10:14:09 weixin_33994429 阅读数 261

OpenCV曲线拟合与圆拟合

使用OpenCV做图像处理与分析的时候,经常会遇到需要进行曲线拟合与圆拟合的场景,很多OpenCV开发者对此却是一筹莫展,其实OpenCV中是有现成的函数来实现圆拟合与直线拟合的,而且还会告诉你拟合的圆的半径是多少,简直是超级方便,另外一个常用到的场景就是曲线拟合,常见的是基于多项式拟合,可以根据设定的多项式幂次生成多项式方程,然后根据方程进行一系列的点生成,形成完整的曲线,这个车道线检测,轮廓曲线拟合等场景下特别有用。下面就通过两个简单的例子来分别学习一下曲线拟合与圆拟合的应用。

一:曲线拟合与应用

基于Numpy包的polyfit函数实现,其支持的三个参数分别是x点集合、y点集合,以及多项式的幂次。得到多项式方程以后,就可以完整拟合曲线,图中有如下四个点:
OpenCV曲线拟合与圆拟合

调用polyfit生成的二阶多项式如下:
OpenCV曲线拟合与圆拟合

拟合结果如下:
OpenCV曲线拟合与圆拟合

使用三阶多项式拟合,调用polyfit生成的多项式方程如下:
OpenCV曲线拟合与圆拟合

生成的拟合曲线如下:
OpenCV曲线拟合与圆拟合

使用polyfit进行曲线拟合时候需要注意的是,多项式的幂次最大是数据点数目N - 1幂次多项式,比如有4个点,最多生成3阶多项式拟合。上述演示的完整代码实现如下:

def circle_fitness_demo():
    image = np.zeros((400, 400, 3), dtype=np.uint8)
    x = np.array([30, 50, 100, 120])
    y = np.array([100, 150, 240, 200])
    for i in range(len(x)):
        cv.circle(image, (x[i], y[i]), 3, (255, 0, 0), -1, 8, 0)
    cv.imwrite("D:/curve.png", image)

    poly = np.poly1d(np.polyfit(x, y, 3))
    print(poly)
    for t in range(30, 250, 1):
        y_ = np.int(poly(t))
        cv.circle(image, (t, y_), 1, (0, 0, 255), 1, 8, 0)
    cv.imshow("fit curve", image)
    cv.imwrite("D:/fitcurve.png", image)

二:圆拟合与应用

圆的拟合是基于轮廓发现的结果,对发现的近似圆的轮廓,通过圆拟合可以得到比较好的显示效果,轮廓发现与拟合的API分别为findContours与fitEllipse,
有图像如下:
OpenCV曲线拟合与圆拟合

使用轮廓发现与圆拟合处理结果如下:
OpenCV曲线拟合与圆拟合

红色表示拟合的圆,蓝色是圆的中心位置
上述完整的演示代码如下:

def circle_fitness_demo():
    src = cv.imread("D:/javaopencv/c2.png")
    cv.imshow("input", src)
    src = cv.GaussianBlur(src, (3, 3), 0)
    gray = cv.cvtColor(src, cv.COLOR_BGR2GRAY)
    ret, binary = cv.threshold(gray, 0, 255, cv.THRESH_BINARY | cv.THRESH_OTSU)
    cv.imshow("binary", binary)

    image, contours, hierachy = cv.findContours(binary, cv.RETR_EXTERNAL, cv.CHAIN_APPROX_SIMPLE)
    for i in range(len(contours)):
        rrt = cv.fitEllipse(contours[i])
        cv.ellipse(src, rrt, (0, 0, 255), 2, cv.LINE_AA)
        x, y = rrt[0]
        cv.circle(src, (np.int(x), np.int(y)), 4, (255, 0, 0), -1, 8, 0)
    cv.imshow("fit circle", src)
    cv.imwrite("D:/fitcircle.png", src)

学习OpenCV Python 系列课程
OpenCV Python零基础入门视频教程
OpenCV Python图像处理进阶视频教程


吾心信其可行,则移山填海之难,终有成功之日;
吾心信其不可行,则反掌折枝之易,亦无收效之期也


转载于:https://blog.51cto.com/gloomyfish/2128313

matlab曲线拟合

阅读数 19778

没有更多推荐了,返回首页