大家一起讨论一下图象拼接问题！

摘要：相位相关度是一种使用傅立叶变换的全局图像匹配方法，如果直接按照傅立叶变换公式进行计算，每对一行长度为N的像素进行变换，需要N2次乘法和N*(N-1)次加法，对于N较大时，所需的计算量相当大。本文使用一种快速傅立叶变换及编码的方法，通过对N进行因子拆分来加快图像的处理速度，从而达到实用的目的。在得到图像偏移量之后，我们用寻找最佳缝合线的方法来合成拼图，以消除拼接过程中可能出现的鬼影现象。
关键词图像拼接，相位相关度，快速傅立叶变换，最佳缝合线1 引言
全景拼图（Panorama Mosaic）是一种应用较广的基于图像的绘制技术（Image Based Rendering）。使用这种技术，可以把数码相机（甚至普通相机）围绕相机光心旋转拍摄的一系列照片合成为新的全景图。全景图的生成过程分为两个主要的步骤：图像配准和最终的图像缝合。
图像配准技术主要有模板匹配、相位相关度、以及基于运动（Motion-Based）的方法[6]等。模板匹配算法的基本思路是：首先由图像像素灰度值信息提取特征区域，然后以一幅图像中的特征区域作为模板，在另一幅图像中搜索最佳匹配。这种算法的关键在于中心灰度值的确定和灰度波动值的选取，这直接关系到特征区域的拼接能否实现，然而不同图像间灰度值的变化范围较大，使得难以精确选取所需要的值。基于运动的方法采用传统的8个未知参数矩阵来描述两幅图像之间的变换关系，并用最小化算法估计该矩阵，然后根据矩阵之间的参数约束关系计算相机焦距，再代入焦距值，进一步求得摄像机坐标间的旋转矩阵。这种方法在迭代求精之前必须均衡图像间的光照差，并且需要比较精确的迭代初值。
相位相关度方法，对具有重叠区域的相邻两幅图像的灰度值进行傅立叶变换以求得图像间的偏移量。与基于运动的方法和模板匹配相比，相位相关度方法对两幅图像间的曝光差异并不敏感，并且没有阈值（模板匹配）和最小二乘法迭代初值（基于运动的方法）难以确定的问题，在匹配的精度方面要优于模板匹配的方法。图像系列先变换到柱面坐标，再用相位相关度求解，就可以得到精确的偏移量。相位相关度方法采用全局匹配方式，对于图像间重叠区域的大小并没有太多限制，其他方法一般要求重叠区域在20%-40%之间，因为计算量会随着搜索区域的增大而增加，相位相关度则不存在这方面的问题。
然而，相位相关度方法的缺点是需要作大量的运算，每拼接两幅图像需要做三次二维傅立叶变换，二维傅立叶变换公式如下：
                    （1）
这里N和M为图像X和Y方向上的像素个数，f(x,y)为图像的灰度值函数，u,v为傅立叶频域坐标。上述二维傅立叶变换可以拆作两个一维傅立叶变换，即先计算上式右边对x的求和，把求和结果作为新的f(x,y)，再对y进行一维傅立叶变换。对于长度为N的一维变
换，共需要N2次乘法和N*(N-1)次加法。当N较大时(一般图像的两个一维值在几百至一千多之间)，所需的运算量相当大，这直接导致拼接速度过慢。本文使用一种快速傅立叶变换来减少运算量，从而达到实用的目的。
在求得图像间的偏移量后，接下来就是合成新的图像，由于场景中同一物体在相邻两幅图像重叠区域中的形状有所不同，或者场景中的移动物体在重叠区域内只存在于一幅图像中，这样在合成的新图像中就会有模糊的影像，称为鬼影(Ghost)。我们采用最佳缝合线的方法[5]，[8]来消除鬼影。
本文接下来将首先介绍图像的柱面映射，然后再阐述所使用的快速傅立叶变换，最后讨论图像缝合及曝光差异的消除。
2.1 柱面坐标映射
由于相机是围绕固定中心点旋转拍摄场景，两幅图像中的同一场景是在不同的视角下拍摄的，而所获得的图像是场景的平面投影，因而在图像的重叠区域内具有不同的位置和形状，这种差异在近景拍摄和室内场景拍摄时尤为明显。我们在使用相位相关度进行图像配准前，先把图像映射到柱面坐标，以减少重叠区域内场景的差异：

                                                (2)
这里x和y是原来的像素平面坐标，s和t是新的柱面坐标，f是相机的焦距。我们在处理图像时，如果不知道相机的焦距，就先对两幅相邻图像在平面坐标下使用相位相关度求取水平偏移量x0，取f=n(N-x0)/2л作为柱面变换的近似焦距，n为旋转一周拍摄到的图像总数目，N为单幅图像x方向上像素的个数。此外我们还做了一个焦距变换菜单，当知道焦距时，可以直接设定焦距。2.2 相位相关度
相位相关度的主要思想是：相邻两幅图像的像素值f1和f2，存在着：f2（x，y）=f1（x-x0，y-y0）, x0,y0是图像之间的平移量（拍摄时，使相机的光心与地面垂直，再把所得到的图像进行柱面映射，那么图像间的偏移就可以看作是平移运动）。设f（x，y）的傅立叶变换为F（u，v），根据傅立叶变换的性质有( 表示函数和其傅立叶变换的对应性)：
                    (3)
两幅图像在傅立叶变换域上的偏移定义为：
                     (4)

其中F2*(u,v)是f2(x,y)傅立叶变换的共轭函数。对上式右边作反向傅立叶变换，并寻找变换后具有实部最大值的像素点，该像素点坐标就是所求的x和y方向的偏移量[4]。2.2 快速傅立叶变换及编码
目前介绍快速傅立叶变换的图像处理书籍，都只涉及到一维长度N=2m的情况：采用蝶形算法或者逐次叠加法，可以把运算量减为N*log2N，然而图像的高度和宽度，通常并不能够满足为2m这种要求，如果对图像截取的太多，就会导致匹配不准。在这里我们采用Cooley-Tukey算法[3]来减少运算量，如果一维傅立叶变换的长度N不是素数，则N=n1*n2，0<n1,n2<N，且为整数。则一维傅立叶变换：
             j, k = 0,1,…, N-1               (5)
改写为：
    (6)
j=n2j1+j2,     j1，k1=0,…n1-1
     k=n1k2+k1       j2，k2=0，…n2-1。
我们注意到，是以M为周期的函数，存在着，0<i<M/2，即与为共轭函数，因此可以把对称的乘法简化为加法来运算。在公式(6)中，如果n1和n2不为素数并且大于10，还可以继续递归拆分，直到是素数或者小于10为止。
我们在程序中，事先把n为2至10的傅立叶变换因子进行运算和编码，对于大于10的素数因子，采用对半拆分的方法来运算。具体步骤为：先对N进行因子拆分，然后根据拆分的结果对像素的存储位置进行重排列，最后根据不同的因子循环调用m次相应的傅立叶变换函数（m为N的因子个数）。如果拆分的因子中有2和8，那么我们把2×8变为4×4，这是因为，对于因子4只需要对相应的像素做加减法即可。在实际处理图像时，我们把图像的宽度和高度截去几个象素，即N=(int)(N/10)*10，以便于对N进行因子拆分，又不影响全局变换。在求得偏移量后加上截去的像素，就是所求的图像位移。3 图像缝合
在取得图像间的位移后，我们用寻找最佳缝合线的方法来消除合成图像中的鬼影[5]，[8]。寻找最佳缝合线的方法：
                   (7)
这里M是重叠区域的高度，Ecolor是相邻两幅图像重叠区域内对应像素的颜色差值，在这里是对相应的像素灰度值进行计算。Egemetry是相应的像素结构差值，这里我们使用两个改进的sobel加权算子Sx和Sy来求取x和y方向上结构差：
  和   求解过程为：把重叠区域内第一行的第N个像素作为第N条缝合线的起始像素，再根据公式(7)对每一条缝合线的当前点以及与该点紧邻的下一行中的三个像素值求取颜色差值和结构差值，分别相加进行比较，从中选取最小的计算结果，并以求得最小值的像素的位置作为新的当前点，继续求解下一行的三个相邻像素。以此类推，直到最后一行，记录求取的像素位置作为一条缝合线。用同样的方法可以求得N条缝合线(N为重叠区域的宽度)，从中选取值最小的一条作为拼接缝合线。
图1说明了生成缝合线的开始两步，左图为缝合线的起始点，右图为对第二行进行搜索：经过第一行每个像素点的缝合线发出了三根线段（对于边缘像素是两根）指向下一行紧邻的三个像素点，实线段则表示该缝合线的最小值的扩展方向。

             左                            右
                 图1  最佳缝合线搜索
在找到拼接缝合线后，我们用左，右两边的图像分别填充缝合线左，右两边的区域。如果图像间存在着明显的曝光差异，那么在合成后的图像缝合线两边也存在着这种差异。解决方法有，对重叠区域进行加权计算[6]，以形成平滑过渡，这种方法简单、运算少但效果并不理想。Bert[7] 提出一种多分辨率方法，能较好的解决曝光差异问题，并且可以对前面的配准结果作进一步的校正，然而所需的计算量要比前一种方法大得多，而且由于是对原有像素在x方向进行扩展，因此会在合成图像的个别区域内出现失真的情况。
在这里我们求取相邻两幅图像在重叠区域内对应像素灰度值和的比值：
                                                (8)
S为重叠区域内的所有像素；f1(x,y)，f2(x,y)为左右相邻图像的灰度值函数。如果N>1，则对左图的所有像素的R,G,B值除以N，如果N<1，则对右图进行相应的处理。由于相邻两幅图像在经过配准之后，重叠区域内的像素是一一对应的，所以比值N在某种程度上能够反映两幅图像间的曝光差异。  4 实验结果
   对于12幅422×288的图像，我们在一台P4  1.7G，512M内存的机子上运行，图像柱面映射加上拼接时间共花费30秒，8幅1024×768的图像共需1分10秒，拼接效果：绝大部分图像的拼接，在自动拼接完成后，不需要再进行手工校正。图2显示了室内场景未作柱面坐标映射直接拼接的效果图，绿色边框显示源图的区域。图3显示了室内场景作柱面坐标映射后的效果图。图4显示了昆明世博园的一个场景拼接后的效果图。图5显示了最后合成的图像。

                      图2  未作柱面坐标变换的室内场景拼接结果

                      图3  作柱面坐标变换后的室内场景拼接结果
                        图4   昆明世博园一场景拼接后的图像

                        图5   合成的图像5 结论
本文提出了一种基于相位相关度的快速图像拼接方法，大大加快了傅立叶变换的图像处理速度，从而达到实用的目的。我们在拼接前对图像进行裁剪，以便于拆分素数因子，加快拼接速度，但这样就导致第i和i+1幅图像拼接完成后，对第i+1和i+2 幅图像进行处理时，需要重新计算一次第i+1幅图像的傅立叶变换。我们以后的工作将进一步优化傅立叶编码，以消除这种重复计算。
在求得图像间的偏移量后，我们寻找最佳缝合线来消除合成图像时可能出现的鬼影现象，实验结果显示了这种方法具有良好的合成效果。对于图像间曝光差异所带来的问题，也是我们今后需要进一步研究的课题，以求找到一种既快速又有良好效果的方法。参考文献
听说IPP是用傅立叶变换做的,不知是真是假.我以前写了一些,没有成功,以下是我的代码,看看给各位有没有帮助.如果大家调好了记得发上来呀.
/*
void CCh1_1App::OnPictureLink()
{
int i,j,ii;
int nDocCount=0;//记录文档的数量
CCh1_1Doc *pDoc[5];//文档指针
HDIB  hDIB[5];//图像句柄
LPSTR lpDIB[5];//指向DIB的指针
LPSTR lpDIBBits[5];//指向DIB数据区的指针
for(i=0;i<5;i++)
{
pDoc[i]=NULL;
hDIB[i]=NULL;
lpDIB[i]=NULL;
lpDIBBits[i]=NULL;
}
POSITION position=m_pDocTemplate->GetFirstDocPosition( );
for(i=0;position!=NULL;i++,nDocCount++)
{
pDoc[i]=(CCh1_1Doc *)m_pDocTemplate->GetNextDoc(position) ;
hDIB[i]=pDoc[i]->GetHDIB();
lpDIB[i]=(LPSTR) ::GlobalLock(hDIB[i]);
lpDIBBits[i] = ::FindDIBBits(lpDIB[i]);
}
if(nDocCount<=1)
{
AfxMessageBox("请打开两幅以上大小相等的图像");
return;
}
if(nDocCount>5)
{
AfxMessageBox("最多只能打开五幅大小相等的图像");
return;
}
LPSTR lpNew;
LPSTR lpData[5],lpDest;
char * upNewBMP;
// 图像每行的字节数
LONG lLineBytes;
int lWidth,lHeight;
lWidth=::DIBWidth(lpDIB[0]);
lHeight=::DIBHeight(lpDIB[0]);
// 计算图像每行的字节数
lLineBytes = WIDTHBYTES(lWidth * 8);//*PictureBits);
m_hNewDIB =  ::GlobalAlloc(GHND,lWidth*lHeight+*(LPDWORD)lpDIB[0] + ::PaletteSize(lpDIB[0]));
upNewBMP=( char * )::GlobalLock(m_hNewDIB);
memcpy( upNewBMP,lpDIB[0], (lpDIBBits[0]-lpDIB[0])*8 );
::GlobalUnlock(m_hNewDIB);
lpNew=::FindDIBBits(upNewBMP);
// 更改光标形状
BeginWaitCursor();
// 循环控制变量
int y;
int x;

// 临时变量
double dTmpOne;
double  dTmpTwo;

// 傅立叶变换竖直方向点数
int nTransHeight ; // 傅立叶变换水平方向点数
int nTransWidth  ;

// 计算进行傅立叶变换的点数（2的整数次幂）
dTmpOne = log(lWidth)/log(2);
dTmpTwo = ceil(dTmpOne)    ;
dTmpTwo = pow(2,dTmpTwo)    ;
nTransWidth = (int) dTmpTwo    ;

// 计算进行傅立叶变换的点数（2的整数次幂）
dTmpOne = log(lHeight)/log(2);
dTmpTwo = ceil(dTmpOne)    ;
dTmpTwo = pow(2,dTmpTwo)    ;
nTransHeight = (int) dTmpTwo   ; // 计算图象数据存储每行需要的字节数
// BMP文件的每行数据存储是DWORD对齐的
int nSaveWidth;
nSaveWidth = ( (lWidth << 3) + 31)/32 * 4 ;
// 图象象素值
unsigned char unchValue; complex<double> *TD = new complex<double>[nTransWidth * nTransHeight];
complex<double> *FD = new complex<double>[nTransWidth * nTransHeight];
complex<double> *TD1 = new complex<double>[nTransWidth * nTransHeight];
complex<double> *FD1 = new complex<double>[nTransWidth * nTransHeight];
complex<double> *TD3 = new complex<double>[nTransWidth * nTransHeight];
// 初始化
// 图象数据的宽和高不一定是2的整数次幂，所以pCTData
// 有一部分数据需要补0
for(y=0; y<nTransHeight; y++)
{
for(x=0; x<nTransWidth; x++)
{
TD[y*nTransWidth + x]=complex<double>(0,0);
TD1[y*nTransWidth + x]=complex<double>(0,0);
}
} // 把图象数据传给pCTData
for(y=0; y<lHeight; y++)
{
for(x=0; x<lWidth; x++)
{
unchValue = lpDIBBits[0][y*nSaveWidth +x];
TD[y*nTransWidth + x]=complex<double>(unchValue,0);
unchValue = lpDIBBits[0][y*nSaveWidth +x];
TD1[y*nTransWidth + x]=complex<double>(unchValue,0);
}
} // 傅立叶正变换
DIBFFT_2D(TD, lWidth, lHeight, FD) ;
// 傅立叶正变换
DIBFFT_2D(TD1, lWidth, lHeight, FD1) ;

// 中间变量
double dTemp;

// 进行付立叶变换的宽度和高度（2的整数次方）
LONG w;
LONG h;

int wp;
int hp;

// 赋初值
w = 1;
h = 1;
wp = 0;
hp = 0;

// 计算进行付立叶变换的宽度和高度（2的整数次方）
while(w * 2 <= lWidth)
{
w *= 2;
wp++;
}

while(h * 2 <= lHeight)
{
h *= 2;
hp++;
}

//以下计算偏移量
CString ss;
int nDpos;
double rcps, icps, mag;
for(i = 0; i < h; i++)
{
for(j = 0; j < w; j++)
{
nDpos=j+i*w;
rcps=FD[nDpos].real()*(FD1[nDpos].real()) + FD[nDpos].imag()*(FD1[nDpos].imag());
icps=FD[nDpos].imag()*(FD1[nDpos].real()) - (FD1[nDpos].imag())*(FD[nDpos].real());

mag=sqrt(rcps*rcps+icps*icps);

if(mag==0)
{
rcps=0; icps = 0;   // Just in case
}
else
{
rcps /= mag;
icps /= mag;
}

FD[nDpos]=complex<double>( rcps, icps);
}
}
IFFT_2D(FD,TD3,lWidth,lHeight);
int peakx = 0, peaky = 0;      // Position of the peak on the surface
float height = 0;
float peakh = 0;
int dpos;
CString s;
for(i = 0; i < h; i++)
{
for(j = 0; j < w; j++)
{
dpos = j+i*w;
height = TD3[dpos].real();
lpData[1]=lpDIBBits[1]+lLineBytes*( i) +j; *lpData[1]=height;
if(height > peakh)
{
peakx = j; peaky = i;
peakh = height;

}

}
}
//for (i = 0; i <lHeight; i ++)
{

// lpData[1]=lpDIBBits[1]+lLineBytes*(lHeight - 1 -  i) + peakx;
// *lpData[1]=0;

}

s.Format("%d ,%d ,%f",peakx,peaky,peakh);
AfxMessageBox(s);
for(i=0;i<nDocCount;i++)
::GlobalUnlock(hDIB[i]);
// ((CMainFrame *)m_pMainWnd)->CreateBMP(m_hNewDIB);
delete TD,FD,TD1,FD1,TD3;

}
*/

调试易

大家一起讨论一下图象拼接问题！

解决方案 »