新建 文本文档1.txt

这个程序可以读取卫星影像外方位元素文件并保存为制定格式

void CDibApiTestView::Onransac3() 
{
	// TODO: Add your command handler code here
	CFileDialog dlg(TRUE,_T("TXT"),_T("*.TXT"),OFN_HIDEREADONLY|OFN_OVERWRITEPROMPT,_T("文本文件(*.TXT)|*.TXT|"));	
    CString wjname;
	if(IDOK==dlg.DoModal())
	wjname.Format("%s",dlg.GetPathName());
    CFile file(wjname,CFile::modeRead);
	CArchive ar(&file,CArchive::load);
	int i,num0=426;//451;//333;
	double r[2][1000],c[2][1000];
	for (i=0;i<num0;i++)
	{
		ar>>c[0][i];
		ar>>r[0][i];
		ar>>c[1][i];
		ar>>r[1][i];
	}
	typedef struct  
	{
		int i[3];//行
		int j[3];//列
	}sanjiaoxing;
	sanjiaoxing sjx;//三角形
		/*产生随机序列*/
	int temp;
	int ind[3];
re:		srand((unsigned)time(NULL)); 
	    for(i=0;i<3;i++)
		{
		   temp = rand();
		   ind[i] = ((double)(temp%1000))/1000*num0;
		}
	for (i=0;i<3;i++)
	{
		sjx.i[i]=r[0][ind[i]];
		sjx.j[i]=c[0][ind[i]];
	}
	double d[3];//三角形三边长
    d[0]=sqrt((sjx.i[0]-sjx.i[1])*(sjx.i[0]-sjx.i[1])+(sjx.j[0]-sjx.j[1])*(sjx.j[0]-sjx.j[1]));
	d[1]=sqrt((sjx.i[0]-sjx.i[2])*(sjx.i[0]-sjx.i[2])+(sjx.j[0]-sjx.j[2])*(sjx.j[0]-sjx.j[2]));
	d[2]=sqrt((sjx.i[1]-sjx.i[2])*(sjx.i[1]-sjx.i[2])+(sjx.j[1]-sjx.j[2])*(sjx.j[1]-sjx.j[2]));
	double mj,midd,midd2,pc[3];
	midd=(d[0]+d[1]+d[2])/2;
	midd2=(d[0]+d[1]+d[2])/3;
	for (i=0;i<3;i++)
	{
		pc[i]=fabs(d[i]-midd2)/midd2;
	}
	mj=sqrt(fabs(midd*(midd-d[0])*(midd-d[1])*(midd-d[2])));
	//锐角三角形条件：三边长不能相差太大

	if (mj<2500||mj>16000||pc[0]>0.4||pc[1]>0.4||pc[2]>0.4)
	{
		goto re;
	}
	//找到了合适面积的三角形，找出三角形内的点，锐角三角形？或圆？
	//取重心，求临近圆域
    double center[2],R;
	center[0]=(c[0][ind[0]]+c[0][ind[1]]+c[0][ind[2]]);
	center[1]=(r[0][ind[0]]+r[0][ind[1]]+r[0][ind[2]]);
	R=(sqrt((c[0][ind[0]]-center[0])*(c[0][ind[0]]-center[0])+(r[0][ind[0]]-center[1])*(r[0][ind[0]]-center[1]))+\
        sqrt((c[0][ind[1]]-center[0])*(c[0][ind[1]]-center[0])+(r[0][ind[1]]-center[1])*(r[0][ind[1]]-center[1]))+\
		sqrt((c[0][ind[2]]-center[0])*(c[0][ind[2]]-center[0])+(r[0][ind[2]]-center[1])*(r[0][ind[2]]-center[1])))/3;
	//统计位于这一圆域内的所有点
	double rhx[200],chx[200];
	int hxdsh=0,hxdh[200];
	for (i=0;i<num0;i++)
	{
		d[0]=sqrt((c[0][i]-center[0])*(c[0][i]-center[0])+(r[0][i]-center[1])*(r[0][i]-center[1]));
		if (d[0]<R)
		{
            rhx[hxdsh]=r[0][i];
			chx[hxdsh]=c[0][i];
			hxdh[hxdsh++]=i;
		}
	}
	if (hxdsh<5)//圆域包含3+5个点
	{
		goto re;
	}
	//计算仿射变换系数
    double x[2][1000],y[2][1000];
	x[0][0]=c[0][ind[0]],y[0][0]=r[0][ind[0]],x[1][0]=c[1][ind[0]],y[1][0]=r[1][ind[0]];
	x[0][1]=c[0][ind[1]],y[0][1]=r[0][ind[1]],x[1][1]=c[1][ind[1]],y[1][1]=r[1][ind[1]];
	x[0][2]=c[0][ind[2]],y[0][2]=r[0][ind[2]],x[1][2]=c[1][ind[2]],y[1][2]=r[1][ind[2]];
	double a[2][5];
	double add[2];
	add[0]=add[1]=0;
	for (i=0;i<3;i++)
	{
        add[0]+=x[1][i]-x[0][i];
		add[1]+=y[1][i]-y[0][i];
	}
	a[0][0]=add[0]/3.0;
	a[1][0]=add[1]/3.0;
	for (i=1;i<3;i++)
	{
		a[0][i]=a[1][i]=0;
	}
	a[0][1]=a[1][2]=1;
	int num=3,xhcsh=0;
	double b[1000],bt[1000],btb[25],btl[5],l[1000],da[5];
	bool fg=0;
	do 
	{
        for (i=0;i<num;i++)
		{
			b[i*3]=1;
		    b[i*3+1]=x[0][i];
			b[i*3+2]=y[0][i];
		}
		for (i=0;i<num;i++)
		{
			if (fg==0)
			{
			    l[i]=x[1][i]-(a[fg][0]+a[fg][1]*x[0][i]+a[fg][2]*y[0][i]);
			}
			else
			{
			    l[i]=y[1][i]-(a[fg][0]+a[fg][1]*x[0][i]+a[fg][2]*y[0][i]);
			}
		}
		convert(b,num,3,bt);
		matrix(bt,3,b,3,num,btb);
		matrix(bt,3,l,1,num,btl);
		Inv(btb,3);
		matrix(btb,3,btl,1,3,da);
		for (i=0;i<3;i++)
		{
			a[fg][i]+=da[i];
		}
		if (fg==0)
		{
			fg=1;
		}
		else if (fg==1)
		{
			fg=0;
		}
		xhcsh++;
	} while(xhcsh<10);
    //统计适合这一仿射变换的点对数，大于一定阈值则认为可以接受（包括了三角形顶点）
	double yzx=5,yzy=5;
	double jsz[2];
	int No[500],No2[500];
	int tj=0;
	for (i=0;i<hxdsh;i++)
	{
        jsz[0]=c[1][hxdh[i]]-(a[0][0]+a[0][1]*chx[i]+a[0][2]*rhx[i]);
		jsz[1]=r[1][hxdh[i]]-(a[1][0]+a[1][1]*chx[i]+a[1][2]*rhx[i]);
		if (fabs(jsz[0])<yzx&&fabs(jsz[1])<yzy)
		{
			No[tj]=i;
			No2[tj++]=hxdh[i];
			if (tj>2)
			{
				break;
			}
		}
	}
	//如果符合这一变换的点数小于阈值，重新寻找三角形，否则，用所有点对计算变换参数，并去除其中余差较大的点
	if (tj<=2)
	{
		goto re;
	}
	else
	{
		num=tj;
	    for (i=0;i<num;i++)
	    {
			x[0][i]=chx[No[i]],y[0][i]=rhx[No[i]],x[1][i]=c[1][No2[i]],y[1][i]=r[1][No2[i]];
            x[0][i]=c[0][No2[i]],y[0][i]=r[0][No2[i]];//检验用
	    }
	    add[0]=add[1]=0;
	    for (i=0;i<num;i++)
		{
            add[0]+=x[1][i]-x[0][i];
		    add[1]+=y[1][i]-y[0][i];
		}
	a[0][0]=add[0]/num;
	a[1][0]=add[1]/num;
	for (i=1;i<4;i++)
	{
		a[0][i]=a[1][i]=0;
	}
	a[0][1]=a[1][2]=1;
	xhcsh=0;
	fg=0;
	do 
	{
        for (i=0;i<num;i++)
		{
			b[i*4]=1;
		    b[i*4+1]=x[0][i];
			b[i*4+2]=y[0][i];
			b[i*4+3]=x[0][i]*y[0][i];
		}
		for (i=0;i<num;i++)
		{
			if (fg==0)
			{
			    l[i]=x[1][i]-(a[fg][0]+a[fg][1]*x[0][i]+a[fg][2]*y[0][i]+a[fg][3]*x[0][i]*y[0][i]);
			}
			else
			{
			    l[i]=y[1][i]-(a[fg][0]+a[fg][1]*x[0][i]+a[fg][2]*y[0][i]+a[fg][3]*x[0][i]*y[0][i]);
			}
		}
		convert(b,num,4,bt);
		matrix(bt,4,b,4,num,btb);
		matrix(bt,4,l,1,num,btl);
		Inv(btb,4);
		matrix(btb,4,btl,1,4,da);
		for (i=0;i<4;i++)
		{
			a[fg][i]+=da[i];
		}
		if (fg==0)
		{
			fg=1;
		}
		else if (fg==1)
		{
			fg=0;
		}
		xhcsh++;
	} while(xhcsh<10);
	}

    //用新的变换参数检验余差大的点
	yzx=5,yzy=5;
    for (i=0;i<hxdsh;i++)
    {
		jsz[0]=c[1][hxdh[i]]-(a[0][0]+a[0][1]*c[0][hxdh[i]]+a[0][2]*r[0][hxdh[i]]+a[0][3]*r[0][hxdh[i]]*c[0][hxdh[i]]);
		jsz[1]=r[1][hxdh[i]]-(a[1][0]+a[1][1]*c[0][hxdh[i]]+a[1][2]*r[0][hxdh[i]]+a[1][3]*r[0][hxdh[i]]*c[0][hxdh[i]]);
		if (fabs(jsz[0])>yzx||fabs(jsz[1])>yzy)
		{
			for (int j=hxdh[i];j<num0;j++)
			{
				r[0][j]=r[0][j+1];
				c[0][j]=c[0][j+1];
				r[1][j]=r[1][j+1];
				c[1][j]=c[1][j+1];
			}
			i--;//warning后面还有i++
			num0--;
		}
    }
	//已完成一个三角形的检查
	//保存这个三角形顶点，下次检验不用该三角形

	for (i=0;i<num0;i++)
    {
     	p1[i].y=r[0][i];
		p1[i].x=c[0][i];
		p2[i].y=r[1][i];
		p2[i].x=c[1][i];
    }
	fsycl==1;
	ppd=num0;
	///////////////////////////////////////////////////////
	CString strout;
	strout.Format("%d",num0);
	MessageBox(strout);    	
}