www.pudn.com > least-square.rar > Unit1.cpp, change:2007-09-07,size:10109b


//--------------------------------------------------------------------------- 
#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()); 
    } 
    } 
//---------------------------------------------------------------------------