2016-01-02 12:34:16 GarfieldEr007 阅读数 7543

使用Matlab进行拟合是图像处理中线条变换的一个重点内容,本文将详解Matlab中的直线拟合和曲线拟合用法。

关键函数:

fittype

Fit type for curve and surface fitting

Syntax

ffun = fittype(libname)
ffun = fittype(expr)
ffun = fittype({expr1,...,exprn})
ffun = fittype(expr, Name, Value,...)
ffun= fittype({expr1,...,exprn}, Name, Value,...)

/***********************************线性拟合***********************************/

线性拟合公式:

coeff1 * term1 + coeff2 * term2 + coeff3 * term3 + ...
其中,coefficient是系数,term都是x的一次项。

线性拟合Example:

Example1: y=kx+b;

法1:

[csharp] view plaincopy
  1. x=[1,1.5,2,2.5,3];y=[0.9,1.7,2.2,2.6,3];  
  2. p=polyfit(x,y,1);  
  3. x1=linspace(min(x),max(x));  
  4. y1=polyval(p,x1);  
  5. plot(x,y,'*',x1,y1);  
结果:p =    1.0200    0.0400

即y=1.0200 *x+ 0.0400

法2:

[csharp] view plaincopy
  1. x=[1;1.5;2;2.5;3];y=[0.9;1.7;2.2;2.6;3];  
  2. p=fittype('poly1')  
  3. f=fit(x,y,p)  
  4. plot(f,x,y);  
运行结果:

[csharp] view plaincopy
  1.  x=[1;1.5;2;2.5;3];y=[0.9;1.7;2.2;2.6;3];  
  2. p=fittype('poly1')  
  3. f=fit(x,y,p)  
  4. plot(f,x,y);  
  5.   
  6. p =   
  7.   
  8.      Linear model Poly1:  
  9.      p(p1,p2,x) = p1*x + p2  
  10.   
  11. f =   
  12.   
  13.      Linear model Poly1:  
  14.      f(x) = p1*x + p2  
  15.      Coefficients (with 95% confidence bounds):  
  16.        p1 =        1.02  (0.7192, 1.321)  
  17.        p2 =        0.04  (-0.5981, 0.6781)  

Example2:y=a*x + b*sin(x) + c

法1:

[csharp] view plaincopy
  1. x=[1;1.5;2;2.5;3];y=[0.9;1.7;2.2;2.6;3];  
  2. EXPR = {'x','sin(x)','1'};  
  3. p=fittype(EXPR)  
  4. f=fit(x,y,p)  
  5. plot(f,x,y);  

运行结果:

[csharp] view plaincopy
  1.  x=[1;1.5;2;2.5;3];y=[0.9;1.7;2.2;2.6;3];  
  2. EXPR = {'x','sin(x)','1'};  
  3. p=fittype(EXPR)  
  4. f=fit(x,y,p)  
  5. plot(f,x,y);  
  6.   
  7. p =   
  8.   
  9.      Linear model:  
  10.      p(a,b,c,x) = a*x + b*sin(x) + c  
  11.   
  12. f =   
  13.   
  14.      Linear model:  
  15.      f(x) = a*x + b*sin(x) + c  
  16.      Coefficients (with 95% confidence bounds):  
  17.        a =       1.249  (0.9856, 1.512)  
  18.        b =      0.6357  (0.03185, 1.24)  
  19.        c =     -0.8611  (-1.773, 0.05094)  

法2:
[csharp] view plaincopy
  1. x=[1;1.5;2;2.5;3];y=[0.9;1.7;2.2;2.6;3];  
  2.  p=fittype('a*x+b*sin(x)+c','independent','x')  
  3. f=fit(x,y,p)  
  4. plot(f,x,y);  
运行结果:
[csharp] view plaincopy
  1. x=[1;1.5;2;2.5;3];y=[0.9;1.7;2.2;2.6;3];  
  2.  p=fittype('a*x+b*sin(x)+c','independent','x')  
  3. f=fit(x,y,p)  
  4. plot(f,x,y);  
  5.   
  6. p =   
  7.   
  8.      General model:  
  9.      p(a,b,c,x) = a*x+b*sin(x)+c  
  10. Warning: Start point not provided, choosing random start  
  11. point.   
  12. > In fit>iCreateWarningFunction/nThrowWarning at 738  
  13.   In fit>iFit at 320  
  14.   In fit at 109   
  15.   
  16. f =   
  17.   
  18.      General model:  
  19.      f(x) = a*x+b*sin(x)+c  
  20.      Coefficients (with 95% confidence bounds):  
  21.        a =       1.249  (0.9856, 1.512)  
  22.        b =      0.6357  (0.03185, 1.24)  
  23.        c =     -0.8611  (-1.773, 0.05094)  


/***********************************非线性拟合***********************************/

Example:y=a*x^2+b*x+c

法1:

  1. x=[1;1.5;2;2.5;3];y=[0.9;1.7;2.2;2.6;3];  
  2.  p=fittype('a*x.^2+b*x+c','independent','x')  
  3. f=fit(x,y,p)  
  4. plot(f,x,y);  

运行结果:

[csharp] view plaincopy
  1. p =   
  2.   
  3.      General model:  
  4.      p(a,b,c,x) = a*x.^2+b*x+c  
  5. Warning: Start point not provided, choosing random start  
  6. point.   
  7. > In fit>iCreateWarningFunction/nThrowWarning at 738  
  8.   In fit>iFit at 320  
  9.   In fit at 109   
  10.   
  11. f =   
  12.   
  13.      General model:  
  14.      f(x) = a*x.^2+b*x+c  
  15.      Coefficients (with 95% confidence bounds):  
  16.        a =     -0.2571  (-0.5681, 0.05386)  
  17.        b =       2.049  (0.791, 3.306)  
  18.        c =       -0.86  (-2.016, 0.2964)  



法2:

[csharp] view plaincopy
  1. x=[1;1.5;2;2.5;3];y=[0.9;1.7;2.2;2.6;3];  
  2. %use c=0;  
  3. c=0;  
  4. p1=fittype(@(a,b,x) a*x.^2+b*x+c)  
  5. f1=fit(x,y,p1)  
  6. %use c=1;  
  7. c=1;  
  8. p2=fittype(@(a,b,x) a*x.^2+b*x+c)  
  9. f2=fit(x,y,p2)  
  10. %predict c  
  11. p3=fittype(@(a,b,c,x) a*x.^2+b*x+c)  
  12. f3=fit(x,y,p3)  
  13.   
  14. %show results  
  15. scatter(x,y);%scatter point  
  16. c1=plot(f1,'b:*');%blue  
  17. hold on  
  18. plot(f2,'g:+');%green  
  19. hold on  
  20. plot(f3,'m:*');%purple  
  21. hold off  


from: http://blog.csdn.net/abcjennifer/article/details/7684836

2018-03-19 15:47:39 Sirius_007 阅读数 6298

1.最小二乘拟合

    最简单,但是受离群点的影响比较大,鲁棒性不强。

2.梯度下降法

3.高斯牛顿、列-马算法

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 阅读数 56273

在数据处理和绘图中,我们通常会遇到直线或曲线的拟合问题,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

2018-03-09 20:13:11 qq_34777600 阅读数 7374

数据拟合: 直线拟合--多项式拟合


1.问题概述

在实际问题中,常常需要从一组观察数据

       (xi,yi)  i=1,2,,..,n
去预测函数 y=f(x) 的表达式,从几何角度来说,这个问题就是要由给定的一组数据点(xi,yi)去描绘曲线 y=f(x) 的近似图像。插值方法是处理这类问题的一种数值方法。不过,由于插值曲线要求严格通过所给的每一个数据点,这种限制会保留所给数据的误差。如果个别数据的误差很大,那么插值效果显然是不理想的。
现在面对的问题具有这样的特点:所给数据本事不一定可靠,个别数据的误差甚至可能很大,但给出的数据很多。曲线拟合方法所研究的课题是:从给出的一大堆看上去杂乱无章的数据中找出规律来,就是说,设法构造一条曲线,即所谓拟合曲线,反映所给数据点总的趋势,以消除所给数据的局部波动。

2.理论与方法

假设所给数据点 (xi,yi)  i=1,2,,..,n 的分布大致成一条直线。虽然不能要求所作的拟合直线
                y = a + bx
严格的通过所有数据点(xi,yi)但总希望它尽可能地从所给数据点附近通过,就是说,要求近似地成立
              
这里,数据点的数目通常远远大于待定系数的数目,即N2,因此,拟合直线的构造本质上是个解超定方程组的代数问题。

          
表示按拟合直线y=a+bx求得的近似值,一般来说,它不同于实测值yi,两者之差称作残差。显然,残差的大小是衡量拟合好  坏的重要标志。具体地说,构造拟合曲线可以采用下列三种准则之一:

    使残差的最大绝对值为最小,即

    使残差的绝对值之和为最小,即

          

    使残差的平方和为最小,即

      

分析以上三种准则,前两种提法比较自然,但由于含有绝对值运算,不便于实际应用;基于第三种准则来选取拟合曲线的方法则称作曲线拟合的最小二乘法。


3.算法与设计

1)  直线拟合

对于给定的数据点(xi,yi),i=1,2,...n求作一次式 y=a+bx ,使总误差
为最小。使Q达到极值的参数a,b应满足

即成立

式中,求和表示关于下标i1N求和。


2)  多项式拟合
有时所给数据点用直线拟合并不合适,这时可以考虑用多项式拟合。
对于给定的数据点(xi,yi),i=1,2,...,n求作mm<<n)次多项式

使总误差

为最小。由于Q可以视作关于aj(j = 0,1,...,m)的多元函数,故上述拟合多项式的构造问题可归结为多元函数的极值问题。


即得到关于系数aj的线性方程组通常称作正则方程组。


4.案例分析

Example:

Data set1

Memory Capacity in GBytes

Price in US dollars

2

9.99

4

10.99

8

19.99

16

29.99

1    What's therelationship between memory capacity and cost? Pleas fitting a linear and higher polynomial function model to the data.

由实验结果可以得出在多项式拟合曲线的过程中,不一定是次数越高,结果越准确,结合实际问题与残差值得衡量,可以找到对于具体问题几阶的拟合效果更好。

汇总结果函数图像为:

根据函数图像可以得出表格:

二阶:

Memory Capacity in GBytes

Price in US dollars

2

9.54

4

12.52

8

18.49

16

30.42

三阶:

Memory Capacity in GBytes

Price in US dollars

2

9.06

4

12.61

8

19.17

16

30.1

四阶:

Memory Capacity in GBytes

Price in US dollars

2

9.99

4

11

8

19.98

16

30

通过与原表格的数据进行分析:

不难看出四阶的函数表达式最接近原题。

即:y=13.0376+-2.75*x+0.666667*x^2+-0.0267857*x^3

但是由函数的图像观察到,x在[2,15]的区间内比较符合题意,超出区间后函数的走势开始慢慢趋于不符合实际情况。所以我们考虑在15以后采用二阶的函数。

5.代码

void _2()
{
	//一阶
	double mat[105][2];
	double b[105];
	ifstream ifile;
	ifile.open("e:\\sat.txt");
	for (int i = 0; i < 105; i++)
	{
		mat[i][0] = 1;
		ifile >> mat[i][1];

	//	cout << mat[i][0] << "  " << mat[i][1] << endl;
	}
	for (int i = 0; i < 105; i++)
	{
		ifile >> b[i];

	//	cout << b[i] << endl;
	}
	ifile.close();
/*
	double t = 0, t2 = 0, bi = 0, tb = 0;

	for (int i = 0; i < 105; i++)
	{
		bi += b[i];
		for (int j = 0; j < 2; j++)
		{
			if (j == 1)
			{
				t += mat[i][j];
				t2 += mat[i][j] * mat[i][j];
				tb += mat[i][j] * b[i];
			}
		}
	}

	double x1 = (t2*bi - tb*t) / (105 * t2 - t*t);
	double x2 = (105 * tb - t*bi) / (105 * t2 - t*t);

	//cout << x1 << " " << x2 << endl;
	cout << "拟合方程(关于high_gpa):";
	cout << endl << "y = " << x1 << " + " << x2 << "x" << endl;
	//一阶
*/
	//4阶
	double mat21[105][5];
	ifstream file;
	file.open("e:\\sat.txt");
	for (int i = 0; i < 105; i++)
	{
		mat21[i][0] = 1;
		file >> mat21[i][1];
		//	cout << mat[i][0] << "  " << mat[i][1] << endl;
	}
	for (int i = 0; i < 105; i++) file >> mat21[i][2];
	for (int i = 0; i < 105; i++) file >> mat21[i][3];
	for (int i = 0; i < 105; i++) file >> mat21[i][4];

	for (int i = 0; i < 105; i++)
	{
		file >> b[i];
		//	cout << b[i] << endl;
	}
	file.close();

	double mat22[5][105];
	memset(mat22, 0, sizeof(mat22));
	for (int i = 0; i < 5; i++)
	{
		for (int j = 0; j < 105; j++)
		{
			mat22[i][j] = mat21[j][i];
		}
	}
	double bb2[5][5];
	memset(bb2, 0, sizeof(bb2));
	for (int i = 0; i < 5; i++)
	{
		for (int j = 0; j < 5; j++)
		{
			for (int k = 0; k < 105; k++)
			{
				bb2[i][j] += mat22[i][k] * mat21[k][j];
			}
		}
	}
	double yy2[5];
	memset(yy2, 0, sizeof(yy2));
	for (int i = 0; i < 5; i++)
	{
		for (int k = 0; k < 105; k++)
		{
			yy2[i] += mat22[i][k] * b[k];
		}
	}

	double x21[5], x22[5];
	memset(x21, 0, sizeof(x21));
	memset(x22, 0, sizeof(x22));
	double  temp;
	for (int k = 0; k < 5; k++)
	{
		for (int i = k + 1; i < 5; i++)
		{
			temp = bb2[i][k] /bb2[k][k];
			for (int j = k + 1; j < 5; j++)
			{
				bb2[i][j] -= temp*bb2[k][j];
			}
			yy2[i] -= temp*yy2[k];
		}
	}
	for (int i = 4; i >= 0; i--)
	{
		x21[i] = b[i];
		for (int j = 4; j >= 0; j--)
		{
			if (i != j)
			{
				x21[i] -= bb2[i][j] * x21[j];
			}
		}
		x21[i] /= bb2[i][i];
	}
	cout << endl;
	cout << "拟合方程(关于high_gpa,math_sat,verb_sat,comp_gpa):";
	cout << endl << "y = " << x21[0] << " + " << x21[1] << "*high_gpa" << " + " << x21[2] << "math_sat" << " + " << x21[3] << "verb_sat" << " + " << x21[4] << "comp_gpa"<<endl;
	return;
}

End

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