unit1.~cpp

最小二乘法的曲面拟合方法代码

//---------------------------------------------------------------------------
#include <vcl.h>
#pragma hdrstop
#include "cv.h"
#include "highgui.h"
#include "math.h"
#include <stdio.h>
#include "Unit1.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
#define e 4
#define f 4
int handle;
float a[400];
TForm1 *Form1;
//---------------------------------------------------------------------------
__fastcall TForm1::TForm1(TComponent* Owner)
        : TForm(Owner)
{
}
//---------------------------------------------------------------------------

void __fastcall TForm1::FormCreate(TObject *Sender)
{
    String str1 = " ";
    handle = FileOpen("F:\\毕业论文\\朱暌\\曲面拟合\\result.txt",2);
    int FileSize = FileSeek(handle,0,2);  //gain the size of the file
    FileSeek(handle,0,0);
    if(FileSize != 0)     //clear the file if the file is not null
    {
        for(int i = 0; i < FileSize; i++)
        {
            FileWrite(handle,str1.c_str(),str1.Length());
        }
    }
    FileSeek(handle,0,0);
}
//---------------------------------------------------------------------------
 void hpir2(float x[],float y[],float z[],int n,int m,float a[],int p,int q)
 //求多项式系数；
  {
    int i,j,k,l,kk;
    float apx[20],apy[20],bx[20],by[20],u[20][20];
    float t[20],t1[20],t2[20],xx,yy,d1,d2,g,g1,g2;
    float x2,dd,y1,x1,v[1000];
    d1=1.0*n; apx[0]=0.0;
    for (i=0; i<=n-1; i++)
      apx[0]=apx[0]+x[i];
    apx[0]=apx[0]/d1;
    for (j=0; j<=m-1; j++)
      { v[j]=0.0;
        for (i=0; i<=n-1; i++)
          v[j]=v[j]+z[i*m+j];
        v[j]=v[j]/d1;
      }
    if (p>1)
      { d2=0.0; apx[1]=0.0;
        for (i=0; i<=n-1; i++)
          { g=x[i]-apx[0];
            d2=d2+g*g;
            apx[1]=apx[1]+x[i]*g*g;
          }
        apx[1]=apx[1]/d2;
        bx[1]=d2/d1;
        for (j=0; j<=m-1; j++)
          { v[m+j]=0.0;
            for (i=0; i<=n-1; i++)
              { g=x[i]-apx[0];
                v[m+j]=v[m+j]+z[i*m+j]*g;
              }
            v[m+j]=v[m+j]/d2;
          }
        d1=d2;
      }
    for (k=2; k<=p-1; k++)
      { d2=0.0; apx[k]=0.0;
        for (j=0; j<=m-1; j++) v[k*m+j]=0.0;
        for (i=0; i<=n-1; i++)
          { g1=1.0; g2=x[i]-apx[0];
            for (j=2; j<=k; j++)
              { g=(x[i]-apx[j-1])*g2-bx[j-1]*g1;
                g1=g2; g2=g;
              }
            d2=d2+g*g;
            apx[k]=apx[k]+x[i]*g*g;
            for (j=0; j<=m-1; j++)
              v[k*m+j]=v[k*m+j]+z[i*m+j]*g;
          }
        for (j=0; j<=m-1; j++)
          v[k*m+j]=v[k*m+j]/d2;
        apx[k]=apx[k]/d2;
        bx[k]=d2/d1;
        d1=d2;
      }
    d1=m; apy[0]=0.0;
    for (i=0; i<=m-1; i++)
      apy[0]=apy[0]+y[i];
    apy[0]=apy[0]/d1;
    for (j=0; j<=p-1; j++)
      { u[j][0]=0.0;
        for (i=0; i<=m-1; i++)
	  u[j][0]=u[j][0]+v[j*m+i];
	u[j][0]=u[j][0]/d1;
      }
    if (q>1)
      { d2=0.0; apy[1]=0.0;
        for (i=0; i<=m-1; i++)
          { g=y[i]-apy[0];
            d2=d2+g*g;
            apy[1]=apy[1]+y[i]*g*g;
          }
        apy[1]=apy[1]/d2;
        by[1]=d2/d1;
        for (j=0; j<=p-1; j++)
	  { u[j][1]=0.0;
            for (i=0; i<=m-1; i++)
              { g=y[i]-apy[0];
		u[j][1]=u[j][1]+v[j*m+i]*g;
              }
	    u[j][1]=u[j][1]/d2;
          }
        d1=d2;
      }
    for (k=2; k<=q-1; k++)
      { d2=0.0; apy[k]=0.0;
	for (j=0; j<=p-1; j++) u[j][k]=0.0;
        for (i=0; i<=m-1; i++)
          { g1=1.0;
            g2=y[i]-apy[0];
            for (j=2; j<=k; j++)
              { g=(y[i]-apy[j-1])*g2-by[j-1]*g1;
                g1=g2; g2=g;
              }
            d2=d2+g*g;
            apy[k]=apy[k]+y[i]*g*g;
            for (j=0; j<=p-1; j++)
	      u[j][k]=u[j][k]+v[j*m+i]*g;
          }
        for (j=0; j<=p-1; j++)
	  u[j][k]=u[j][k]/d2;
        apy[k]=apy[k]/d2;
        by[k]=d2/d1;
        d1=d2;
      }
    v[0]=1.0; v[m]=-apy[0]; v[m+1]=1.0;
    for (i=0; i<=p-1; i++)
      for (j=0; j<=q-1; j++)
        a[i*q+j]=0.0;
    for (i=2; i<=q-1; i++)
      { v[i*m+i]=v[(i-1)*m+(i-1)];
        v[i*m+i-1]=-apy[i-1]*v[(i-1)*m+i-1]+v[(i-1)*m+i-2];
        if (i>=3)
          for (k=i-2; k>=1; k--)
            v[i*m+k]=-apy[i-1]*v[(i-1)*m+k]+
                     v[(i-1)*m+k-1]-by[i-1]*v[(i-2)*m+k];
        v[i*m]=-apy[i-1]*v[(i-1)*m]-by[i-1]*v[(i-2)*m];
      }
    for (i=0; i<=p-1; i++)
      { if (i==0) { t[0]=1.0; t1[0]=1.0;}
        else
          { if (i==1)
              { t[0]=-apx[0]; t[1]=1.0;
                t2[0]=t[0]; t2[1]=t[1];
              }
            else
              { t[i]=t2[i-1];
                t[i-1]=-apx[i-1]*t2[i-1]+t2[i-2];
                if (i>=3)
                  for (k=i-2; k>=1; k--)
                    t[k]=-apx[i-1]*t2[k]+t2[k-1]
                         -bx[i-1]*t1[k];
                t[0]=-apx[i-1]*t2[0]-bx[i-1]*t1[0];
                t2[i]=t[i];
                for (k=i-1; k>=0; k--)
                  { t1[k]=t2[k]; t2[k]=t[k];}
              }
          }
        for (j=0; j<=q-1; j++)
          for (k=i; k>=0; k--)
            for (l=j; l>=0; l--)
	      a[k*q+l]=a[k*q+l]+u[i][j]*t[k]*v[j*m+l];
      }
    return;
  }
//---------------------------------------------------------------------------
float jmaxnf(float x[],int n)//分别对x,y求导；
{
    int i,j;
    float w[4]={0.0,0.0,0.0,0.0},y;
    for(int i=0;i<e;i++)
    {
        for(int j=0;j<f;j++)
        {
                w[0]=a[i*f+j]*(float)pow(double(x[0]),double(i))*(float)pow(double(x[1]),double(j));
                w[1]=w[1]+w[0];
                w[2]=w[2]+i*w[0]/x[0];
                w[3]=w[3]+j*w[0]/x[1];
        }
    }
    switch(n)
      { case 0: y=w[1];break;
        case 1: y=w[2]; break;
        case 2: y=w[3]; break;
        default: {}
      }
    return(y);
}
//---------------------------------------------------------------------------
void jmaxn(float x[],float eps,int k,int js[])  //求解二元多次不等式；求极点
{
    int i,j,m,l,jt,il;
    float y[10],b[10];
    float p,z,t,h1,h2,dx,d;
    js[0]=0;
    jt=1;
    h2=0.0;
    while (jt==1)
      { t=0.0; js[0]=js[0]+1;
        for (i=1; i<=2; i++)
          { d=jmaxnf(x,i);
            t=t+fabs(d);
          }
        if (t<eps) jt=0;
        else
          { for (i=0; i<=1; i++)
              { il=5;
                while (il!=0)
                  { j=0; t=x[i]; il=il-1;
                    while (j<=7)
                      { if (j<=2) z=t+j*0.1;
                        else z=h2;
                        x[i]=z;
                        d=jmaxnf(x,i+1);
                        if (fabs(d)+1.0==1.0)
                          { j=10; il=0;}
                        else
                          { h1=d; h2=z;
                            if (j==0)
                              { y[0]=h1; b[0]=h2;}
                            else
                              { y[j]=h1; m=0; l=0;
                                while ((m==0)&&(l<=j-1))
                                    { p=h2-b[l];
                                      if (fabs(p)+1.0==1.0) m=1;
                                      else h2=(h1-y[l])/p;
                                      l=l+1;
                                    }
                                b[j]=h2;
                                if (m!=0) b[j]=1.0e+35;
                                h2=0.0;
                                for (l=j-1; l>=0; l--)
                                  h2=-y[l]/(b[l+1]+h2);
                                h2=h2+b[0];
                              }
                            j=j+1;
                          }
                      }
                    x[i]=h2;
                  }
                x[i]=z;
              }
            if (js[0]==k) jt=0;
          }
      }
    js[1]=1;
    d=jmaxnf(x,0); x[2]=d;
    dx=0.0001; t=x[0];
    x[0]=t+dx; h1=jmaxnf(x,0);
    x[0]=t-dx; h2=jmaxnf(x,0);
    t=h1+h2-2.0*d;
    if (t>0.0) js[1]=-1;
    j=1; jt=1;
    while (jt==1)
      { j=j+1; dx=0.0001; jt=0;
        t=x[j-1];
        x[j-1]=t+dx; h2=jmaxnf(x,0);
        x[j-1]=t-dx; h1=jmaxnf(x,0);
        t=h1+h2-2.0*d;
        if ((t*js[1]<0.0)&&(j<2)) jt=1;
      }
    if (t*js[1]>0.0) js[1]=0;

    return;
  }
//---------------------------------------------------------------------------

void __fastcall TForm1::Button1Click(TObject *Sender)
{
    int i,j,m=0;
    float x[20],y[20],z[400],dt[3];
    float r[3]={0,0,0},eps=1.0e-4;
    int js[2];
    IplImage *src;
    OpenDialog1->FileName = "";
    OpenDialog1->Execute();
    if(OpenDialog1->FileName == "") return;
    src = cvLoadImage(OpenDialog1->FileName.c_str(),0);
    for(j = 0; j < src->width; j++) {
    for(i = 0; i < src->height; i++) {
        z[m]= ((uchar*)(src->imageData + src->widthStep*i))[j];
        m++;
        }
        }
    for (i=0; i<src->width; i++)
    {
    x[i]=i;
    }
    for (i=0; i<src->height; i++)
    {
    y[i]=i;
    }
    hpir2(x,y,z,src->width,src->height,a,e,f);
    for(i=0;i<e;i++)
    {
    for(j=0;j<f;j++)
    {
    String str1 = FloatToStr(a[i*f+j]) + " ";
    FileWrite(handle,str1.c_str(),str1.Length());
    }
    }
    String str7 ="\n";
     FileWrite(handle,str7.c_str(),str7.Length());
    r[0]=1.0;r[1]=1.0;
    jmaxn(r,eps,10,js);
    String str3 = IntToStr(js[0]) + ",";
    FileWrite(handle,str3.c_str(),str3.Length());
    if(js[1]<0)
    {
    String str4 ="min:";
    FileWrite(handle,str4.c_str(),str4.Length());
    }
    else
    {
    String str5 ="max:";
    FileWrite(handle,str5.c_str(),str5.Length());
    }
    for(i=0;i<3;i++)
    {
    String str6 = FloatToStr(r[i]) + ",";
    FileWrite(handle,str6.c_str(),str6.Length());
    }
    }
//---------------------------------------------------------------------------