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