📄 register.cpp
字号:
}
//////计算fx，fy，fx2，fy2及fxy
conv2( fx , imWidth , imHeight , prefilt , 1 , 3 );
conv2( fx , imWidth , imHeight , derivfilt , 3 , 1 );
conv2( fy , imWidth , imHeight , prefilt , 3 , 1 );
conv2( fy , imWidth , imHeight , prefilt , 1 , 3 );
ArraydotMultiArray( fx2 , fx , fx , imWidth , imHeight);
conv2( fx2 , imWidth , imHeight , blur , 1 , 3 );
conv2( fx2 , imWidth , imHeight , blur , 3 , 1 );
ArraydotMultiArray( fy2 , fy , fy , imWidth , imHeight);
conv2( fy2 , imWidth , imHeight , blur , 1 , 3 );
conv2( fy2 , imWidth , imHeight , blur , 3 , 1 );
ArraydotMultiArray( fxy , fx , fy , imWidth , imHeight);
conv2( fxy , imWidth , imHeight , blur , 1 , 3 );
conv2( fxy , imWidth , imHeight , blur , 3 , 1 );
delete[] fx;
delete[] fy;
float * mean_array = new float [lineByte * imHeight];
//detOfHessian为每个像素点处的特征矩阵
float * detOfHessian = new float [lineByte * imHeight];
float * norm_array = new float [lineByte * imHeight];
//qualityOfcorner为每个像素的特征的品质因素，值越大越代表是角点
float * qualityOfcorner = new float [lineByte * imHeight];
float maxquality = 0;
//thresh为阈值系数，它乘以最大品质因素即为阈值
float thresh = 0.05; // 0.005;
//距离边界小于border的像素不考虑
int border = 20; // do not consider the boundary pixels
float t1, t2;
for ( i = 0; i < imWidth ; i++){
for (j = 0 ; j < imHeight ; j++){
mean_array[j * lineByte + i]
=( fx2[j * lineByte + i] + fy2[j * lineByte + i])/2;
detOfHessian[j * lineByte + i]
= fx2[j * lineByte + i] * fy2[j * lineByte + i]
- pow(fxy[j * lineByte + i] , 2);
t1=pow(mean_array[j*lineByte+i] , 2);
t2=detOfHessian[j * lineByte + i];
norm_array[j*lineByte+i] = sqrt(t1-t2);
//边界点不予考虑
if (i < border || i > imWidth - border || j < border
|| j > imHeight - border)
{
qualityOfcorner[j * lineByte + i] = 0;
}
else
{
t1=abs(mean_array[j*lineByte+i] - norm_array[j*lineByte+i]);
t2=abs(mean_array[j*lineByte+i] + norm_array[j*lineByte+i]);
if(t1<t2)
qualityOfcorner[j * lineByte + i]= t1;
else
qualityOfcorner[j * lineByte + i]= t2;
}
///求最大的特征值
if (qualityOfcorner[j * lineByte + i] > maxquality)
maxquality = qualityOfcorner[j * lineByte + i];
}
}
delete[] fx2;
delete[] fy2;
delete[] fxy;
delete[] mean_array;
delete[] detOfHessian;
delete[] norm_array;
//thq为阈值
float thq = thresh * maxquality ;
//supneibor为取最大的邻域范围
int supneibor = 20;
//求局部最大
for ( i = supneibor/2 ; i < imWidth-supneibor/2 ; i++)
{
for ( j = supneibor/2 ; j < imHeight-supneibor/2 ; j++)
{
if ( qualityOfcorner[j * lineByte + i] > thq &&
LocalMaximum(qualityOfcorner, imWidth, imHeight,i, j,supneibor))
{
//满足不在边界范围和局部最大两个条件就认为是角点
pCornerpos->Add(CPoint(i,j));
}
}
}
delete[] qualityOfcorner;
//进一步根据距离进行抑制，两角点相距不能太近
SpaceSuppress(pCornerpos);
}
/***********************************************************************
* 函数名称：
* DrawCross()
*
*函数参数：
* unsigned char *imgBufIn, 待提取角点的图像数组
* int width, 图像的宽
* int height, 图像的高
* ptArray * pCornerpos ,存放特征点的矩阵
*返回值：
* 无
*
*说明：在图像的特征点的位置画十字，十字的宽为5，高为5
***********************************************************************/
void Register::DrawCross(unsigned char * imgBufIn , int imWidth ,
int imHeight , ptArray * pCornerpos )
{
int NumOfCorner = pCornerpos->GetSize();
CPoint pt;
for ( int m = 0 ; m< NumOfCorner ; m++)
{
//获得特征点
pt=pCornerpos->GetAt(m);
int x,y;
int i,j;
x = pt.x;
y = pt.y;
int lineByte=(imWidth+3)/4*4;
//画十字，即像素值设为255
for ( j = -2 ; j <= 2 ; j++ )
{
imgBufIn[(y+j) * lineByte + x] = 255;
}
for ( i = -2 ; i <= 2 ; i++ )
{
imgBufIn[y * lineByte + x+i] = 255;
}
}
}
/***********************************************************************
* 函数名称：
* HarisConner()
*
*函数参数：
* 无
*返回值：
* 无
*
*说明：调用HarisConnerDetect()函数检测图像1的角点，并DrawCross画出十字，
* 结果输出到m_pImgDataOut中，对外的接口
***********************************************************************/
void Register::HarisConner()
{
if (m_pImgData1==NULL)
{
return;
}
HarisConnerDetect(m_pImgData1,m_imgWidth1,m_imgHeight1,&m_cornerpos1);
if(m_pImgDataOut!=NULL){
delete []m_pImgDataOut;
m_pImgDataOut=NULL;
}
if(m_lpColorTableOut!=NULL){
delete []m_lpColorTableOut;
m_lpColorTableOut=NULL;
}
//输出图像的位数与输入图像相同
m_nBitCountOut=m_nBitCount1;
m_nColorTableLengthOut=ComputeColorTabalLength(m_nBitCount1);
//输出图像的宽和高，与输入图像相同
m_imgWidthOut=m_imgWidth1;
m_imgHeightOut=m_imgHeight1;
int lineByte=(m_imgWidthOut*m_nBitCountOut/8+3)/4*4;
int imgBufSize=m_imgHeightOut*lineByte;
m_pImgDataOut=new BYTE[imgBufSize];
memcpy(m_pImgDataOut, m_pImgData1, imgBufSize);
//如果有颜色表，分配颜色表所需要的空间，并拷贝颜色表
if(m_nColorTableLengthOut!=0){
m_lpColorTableOut= new RGBQUAD[m_nColorTableLengthOut];
memcpy(m_lpColorTableOut,m_lpColorTable1,
sizeof(RGBQUAD)*m_nColorTableLengthOut);
}
//在输出图像上画十字
DrawCross(m_pImgDataOut,m_imgWidthOut,m_imgHeightOut,&m_cornerpos1);
}
/***********************************************************************
* 函数名称：
* RegisterBasedConner()
*
*函数参数：
* unsigned char *imgBuf1, 图像1
* unsigned char *imgBuf2, 图像2
* int width, 图像的宽
* int height, 图像的高
* ptArray *pCornerpos1, 图像1单独提取的角点
* ptArray *pCornerpos2 , 图像2单独提取的角点
* ptArray *pCornerMatch1 , 配准后图像1中与图像2中匹配上的角点
* ptArray *pCornerMatch2 ,配准后图像2中与图像1中匹配上的角点
*返回值：
* 无
*
*说明：根据两幅图像各自提取的角点进行配准
***********************************************************************/
void Register::RegisterBasedConner(unsigned char *imgBuf1,unsigned char *imgBuf2,
int width, int height,ptArray *pCornerpos1,ptArray *pCornerpos2 ,
ptArray *pCornerMatch1 , ptArray *pCornerMatch2)
{
int cnt1 , cnt2;
cnt1 = pCornerpos1->GetSize();
cnt2 = pCornerpos2->GetSize();
int i , j , m ,n ;
int W = 15 ;
double aver1 , aver2;
double sigama1 , sigama2;
double temp;
int x1, y1, x2 , y2;
int lineByte = (width+3)/4*4;
int sigama = (int) (1.0/8*height);
double Cij;
//G为相关矩阵
double * G = new double[cnt1*cnt2];
double r;
//建立相关矩阵
for ( i = 0 ; i < cnt1 ; i++)
{
for ( j = 0 ; j < cnt2 ; j++)
{
//获取角点
x1 = pCornerpos1->GetAt(i).x;
y1 = pCornerpos1->GetAt(i).y;
x2 = pCornerpos2->GetAt(j).x;
y2 = pCornerpos2->GetAt(j).y;
//aver1、aver2为角点邻域像素的均值
aver1 = 0;
aver2 = 0;
//sigama1、sigama2为角点邻域像素的方差
sigama1 = 0;
sigama2 = 0;
temp=0;
//求取均值
for ( m = - (int) (W/2) ; m <= (int) (W/2) ; m++ )
{
for ( n = -(int) (W/2) ; n <= (int) (W/2) ; n++ )
{
aver1 += *(imgBuf1+lineByte*(y1+m)+x1+n);
aver2 += *(imgBuf2+lineByte*(y2+m)+x2+n);
}
}
aver1 /= (W*W);
aver2 /= (W*W);
//求取方差
for ( m = - (int) (W/2) ; m <= (int) (W/2) ; m++ )
{
for ( n = -(int) (W/2) ; n <= (int) (W/2) ; n++ )
{
sigama1 += pow(*(imgBuf1+lineByte*(y1+m)+x1+n)-aver1,2);
sigama2 += pow(*(imgBuf2+lineByte*(y2+m)+x2+n)-aver2,2);
temp += (*(imgBuf1+lineByte*(y1+m)+x1+n)-aver1)
*(*(imgBuf2+lineByte*(y2+m)+x2+n)-aver2);
}
}
sigama1=sqrt(double(sigama1));
sigama2=sqrt(double(sigama2));
Cij=temp/(W*W*sigama1*sigama2);
r=pow(x1-x2,2)+pow(y1-y2,2);
G[i*cnt2+j]=(Cij+1)/2*exp(-r/(2*sigama*sigama));
}
}
double * U = new double [cnt1*cnt1];
double * D = new double [cnt1*cnt2];
double * V = new double [cnt2*cnt2];
double eps = 1e-5;
int ka = max(cnt1,cnt2)+1;
//进行奇异值分解
dluav( G , cnt1 , cnt2 , U , V , eps , ka);
int length = min(cnt1,cnt2);
//将对角线上的值设为1
for ( i = 0 ; i < length ; i++ )
{
G[i*cnt2+i] = 1;
}
//P新的相关矩阵
double * P = new double [cnt1*cnt2];
damul( U , G , cnt1 , cnt1, cnt2 , D );
damul( D , V , cnt1 , cnt2 ,cnt2 , P );
ptArray pt1,pt2;
double maxrow = -1999;
double maxcol = -1999;
int imax,jmax;
for ( i = 0 ; i < cnt1 ; i++)
{
maxrow = -1999;
maxcol = -1999;
//求取该行最大的点所在的列
for ( j = 0 ; j < cnt2 ; j++)
{
if (P[i*cnt2+j] > maxrow)
{
maxrow = P[i*cnt2+j];
jmax = j ;
}
}
//求取jmax列最大的点所在的行
for ( m = 0 ; m < cnt1 ; m++)
{
if (P[m*cnt2+jmax] > maxcol)
{
maxcol = P[m*cnt2+jmax];
imax = m;
}
}
//判断是否既为行最大和列最大
if ( imax == i )
{
pCornerMatch1->Add(pCornerpos1->GetAt(i));
pCornerMatch2->Add(pCornerpos2->GetAt(jmax));
}
}
delete []U;
delete []V;
delete []G;
delete []P;
}
/***********************************************************************
* 函数名称：
* ConnerRegistering()
*
*函数参数：
* 1为成功,0为失败
*
*返回值：
* 无
*
*说明：调用matchBasedConner()函数对图像1和图像2进行角点匹配 ，并将配准后的
* 两幅图合并到一幅图中，并画出角点对的对应关系，对外的接口
***********************************************************************/
BOOL Register::ConnerRegistering()
{
if (m_pImgData1==NULL || m_pImgData2==NULL)
{
return false;
}
if (m_imgHeight1!=m_imgHeight2 || m_imgWidth1!=m_imgWidth2)
{
return false;
}
//图像1提取角点
HarisConnerDetect(m_pImgData1,m_imgWidth1,m_imgHeight1,&m_cornerpos1);
//图2提取角点
HarisConnerDetect(m_pImgData2,m_imgWidth2,m_imgHeight2,&m_cornerpos2);
//角点匹配
RegisterBasedConner(m_pImgData1,m_pImgData2,m_imgWidth2,m_imgHeight2,
&m_cornerpos1, &m_cornerpos2,&m_Matchpt1,&m_Matchpt2);
if(m_pImgDataOut!=NULL){
delete []m_pImgDataOut;
m_pImgDataOut=NULL;
}
if(m_lpColorTableOut!=NULL){
delete []m_lpColorTableOut;
m_lpColorTableOut=NULL;
}
//输出图像的位数与输入图像相同
m_nBitCountOut=m_nBitCount1;
m_nColorTableLengthOut=ComputeColorTabalLength(m_nBitCount1);
int lineByte=(m_imgWidth1*m_nBitCount1/8+3)/4*4;
int imgBufSize=m_imgHeight1*lineByte;
//输出图像的宽为原图的两倍，高与原图的相同
m_imgWidthOut=2*lineByte;
m_imgHeightOut=m_imgHeight1;
m_pImgDataOut=new BYTE[imgBufSize*2];
unsigned char *tempbuf1 = new unsigned char[lineByte*m_imgHeight1];
unsigned char *tempbuf2 = new unsigned char[lineByte*m_imgHeight2];
memcpy(tempbuf1,m_pImgData1,lineByte*m_imgHeight1);
memcpy(tempbuf2,m_pImgData2,lineByte*m_imgHeight1);
//画十字
DrawCross(tempbuf1,m_imgWidth1,m_imgHeight1,&m_Matchpt1);
DrawCross(tempbuf2,m_imgWidth2,m_imgHeight2,&m_Matchpt2);
//将匹配的两幅图合并为一幅图，第一幅图在左边，第二幅图在右边
for (int i = 0 ; i<lineByte;i++)
{
for (int j = 0 ; j<m_imgHeightOut;j++)
{
m_pImgDataOut[j*2*lineByte+i]=tempbuf1[j*lineByte+i];
m_pImgDataOut[j*2*lineByte+i+lineByte]=tempbuf2[j*lineByte+i];
}
}
CPoint pt1;
pt1=m_Matchpt1.GetAt(0);
//匹配的点对之间画线
DrawLine(m_pImgDataOut,m_imgWidthOut,m_imgHeightOut,&m_Matchpt1,&m_Matchpt2);
//如果有颜色表，分配颜色表所需要的空间，并拷贝颜色表
if(m_nColorTableLengthOut!=0){
m_lpColorTableOut= new RGBQUAD[m_nColorTableLengthOut];
memcpy(m_lpColorTableOut,m_lpColorTable1, sizeof(RGBQUAD)*m_nColorTableLengthOut);
}
delete[] tempbuf1;
delete[] tempbuf2;
return true;
}
/***********************************************************************
* 函数名称：
* DrawLine()
*
*函数参数：
* unsigned char *imgBufOut, 图像数组
* int widthOut, 图像的宽
* int heightOut, 图像的高
* ptArray *pcornerpos1, 存放合并前的第一幅图像的匹配上的角点, 即在原图中的坐标
* ptArray *pcornerpos2, 存放合并前的第二幅图像的匹配上的角点，即在原图中的坐标
*返回值：
* 无
*
*说明：匹配的点对之间画线，结果在imgBufOut中输出
***********************************************************************/
void Register::DrawLine(unsigned char *imgBufOut, int widthOut, int heightOut,
ptArray *pcornerpos1, ptArray *pcornerpos2)
{
double k;
int num = pcornerpos1->GetSize();
CPoint pt1;
CPoint pt2;
for ( int i=0; i<num ; i++)
{
pt1=pcornerpos1->GetAt(i);
pt2=pcornerpos2->GetAt(i);
pt2.x += widthOut/2;
//k为直线斜率
k = 1.0*(pt2.y-pt1.y)/(1.0*(pt2.x-pt1.x));
//画线，即将线上的像素值设为255
for (int x = pt1.x ; x<=pt2.x ; x++)
{
int y = (int)(pt1.y+k*(x-pt1.x));
imgBufOut[y*widthOut+x] = 255;
}
}
}
//实数矩阵的奇异值分解
//利用Householder变换及变形QR算法
//a-长度为m*n的数组，返回时其对角线依次给出奇异值，其余元素为0
//m-矩阵的行数
//n-矩阵的列数
//u-长度为m*m的数组，返回时存放奇异值分解的左奇异量U
//v-长度为n*n的数组，返回时存放奇异值分解的左奇异量V
//eps-精度要求
//ka-整型变量，其值为max(n,m)+1
//调用函数：dluav(),ppp(),sss()
int Register::dluav(double a[],int m,int n,double u[],double v[],double eps,int ka)
{
int i,j,k,l,it,ll,kk,ix,iy,mm,nn,iz,m1,ks;
double d,dd,t,sm,sm1,em1,sk,ek,b,c,shh,fg[2],cs[2];
double *s,*e,*w;
s=(double *)malloc(ka*sizeof(double));
e=(double *)malloc(ka*sizeof(double));
w=(double *)malloc(ka*sizeof(double));
it=60;
k=n;
if (m-1<n)
{
k=m-1;
}
l=m;
if (n-2<m)
{
l=n-2;
}
if (l<0)
{
l=0;
}
ll=k;
if (l>k)
{
ll=l;
}
if (ll>=1)
{
for (kk=1; kk<=ll; kk++)
{
if (kk<=k)
{
d=0.0;
for (i=kk; i<=m; i++)
{
ix=(i-1)*n+kk-1;
d=d+a[ix]*a[ix];
}
s[kk-1]=sqrt(d);
if (s[kk-1]!=0.0)
{
ix=(kk-1)*n+kk-1;
if (a[ix]!=0.0)
{
s[kk-1]=fabs(s[kk-1]);
if (a[ix]<0.0)
{
s[kk-1]=-s[kk-1];
}
}
for (i=kk; i<=m; i++)
{
iy=(i-1)*n+kk-1;
a[iy]=a[iy]/s[kk-1];
}
a[ix]=1.0+a[ix];
}
s[kk-1]=-s[kk-1];
}
if (n>=kk+1)
{
for (j=kk+1; j<=n; j++)
{
if ((kk<=k)&&(s[kk-1]!=0.0))
{
d=0.0;
for (i=kk; i<=m; i++)
{
ix=(i-1)*n+kk-1;
iy=(i-1)*n+j-1;
d=d+a[ix]*a[iy];
}
d=-d/a[(kk-1)*n+kk-1];
for (i=kk; i<=m; i++)
{
ix=(i-1)*n+j-1;
iy=(i-1)*n+kk-1;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -