最小二乘高手请进~~不怕大家烦，我又来了......

[月小夏]

请帮忙看看法方程运算那里公式有没有问题~~程序整体可以运行，但是就是结果不大对

因为太长，我就分开来帖了

#include <stdlib.h>
#include <malloc.h>
#include <stdio.h>
#include <math.h>


/*==================================================*/
/* 高斯解方程子程序 */
/* a[n1][n1]*x[n1][1]=b[n1][1] AX=B */
/*==================================================*/
int gs(float *a,float *b,int n1,float *x)
{int *js,l,k,i,j,is,p,q;
float d,t;
js=(int *)malloc(n1*sizeof(int));
l=1;
for(k=0;k<=n1-2;k++)
{d=0.0;
for(i=k;i<=n1-1;i++)
for(j=k;j<=n1-1;j++)
{t=(float)fabs(*(a+i*n1+j));
if(t>d) {d=t;js[k]=j;is=i;}
}
if(d+1.0==1.0) l=0;
else
{if(js[k]!=k)
for(i=0;i<=n1-1;i++)
{p=i*n1+k;q=i*n1+js[k];
t=*(a+p);*(a+p)=*(a+q);*(a+q)=t;
}
if(is!=k)
{for(j=k;j<=n1-1;j++)
{p=k*n1+j;q=is*n1+j;
t=*(a+p);*(a+p)=*(a+q);*(a+q)=t;
}
t=*(b+k);*(b+k)=*(b+is);*(b+is)=t;
}
}
if(l==0)
return(0);
d=*(a+k*n1+k);
for(j=k+1;j<=n1-1;j++)
{p=k*n1+j;*(a+p)=*(a+p)/d;}
*(b+k)=*(b+k)/d;
for(i=k+1;i<=n1-1;i++)
{for(j=k+1;j<=n1-1;j++)
{p=i*n1+j;
*(a+p)=*(a+p)-*(a+i*n1+k)*(*(a+k*n1+j));
}
*(b+i)=*(b+i)-*(a+i*n1+k)*(*(b+k));
}
}
d=*(a+(n1-1)*n1+n1-1);
if(fabs(d)+1.0==1.0)
{free(js);
return(0);
}
*(x+n1-1)=*(b+n1-1)/d;
for(i=n1-2;i>=0;i--)
{t=0.0;
for(j=i+1;j<=n1-1;j++)
t=t+*(a+i*n1+j)*(*(x+j));
*(x+i)=*(b+i)-t;
}
js[n1-1]=n1-1;
for(k=n1-1;k>=0;k--)
if(js[k]!=k)
{t=*(x+k);*(x+k)=*(x+js[k]);*(x+js[k])=t;}
free(js);
return(1); //有解
}


这个子程序应该没错，主程序我贴在下面

[月小夏]

void main()
{
// char gx[25],gy[25];
int m=5; //匹配窗口大小m*m=5*5
double min1=0.005;
FILE *fp;
int i,j,k=0; //控制循环变量
int p,q; //控制循环变量
int i1,i2,j1,j2,i3,j3,i4,j4/*,ii1,jj1*//*,ii2*/;

int Height=10; //图像的高
int Width =10; //图像的宽
long size=Height*Width;

unsigned char *Left; //指向左图的指针
unsigned char *Right; //指向右图的指针
Left= (unsigned char *)calloc(size,1);
Right= (unsigned char *)calloc(size,1);

unsigned char Limage[100],Rimage[100]; //左图和右图存放的数组

int k1,k2,k3,k4; //匹配起点(k3,k1),匹配中心(X01,Y01) k1==i k3==j
int X01,Y01,X02,Y02; //搜索起点(k4,k2),搜索中心(X02,Y02) k2==i k4==j

int *Xr,*Yr ; //Xr是X 坐标改正，Yr是Y坐标改正
// Xr = (int *)calloc(m*m,1);
// Yr = (int *)calloc(m*m,1);
Xr = (int *)calloc(m*m,sizeof(int));
Yr = (int *)calloc(m*m,sizeof(int));


unsigned char *F,*G,*GO,*G1;
F =(unsigned char*)calloc(m*m,1); //左匹配窗口
G =(unsigned char*)calloc(4*m*m,1); //右缓存窗口
GO = (unsigned char*)calloc(m*m,1); //右搜索窗口（经过辐射畸变改正的）
G1 = (unsigned char*)calloc(m*m,1); //直接从右缓存窗口读取的小块，辐射改正的中间变量GO[]=h0+h1*G1[]

int Gx,Gy; //差分运算代替求偏导

float A[9],f1[8][8],g1[8][1],x[8][1]; //f1即为C'C; g1即为C'L; x是方程的解:x=[da0,da1,da2,db0,db1,db2,dh0,dh1]
int solut; //高斯方程是否有解
float a[8]={0,1,0,0,0,1,0,1}; //a0,a1,a2,b0,b1,b2,h0,h1

float Xt=0.0,Yt=0.0,Xs=0.0,Ys=0.0; //Xt,Yt计算最佳匹配点位的中间变量，Xs,Ys是同名点坐标
float Xt1=0.0,Xt2=0.0,Yt1=0.0,Yt2=0.0;

int shiftx,shifty;

//读左图全图并存放至数组Limage
fp = fopen("IMG00.bmp","rb");
if (fp==NULL)
{
printf("open file IMG0.raw error!\n");
exit(0);
}
fread(Left,size,1,fp);
fclose(fp); fp = NULL;
for(i=0;i<Height;i++)
for(j=0;j<Width;j++)
Limage[i*Width+j]=Left[i*Width+j];


//读右图全图并存放至数组Rimage
fp = fopen("IMG11.bmp","rb");
if (fp==NULL)
{
printf("open file IMG1.raw error!\n");
exit(0);
}
fread(Right,size,1,fp);
fclose(fp); fp = NULL;
for(i=0;i<Height;i++)
for(j=0;j<Width;j++)
Rimage[i*Width+j]=Right[i*Width+j];

j=0;
while (j<Height/m) //如果高度方向上，左窗口未越界则窗口向下移动m行
{
i=0;
while(i<Width/m) //如果宽度方向上，左窗口未越界则窗口向右移动m列
{
Left = &(Limage[j*m*Width+i*m]); //left指针指向左窗口左上角的那个点，读取m*m大小的数组作为左窗口放在F数组中
for(i1=0;i1<m;i1++)
{
for(j1=0;j1<m;j1++)
F[i1*m+j1]= Left[i1*Width+j1]; //左窗口作为目标图像 !!y*wid+x!!
}
X01=m/2; Y01=m/2; //左窗口匹配中心为(m/2,m/2)，此处开始使用相对坐标
k1 =Y01-m/2; k3 =X01-m/2; //匹配起点(k3,k1),匹配中心(X01,Y01)


//读右窗口，同左窗口方法一样，但是每次读取2m*2m大小的块作为缓存
Right = &(Rimage[j*2*m*Width+i*2*m]);
for(i1=0;i1<2*m;i1++)
for(j1=0;j1<2*m;j1++)
G[i1*2*m+j1]= Right[i1*Width+j1]; //搜索缓存图像

//定搜索窗口（每次从缓存里读出m*m的一块）
shifty=0;
while(shifty<5) //取数字5是未防止右窗口越界
{
shiftx=0;
while(shiftx<5)
{
for(i1=0;i1<m;i1++)//从搜索缓存里读出 G1 是搜索窗口
for(j1=0;j1<m;j1++)
G1[i1*m+j1]=G[(i1+shifty)*2*m+j1+shiftx];

//计算搜索窗口中心位置，使用相对坐标
X02=m/2; Y02=m/2;
k2 =Y02-m/2; k4 =X02-m/2;


//几何改正：根据几何变形改正参数a0,a1,a2,b0,b1,b2将左方影像窗口坐标变化到右方影像窗口
//xr=a0+a1*xL+a2*yL
//yr=b0+b1*xL+b2*yL
for(i1=0;i1<m;i1++)
{
for(j1=0;j1<m;j1++)
{
Xr[i1*m+j1]=(int)(a[0]+a[1]*k3+a[2]*k1);
Yr[i1*m+j1]=(int)(a[3]+a[4]*k3+a[5]*k1);
k3++;
}
k1++;
}
k1 =Y01-m/2; k3 =X01-m/2; //上面使用了k1和k3，故重新初始化一下匹配起点(k3,k1),匹配中心(X01,Y01)

//右方影像窗口进行辐射畸变改正：GO[]=h0+h1*G1[]
for(i2=0;i2<m;i2++)
for(j2=0;j2<m;j2++)
{
k4=Xr[i2*m+j2];
k2=Yr[i2*m+j2];
GO[i2*m+j2]=(unsigned char)(a[6]+a[7]*G1[k2*m+k4]);
}
k2 =Y02-m/2; k4 =X02-m/2;
//初始化：f1即为C'C; g1即为C'L; x是方程的解:x=[da0,da1,da2,db0,db1,db2,dh0,dh1]
for(i2=0;i2<8;i2++)
{
g1[i2][0]=0.;
for(j2=0;j2<8;j2++)
f1[i2][j2]=0.;
}

//判断是否要重解变形参数，即计算变形参数改正值x[]=da0,da1,da2,db0,db1,db2,dh0,dh1
//（判断条件：参数改正值的绝对值小于某个阈值min1）
//重解(min1>.005)：则进行组法方程并解法方程
//不重解(min1<0.005)：计算最佳匹配点位，给出同名点坐标，程序结束

while(min1>=0.005)
{
k=0;
//******法化组成法方程******//
for(i3=0;i3<m;i3++)
{
for(j3=0;j3<m;j3++)
{
i1=(i3-1<0)?i3:i3-1; //防止差分越界，以后无用
i2=(i3+1>m-1)?i3:i3+1;
j1=(j3-1<0)?j3:j3-1;
j2=(j3+1>m-1)?j3:j3+1;
Gx=(GO[i3*m+j2]-GO[i3*m+j1])/2; //差分运算,代替求偏导
Gy=(GO[i2*m+j3]-GO[i1*m+j3])/2;

//******** 矩阵 C‘L*********//
A[0]=-a[7]*Gx; //A[0]=-h1*g'x =-h1*c
A[1]=-a[7]*Gx*Xr[k]; //A[1]=-h1*g'x*x=-h1*c
A[2]=-a[7]*Gx*Yr[k]; //A[2]=-h1*g'x*y=-h1*c
A[3]=-a[7]*Gy; //A[3]=-h1*g'y =-h1*c
A[4]=-a[7]*Gy*Xr[k]; //A[4]=-h1*g'y*x=-h1*c
A[5]=-a[7]*Gy*Yr[k]; //A[5]=-h1*g'y*y=-h1*c
A[6]=-1.0; //A[6]=-h1*1 =-h1*c
A[7]=-GO[i3*m+j3]*1.0f;
A[8]=(float)(F[i3*m+j3]-GO[i3*m+j3]);
k=k+1;

for(p=0;p<8;p++)
{
g1[p][0]=g1[p][0]+A[p]*A[8]; //C'L
for(q=p;q<8;q++) //计算矩阵C'C(step 1)
f1[p][q]=f1[p][q]+A[p]*A[q];
}
//********计算 矩阵 C'C(step 2)*********//
for(p=0;p<8;p++)
for(q=p;q<8;q++)
f1[q][p]=f1[p][q];
}

}

solut=gs(*f1,*g1,8,*x); /* 高斯求解改正数值 */
if(solut) //有解
{
for(i4=0;i4<8;i4++)
a[i4]=a[i4]+x[i4][0];
//对a进行改正
min1=0;
for(i4=0;i4<8;i4++)
if(fabs(x[i4][0])>min1) min1=fabs(x[i4][0]); //求x最大值
}
if(Xr[m*m/2]>= 0.5) X02=X02+(int)(Xr[m*m/2]+0.5);
if(Xr[m*m/2]<=-0.5) X02=X02+(int)(Xr[m*m/2]-0.5);
if(Yr[m*m/2]>= 0.5) Y02=Y02+(int)(Yr[m*m/2]+0.5);
if(Yr[m*m/2]<=-0.5) Y02=Y02+(int)(Yr[m*m/2]-0.5);
}

shiftx++;

}
shifty++;

}


for(i4=0;i4<m;i4++)
{
for(j4=0;j4<m;j4++)
{
i1=(i4-1<0)?i4:i4-1; //防止差分越界，以后无用
i2=(i4+1>m-1)?i4:i4+1;
j1=(j4-1<0)?j4:j4-1;
j2=(j4+1>m-1)?j4:j4+1;
Gx=(F[i4*m+j2]-F[i4*m+j1])/2; //差分运算
Gy=(F[i2*m+j4]-F[i1*m+j4])/2;

Xt1=Xt1+(float)(j4*Gx*Gx);
Xt2=Xt2+(float)(Gx*Gx);
Yt1=Yt1+(float)(i4*Gy*Gy);
Yt2=Yt2+(float)(Gy*Gy);


}

}


Xt=Xt1/Xt2; //对左窗口中心位置进行改正
Yt=Yt1/Yt2;


//*************print****************//
if ((fp=fopen("result.dat","a"))==NULL)
{
printf("Can not open the file result.dat\n");
exit(0);
}
//相对坐标转绝对坐标
float xL,yL;
xL=i*m+Xt;
yL=j*m+Yt;
Xs=a[0]+a[1]*xL+a[2]*yL; //求出同名点坐标
Ys=a[3]+a[4]*xL+a[5]*yL;

fprintf(fp,"(%f %f)\t\t(%f %f)\n",xL,yL,Xs,Ys);
fclose(fp);
fp = NULL;
//************print*****************//


i++;

}
j++;

}
}

[月小夏]

不好意思，c的标准格式我不会调，如果看着很麻烦，希望可以不怕再麻烦一点，粘贴到vc里面看，~(@^_^@)~

[lingyin55]

up!!

[月小夏]

睡了一觉就发现找不到了.......

[YFLK]

解线性方程C++有现成的模块，建议楼主采用已有的模块，这样可减少调试的时间和编程的难度。可靠性也有保证

[月小夏]

楼上的前辈麻烦说清楚一点，我初学的，从两个礼拜前才接触的vc.....呵呵，很多地方都不知道~~~

[月小夏]

不能沉啊不能沉 %>_<%

[杰克包子]

头很晕 帮顶

[月小夏]

O(∩_∩)O谢谢，好几天了没人理，也不知道是什么原因啊.......⊙﹏⊙

[chiwa737]

看着累啊， 你先跟matlab里面试试然后转c代码呢？至少说明算法思想没有错误

[linglongyouzhi]

如果是作业或者研究算法的话建议就照10楼的去做
如果是要完成的工作的话建议用用libgsl很好的开源数值计算库

[月小夏]

MATLAB不熟....比c还外行.....不过谢谢楼上们的建议了！~

[月小夏]

upupup

[sneaper]

好复杂 lz用最小二乘算法解决什么问题呢啊 我也准备研究这玩意

[月小夏]

是进行两幅图像的最小二乘匹配，找出所有同名点的坐标

[xiewhenwe]

SO LONG

[月小夏]

啊....卡住一个月了，师兄也没有帮我解决，有点崩溃了......