核线影像制作--双像三维建模小软件开发实例（六）

一、基本原理

根据相对定向或绝对定向确定了两张相片的相对位置关系之后，可以把原始影像纠正成核线影像，即两张相片的光轴平行，且与基线（相机头影像中心连线）垂直，同时核线影像的行（或列）与基线也保持平行，这时两张相片上的同名点将是行对准或列对准的，因此，寻找同名的过程将被限制到一维搜索。详细内容可参考计算机视觉的“对极几何”理论或摄影测量的“核线采样”理论。

原始影像其实是物方空间点到原始像空间的投影，而核线影像其实是物方空间点到基线坐标系的投影，二者的几何关系如下图所示，图中IMG1和IMG2是原始影像，EP1和EP2为相应的核线影像。

原始影像空间转换到物空间的关系为

物空间转换到核线影像空间的关系为

而原始影像到核线影像的转换关系则为

反过来，核线影像到原始影像的关系则为

二、算法实现

void EpipolarGenerator(double Bx, double By, double Bz, double *Rotationl,double *Rotationr, double psize, double f,   char *FilepathImgL, char *FilepathImgR, char *FilePathEpiImgL, char *FilePathEpiImgR){//读取原始影像CFxImage cImg1,cImg2;int bandMap[] = {1,2,3};cImg1.Open(FilepathImgL);int nWidth_l = cImg1.GetWidth();int nHeight_l = cImg1.GetHeight();cImg2.Open(FilepathImgR);int nWidth_r = cImg2.GetWidth();int nHeight_r = cImg2.GetHeight();if (nWidth_l != nWidth_r || nHeight_l != nHeight_r){return;}int nWidth, nHeight;nWidth = nWidth_l; nHeight = nHeight_l;unsigned char *ImgL = new unsigned char[nHeight*nWidth];unsigned char *ImgR = new unsigned char[nHeight*nWidth];memset(ImgL,0,sizeof(unsigned char)*nHeight*nWidth);memset(ImgR,0,sizeof(unsigned char)*nHeight*nWidth);cImg1.ReadInterleavingImage(0,0,nHeight,nWidth,1,bandMap,ImgL);cImg2.ReadInterleavingImage(0,0,nHeight,nWidth,1,bandMap,ImgR);cImg1.Close();cImg2.Close();//将当前坐标系转换至基线坐标系double B  = sqrt(Bx*Bx + By*By + Bz*Bz);double m_T = atan2(By,Bx);double m_V = asin(Bz/B);double Rtemp[9],Rl[9],Rr[9];GroundToBaseAuxMatrix(m_T,m_V,Rtemp); MultMatrix(Rtemp,Rotationl,Rl,3,3,3); MultMatrix(Rtemp,Rotationr,Rr,3,3,3);int EpipolarHeight, EpipolarWidth;GetEpiImgSize(nWidth,nHeight,f,psize,Rl,Rr,&EpipolarWidth,&EpipolarHeight);    unsigned char* EpiImgL = new unsigned char[EpipolarWidth*EpipolarHeight];unsigned char* EpiImgR = new unsigned char[EpipolarWidth*EpipolarHeight];//开始制作核线影像double u1,v1,x1,y1;double u2,v2,x2,y2;double c1,r1,c2,r2;int i,j;    double xo = nWidth/2;double yo = nHeight/2;for(i=0;i<EpipolarHeight;i++){for(j=0;j<EpipolarWidth;j++){u1 = j*psize-EpipolarWidth*psize/2;v1 = -i*psize+EpipolarHeight*psize/2;x1 = -f*(Rl[0]*u1 + Rl[3]*v1 - Rl[6]*f)/(Rl[2]*u1 + Rl[5]*v1 - Rl[8]*f);y1 = -f*(Rl[1]*u1 + Rl[4]*v1 - Rl[7]*f)/(Rl[2]*u1 + Rl[5]*v1 - Rl[8]*f);c1 = x1/psize + xo;r1 = -y1/psize + yo;if((c1 < 0.0) || (r1 < 0.0) || (r1 > nHeight-1 ) || (c1 > nWidth-1)){EpiImgL[i*EpipolarWidth+j] = 0;}else// 双线性内插{EpiImgL[i*EpipolarWidth+j] = BilinearInterpolation(ImgL,nWidth,c1,r1);}u2 = j*psize-EpipolarWidth/2*psize;v2 = v1;x2 = -f*(Rr[0]*u2 + Rr[3]*v2 - Rr[6]*f)/(Rr[2]*u2 + Rr[5]*v2 - Rr[8]*f);y2 = -f*(Rr[1]*u2 + Rr[4]*v2 - Rr[7]*f)/(Rr[2]*u2 + Rr[5]*v2 - Rr[8]*f);c2 = x2/psize + xo;r2 = -y2/psize + yo;if((c2 < 0.0) || (r2 < 0.0) || (r2 > nHeight-1) || (c2 > nWidth-1)){EpiImgR[i*EpipolarWidth+j] = 0;}else// 双线性内插{EpiImgR[i*EpipolarWidth+j] = BilinearInterpolation(ImgR,nWidth,c2,r2);}}}    int bandmap[] = { 1, 2, 3, 4};CFxImage EpiImgl,EpiImgr;EpiImgl.Create(FilePathEpiImgL,1,EpipolarWidth,EpipolarHeight);EpiImgl.WriteInterleavingImage(0,0,EpipolarHeight,EpipolarWidth,1,bandmap,EpiImgL);EpiImgr.Create(FilePathEpiImgR,1,EpipolarWidth,EpipolarHeight);EpiImgr.WriteInterleavingImage(0,0,EpipolarHeight,EpipolarWidth,1,bandmap,EpiImgR);EpiImgl.Close();EpiImgr.Close();delete []EpiImgL;delete []EpiImgR;delete []ImgL;delete []ImgR;} 

//地面坐标系到基线辅助坐标系的转换矩阵void GroundToBaseAuxMatrix( double t, double v, double matrix[9]){double sint = sin(t), cost = cos(t);double sinv = sin(v), cosv = cos(v);matrix[0] = cost*cosv;      matrix[1] = sint * cosv;         matrix[2] = sinv;matrix[3] = -sint;          matrix[4] = cost;                matrix[5] = 0;matrix[6] = -cost*sinv;     matrix[7] = -sint*sinv;          matrix[8] = cosv;}void GetEpiImgSize(int nWidth,int nHeight,double f, double psize,double *Rl, double *Rr, int *EpiWidth, int *EpiHeight){double uplx,uply,uprx,upry,downlx,downly,downrx,downry;double Oriuplx,Oriuply,Oriuprx,Oriupry,Oridownlx,Oridownly,Oridownrx,Oridownry;int xo = nWidth/2;int yo = nHeight/2;Oriuplx = -xo*psize;Oriuply = (nHeight - yo)*psize;Oriuprx = (nWidth - xo)*psize;Oriupry = (nHeight - yo)*psize;Oridownlx = -xo*psize;Oridownly = -yo*psize;Oridownrx = (nWidth - xo)*psize;Oridownry = -yo*psize;uplx = -f*(Rl[0]*Oriuplx + Rl[1]*Oriuply - Rl[2]*f)/(Rl[6]*Oriuplx+Rl[7]*Oriuply-Rl[8]*f);uply = -f*(Rl[3]*Oriuplx + Rl[4]*Oriuply - Rl[5]*f)/(Rl[6]*Oriuplx+Rl[7]*Oriuply-Rl[8]*f);uprx = -f*(Rl[0]*Oriuprx + Rl[1]*Oriupry - Rl[2]*f)/(Rl[6]*Oriuprx+Rl[7]*Oriupry-Rl[8]*f);upry = -f*(Rl[3]*Oriuprx + Rl[4]*Oriupry - Rl[5]*f)/(Rl[6]*Oriuprx+Rl[7]*Oriupry-Rl[8]*f);downlx = -f*(Rl[0]*Oridownlx + Rl[1]*Oridownly - Rl[2]*f)/(Rl[6]*Oridownlx+Rl[7]*Oridownly-Rl[8]*f);downly = -f*(Rl[3]*Oridownlx + Rl[4]*Oridownly - Rl[5]*f)/(Rl[6]*Oridownlx+Rl[7]*Oridownly-Rl[8]*f);downrx = -f*(Rl[0]*Oridownrx + Rl[1]*Oridownry - Rl[2]*f)/(Rl[6]*Oridownrx+Rl[7]*Oridownry-Rl[8]*f);downry = -f*(Rl[3]*Oridownrx + Rl[4]*Oridownry - Rl[5]*f)/(Rl[6]*Oridownrx+Rl[7]*Oridownry-Rl[8]*f);// 根据x、y的最大值和最小值确定核线影像的范围double a[4] = {uplx,uprx,downlx,downrx};OnBubbleSort(a,4);// 从小到大排序double x_Min = a[0];// 负值double x_Max = a[3];// 正值double b[4] = {uply,upry,downly,downry};OnBubbleSort(b,4);double y_Max = b[3];// 正值double y_Min = b[0];// 负值int EpipolarWidth1 = x_Max>-x_Min?(2*x_Max/psize):(-2*x_Min/psize);int EpipolarHeight1 = y_Max>-y_Min?(2*y_Max/psize):(-2*y_Min/psize);uplx = -f*(Rr[0]*Oriuplx + Rr[1]*Oriuply - Rr[2]*f)/(Rr[6]*Oriuplx+Rr[7]*Oriuply-Rr[8]*f);uply = -f*(Rr[3]*Oriuplx + Rr[4]*Oriuply - Rr[5]*f)/(Rr[6]*Oriuplx+Rr[7]*Oriuply-Rr[8]*f);uprx = -f*(Rr[0]*Oriuprx + Rr[1]*Oriupry - Rr[2]*f)/(Rr[6]*Oriuprx+Rr[7]*Oriupry-Rr[8]*f);upry = -f*(Rr[3]*Oriuprx + Rr[4]*Oriupry - Rr[5]*f)/(Rr[6]*Oriuprx+Rr[7]*Oriupry-Rr[8]*f);downlx = -f*(Rr[0]*Oridownlx + Rr[1]*Oridownly - Rr[2]*f)/(Rr[6]*Oridownlx+Rr[7]*Oridownly-Rr[8]*f);downly = -f*(Rr[3]*Oridownlx + Rr[4]*Oridownly - Rr[5]*f)/(Rr[6]*Oridownlx+Rr[7]*Oridownly-Rr[8]*f);downrx = -f*(Rr[0]*Oridownrx + Rr[1]*Oridownry - Rr[2]*f)/(Rr[6]*Oridownrx+Rr[7]*Oridownry-Rr[8]*f);downry = -f*(Rr[3]*Oridownrx + Rr[4]*Oridownry - Rr[5]*f)/(Rr[6]*Oridownrx+Rr[7]*Oridownry-Rr[8]*f);// 根据x、y的最大值和最小值确定核线影像的范围double c[4] = {uplx,uprx,downlx,downrx};OnBubbleSort(c,4);// 从小到大排序x_Min = c[0];// 负值x_Max = c[3];// 正值double d[4] = {uply,upry,downly,downry};OnBubbleSort(d,4);y_Max = d[3];// 正值y_Min = d[0];// 负值    int EpipolarWidth2 = x_Max>-x_Min?(2*x_Max/psize):(-2*x_Min/psize);int EpipolarHeight2 = y_Max>-y_Min?(2*y_Max/psize):(-2*y_Min/psize);*EpiWidth = EpipolarWidth1>EpipolarWidth2?EpipolarWidth1:EpipolarWidth2;    *EpiHeight = EpipolarHeight1>EpipolarHeight2?EpipolarHeight1:EpipolarHeight2;}

三、实验结果

以两幅航空影像为实验数据，测试本程序。

如下图所示为原始航空影像，由图可见，同名像点存在上下和左右视差。

如下图所示，为进过核线纠正后的影像，由图可见，同名像点之间只存在左右视差，而不存在上下时差，即同名点的行坐标是一致的。