基于相关系数影像匹配实习报告

基于相关系数影像匹配实习报告一、实习内容与目的通过上机编程实现基于相关系数算法的影像匹配，从而更加熟练的掌握书本上关于影像匹配及相关系数算法的理论知识，并熟悉上机编程的操作。

二、实习原理影像匹配实质上是在两幅或是多幅影像之间识别同名点，其中基于相关系数的算法是实现影像匹配的基本算法之一。

相关系数是标准化的协方差函数，协方差函数除以两信号的方差即得相关系数。

计算相关系数的公式如下：g(x,y)与g ’(x ’,y ’)的相关系数为：),(),(),(q p C C q p C q p g g gg ''=ρ其中：⎰⎰∈-=Dy x gg dxdy y x g E y x g C ),(2)]},([),({ ⎰⎰∈''++'-++'=Dy x g g dxdy q y p x g E q y p x g q p C ),(2)]},([),({),( 若(p0, q0)>(p, q)，(p ≠p0,q ≠q0)，则p0,q0为搜索区影像相对于目标区影像的位移参数。

由离散灰度数据对相关系数的估计为∑∑∑∑∑∑====++==++⋅'-'⋅-'-'-=m i njm i njc rcj r i j i m i njc j r i j i g g g g g g g g r c 11112,,2,11,,)()())((),(ρ∑∑==++'⋅='m i njcj r i rc g nm g 11,,1，∑∑==⋅=m i njj i g nm g 11,1考虑到计算工作量，相关系数的实用公式为∑∑∑∑∑∑∑∑∑∑∑∑∑∑======++++======++==++'⋅-'⋅-'⋅-⋅=m i nj m i n j m i nj c j r ic j r i m i n j j i ji m i njm i n j c j r im i n j j i c j r i j i g n m g g n m gg g n m g g r c 1111211,,2211,,21111,11,,,])(1][)(1[)()(1)(),(ρ三、实习步骤及相关代码在第一次点特征点实习的基础上我们组进行了第二次基于相关系数的影像匹配实习。

所以程序实现的部分操作直接利用前次实习的中的程序。

（1）读入左右影像 FILE* pSrcFileLeft = NULL;pSrcFileLeft = fopen(srcFileLeft,"rb");if(pSrcFileLeft == NULL){printf("原始左片影像打开失败\n");return FAILURE;}FILE* pSrcFileRight = NULL;pSrcFileRight = fopen(srcFileRight,"rb");if(pSrcFileRight == NULL){printf("原始右片影像打开失败\n");return FAILURE;}（2）确定目标窗口大小通过老师课堂上的讲授和我们在编程过程中的不断尝试，我们默认使用的窗口的尺寸为11*11。

（3）确定目标点的位置由前次的实习中得到左相片中的特征点作为这次实习的目标点。

（4）预测右影像搜索范围本程序中也在左相片中计算出了特征点，然后将左右两相片中的特征点进行计算相关系数。

（5）逐窗口计算相关系数并保存相应的储存单元将左相片中的目标点与右相片中的特征点进行相关系数的计算。

（6）比较相关系数，取最大相关系数作为匹配点对于分别比较上一步骤中计算出的相关系数，选取出最大的一个相关系数作为匹配点。

程序操作如下：在程序主界面中点击Extend中选取ImageMatch选项，然后进入参数设置界面。

设置读取和保存影像路径和窗口大小等参数后，点击Progress进行运算。

四、实习结果及分析作为数据源的左右影像如些图所示：左相片右相片经过计算后得到的目标点左右图像分别如下所示：左相片右相片显示匹配点号：通过老师得到的图片可以看出在左相片中得到的五个目标点中有四个在右相片中找到了同名点，由此可以看出本程序基本满足实习要求，并得到理想的结果。

不足的地方就是选取的特征点的数量不够多，这样就无法判断本程序的“误判”的情况如何。

五、实习体会在学习到相关系数的时候，一看到计算相关系数的公式的时候我就有种被折服的感觉，那个公式实在是太复杂，它使我感觉我怎么也不能记住这个复杂的公式。

不过经过这次实习后，我觉得这个公式又不是当初认为的那样无法记忆。

我想这次实习我最大的收获就是对于之前在学习过程中遇到不懂知识以及难以记忆的知识都得到了加强，方便了我对他们的理解与记忆。

其次，在编程过程中更加锻炼自己动手编程的能力，使得我编程能力较之前有着一定的提高。

每次编程实习都是一个积累经验的过程，然后通过一次次的实习后，我最终的能力就能够得到很好的体现。

附件：void Moravec(BYTE* pSrcBits, vector<FEATUREPOINT> *FeaturePoint,int m_threshold, int m_window1, int m_window2, int width, int height, int* Count1){。

此部分程序与点特征程序中代码一样，因此不再重复粘贴}vector<FEATUREPOINT> RemoveReplicative(vector<FEATUREPOINT> v){vector<FEATUREPOINT> ret;vector<FEATUREPOINT>::iterator iter = v.begin();ret.clear();ret.push_back(*iter);vector<FEATUREPOINT>::iterator iter2;BOOL b(FALSE);for (iter = v.begin(); iter != v.end(); ++iter){b = false;for (iter2 = ret.begin(); iter2 != ret.end(); ++iter2){if (iter->x == iter2->x && iter->y == iter2->y){b = TRUE; //存在break;}}if (b == FALSE){ret.push_back(*iter);}}return ret;}void SaveBand(int width, int height, int byteCount, int biBitCount, LPBYTE pBits, CString SavePath){BITMAPFILEHEADER bmfh;BITMAPINFOHEADER bmih;bmfh.bfType = 0x4d42; // 0x42 = "B" 0x4d = "M"bmfh.bfReserved1 = 0;bmfh.bfReserved2 = 0;if (biBitCount == 8){bmfh.bfOffBits = sizeof(BITMAPFILEHEADER) + sizeof(BITMAPINFOHEADER) + 256*4;bmfh.bfSize = sizeof(BITMAPFILEHEADER) + sizeof(BITMAPINFOHEADER) + 256*4 + byteCount;}else if (biBitCount == 24){bmfh.bfOffBits = sizeof(BITMAPFILEHEADER) + sizeof(BITMAPINFOHEADER);bmfh.bfSize = sizeof(BITMAPFILEHEADER) + sizeof(BITMAPINFOHEADER) + byteCount;}else{MessageBox(NULL, "8bit or 24bit accepted", "错误", MB_ICONERROR);return;}bmih.biBitCount = biBitCount;bmih.biWidth = width;bmih.biHeight = height;bmih.biSize = 40;bmih.biPlanes = 1;bmih.biCompression = BI_RGB;bmih.biSizeImage = bmfh.bfSize;bmih.biClrImportant = 0;bmih.biClrUsed = 0;bmih.biXPelsPerMeter = 0;bmih.biYPelsPerMeter = 0;CFile f;//保存位图if (f.Open(SavePath, CFile::modeCreate | CFile::modeWrite)){f.Write(&bmfh, sizeof(bmfh));f.Write(&bmih, sizeof(bmih));if (biBitCount == 8){RGBQUAD rgb[256];int i(0);for (; i < 256; ++i){rgb[i].rgbBlue = i;rgb[i].rgbGreen = i;rgb[i].rgbRed = i;rgb[i].rgbReserved = 0;}//修改颜色表255对应的颜色rgb[255].rgbBlue = 0;rgb[255].rgbGreen = 0;rgb[255].rgbRed = 255;f.Write(rgb, sizeof(RGBQUAD)*256);}f.Write(pBits, byteCount);f.Close();}else{MessageBox(NULL, "保存图像失败", "错误", MB_ICONERROR);return;}}void DrawCross(BYTE* pBits, vector<FEATUREPOINT> FeaturePoint, int _width, int byteCount){vector<FEATUREPOINT>::iterator iter;int i(0), j(0);for (i = 0; i < byteCount; ++i){if (pBits[i] == 255)pBits[i] = 254;}for (iter = FeaturePoint.begin(); iter != FeaturePoint.end(); ++iter){i = iter->y;j = iter->x;pBits[i*_width+j] = 255;pBits[i*_width+j+1] = 255;pBits[i*_width+j-1] = 255;pBits[(i+1)*_width+j] = 255;pBits[(i-1)*_width+j] = 255;pBits[i*_width+j+2] = 255;pBits[i*_width+j-2] = 255;pBits[(i+2)*_width+j] = 255;pBits[(i-2)*_width+j] = 255;}}void DrawCross(BYTE* pBits1, BYTE* pBits2, vector<MATCHEDPOINTS> MatchedPoints, int _width1, int _width2, int byteCount1, int byteCount2){vector<MATCHEDPOINTS>::iterator iter;int i(0), j(0);for (i = 0; i < byteCount1; ++i){if (pBits1[i] == 255)pBits1[i] = 254;}for (j = 0; j < byteCount2; ++j){if (pBits2[j] == 255)pBits2[j] = 254;}for (iter = MatchedPoints.begin(); iter != MatchedPoints.end(); ++iter){i = iter->y1;j = iter->x1;pBits1[i*_width1+j] = 255;pBits1[i*_width1+j+1] = 255;pBits1[i*_width1+j-1] = 255;pBits1[(i+1)*_width1+j] = 255;pBits1[(i-1)*_width1+j] = 255;pBits1[i*_width1+j+2] = 255;pBits1[i*_width1+j-2] = 255;pBits1[(i+2)*_width1+j] = 255;pBits1[(i-2)*_width1+j] = 255;i = iter->y2;j = iter->x2;pBits2[i*_width2+j] = 255; //pBits2[i*_width2+j+1] = 255;pBits2[i*_width2+j-1] = 255;pBits2[(i+1)*_width2+j] = 255;pBits2[(i-1)*_width2+j] = 255;pBits2[i*_width2+j+2] = 255;pBits2[i*_width2+j-2] = 255;pBits2[(i+2)*_width2+j] = 255;pBits2[(i-2)*_width2+j] = 255;}}void SaveList(vector<FEATUREPOINT> FeaturePoint, CString txtFile)FILE* pTxtFile= fopen(txtFile, "wt");if (!pTxtFile)return;fprintf(pTxtFile, "#\t特征点数: %d\n\n", FeaturePoint.size());fprintf(pTxtFile, "坐标X\t坐标Y\t兴趣值\n");vector<FEATUREPOINT>::iterator iter;for (iter = FeaturePoint.begin(); iter != FeaturePoint.end(); ++iter)fprintf(pTxtFile, "%d\t%d\t%d\n", iter->x, iter->y, (int)iter->IV);fclose(pTxtFile);}void SaveList(vector<MATCHEDPOINTS> MatchedPoints, CString txtFile){FILE* pTxtFile= fopen(txtFile, "wt");if (!pTxtFile)return;fprintf(pTxtFile, "#同名点对数: %d\n\n",MatchedPoints.size());vector<MATCHEDPOINTS>::iterator iter;for (iter = MatchedPoints.begin(); iter != MatchedPoints.end(); ++iter)fprintf(pTxtFile, "%d\t%d\t%d\t%d\t%f\n",iter->x1, iter->y1, iter->x2, iter->y2, iter->coefficient);fclose(pTxtFile);}int _tmain(int argc, TCHAR* argv[], TCHAR* envp[]){if (argc != 15){::MessageBox(NULL,"启动参数不正确,请从GUI.exe或命令行参数启动","",MB_ICONASTERISK);return FAILURE;}system("Mode con: COLS=50 LINES=150");system("title Image_Match");int nRetCode = 0;if (!AfxWinInit(::GetModuleHandle(NULL), NULL, ::GetCommandLine(), 0)){cerr << _T("Fatal Error: MFC initialization failed") << endl;nRetCode = 1;return nRetCode;}CString srcFileLeft, srcFileRight, dstFileLeft, dstFileRight, txtFile;CString threshold_Left,window1_Left,window2_Left;CString threshold_Right,window1_Right,window2_Right;CString CoRelated,window3,flag1; //启动参数srcFileLeft = argv[1];srcFileRight = argv[2];dstFileLeft = argv[3];dstFileRight = argv[4];txtFile = argv[5];threshold_Left = argv[6]; //阈值window1_Left = argv[7];window2_Left = argv[8];threshold_Right = argv[9];window1_Right = argv[10];window2_Right = argv[11];CoRelated = argv[12];window3 = argv[13];flag1 = argv[14]; //判断是否在创建结果文件后显示int m_threshold_L = atoi(threshold_Left);int m_window1_L = atoi(window1_Left);int m_window2_L = atoi(window2_Left);int m_threshold_R = atoi(threshold_Right);int m_window1_R = atoi(window1_Right);int m_window2_R = atoi(window2_Right);double m_CoRelated = (double)atoi(CoRelated) / 100;int m_window3 = atoi(window3);int m_flag1 = atoi(flag1);if (GetFileAttributes(srcFileLeft) == -1){CString err;err.Format("文件:%s不存在", srcFileLeft);::MessageBox(NULL,err, "Error", MB_ICONERROR);return FAILURE;}if (GetFileAttributes(srcFileRight) == -1){CString err;err.Format("文件:%s不存在", srcFileRight);::MessageBox(NULL,err, "Error", MB_ICONERROR);return FAILURE;}//打开2幅影像FILE* pSrcFileLeft = NULL;pSrcFileLeft = fopen(srcFileLeft,"rb");if(pSrcFileLeft == NULL){printf("原始左片影像打开失败\n");return FAILURE;}FILE* pSrcFileRight = NULL;pSrcFileRight = fopen(srcFileRight,"rb");if(pSrcFileRight == NULL){printf("原始右片影像打开失败\n");return FAILURE;}//读取左片数据BITMAPFILEHEADER bmfhLeft;BITMAPINFOHEADER bmihLeft;fread(&bmfhLeft, sizeof(BITMAPFILEHEADER), 1, pSrcFileLeft);fread(&bmihLeft, sizeof(BITMAPINFOHEADER), 1, pSrcFileLeft);if(bmihLeft.biBitCount != 8){printf("只支持8位位图: 左片\n");return FAILURE;}DWORD widthLeft = bmihLeft.biWidth;DWORD heightLeft = bmihLeft.biHeight;DWORD _widthLeft = (widthLeft*8 + 31)/32*4; //saved widthLeft in files DWORD pixelCountLeft = _widthLeft*heightLeft;fseek(pSrcFileLeft, 256*sizeof(RGBQUAD), SEEK_CUR);BYTE* pSrcBitsLeft = new BYTE[pixelCountLeft]; //new 1fread(pSrcBitsLeft, sizeof(BYTE), pixelCountLeft, pSrcFileLeft);fclose(pSrcFileLeft);pSrcFileLeft = NULL;//读取右片数据BITMAPFILEHEADER bmfhRight;BITMAPINFOHEADER bmihRight;fread(&bmfhRight, sizeof(BITMAPFILEHEADER), 1, pSrcFileRight);fread(&bmihRight, sizeof(BITMAPINFOHEADER), 1, pSrcFileRight);if(bmihRight.biBitCount != 8){printf("只支持8位位图: 右片\n");return FAILURE;}DWORD widthRight = bmihRight.biWidth;DWORD heightRight = bmihRight.biHeight;DWORD _widthRight = (widthRight*8 + 31)/32*4; //saved widthRight in files DWORD pixelCountRight = _widthRight*heightRight;fseek(pSrcFileRight, 256*sizeof(RGBQUAD), SEEK_CUR);BYTE* pSrcBitsRight = new BYTE[pixelCountRight]; //new 2fread(pSrcBitsRight, sizeof(BYTE), pixelCountRight, pSrcFileRight);fclose(pSrcFileRight);pSrcFileRight = NULL;//开始计算左片特征点//计算左片的兴趣点存储于该vector中vector<FEATUREPOINT> FeaturePointLeft;vector<FEATUREPOINT> FeaturePointRight;int Count1(0);int Count2(0);Moravec(pSrcBitsLeft, &FeaturePointLeft, m_threshold_L, m_window1_L, m_window2_L, widthLeft, heightLeft, &Count1);Moravec(pSrcBitsRight, &FeaturePointRight, m_threshold_R, m_window1_R, m_window2_R,widthRight, heightRight, &Count2);cout << FeaturePointLeft.size() << endl;cout << FeaturePointRight.size() << endl;FeaturePointLeft = RemoveReplicative(FeaturePointLeft);FeaturePointRight = RemoveReplicative(FeaturePointRight);cout << "FeaturePointLeft.size() == " << FeaturePointLeft.size() << endl;cout << "FeaturePointRight.size() == " << FeaturePointRight.size() << endl;SaveList(FeaturePointLeft, "d:\\list1.txt");SaveList(FeaturePointRight, "d:\\list2.txt");//特征点计算完毕，存储于vector中//遍历vector, 在右片中搜寻相关系数最大的窗口vector<FEATUREPOINT>::iterator iter1;vector<FEATUREPOINT>::iterator iter2;int k = m_window3 / 2;int i(0), j(0), m(0), n(0);BYTE* pTemplate = new BYTE[m_window3*m_window3]; //new 5ZeroMemory(pTemplate, sizeof(BYTE)*m_window3*m_window3);BYTE* pTarget = new BYTE[m_window3*m_window3]; //new 6ZeroMemory(pTarget, sizeof(BYTE)*m_window3*m_window3);vector<MATCHEDPOINTS> MatchedPoints;float MaxR(0);float R(0);double A(0), B(0), C(0), D(0), E(0); //计算相关系数的中间量int lMaxWidth;int lMaxHeight;int Count3(0), Count4(0);for (iter1 = FeaturePointLeft.begin(); iter1 != FeaturePointLeft.end(); ++iter1){cout << endl << "Left point ID: " << Count2++ << endl;//计算该点对应的特征点B = 0.0;D = 0.0;MaxR = 0.0;R = 0.0;if (iter1->x < k || iter1->y < k || iter1->x > heightLeft - k || iter1->y > widthLeft - k) {cout << "Left point abort: " << Count3++ << endl;continue;}//逐个搜索//将左片小窗口内的数据拷贝//可先求出B,Dfor (i = 0; i < m_window3; ++i){for(j = 0; j < m_window3 ; ++j){//B = 小窗口所有像素灰度值求和//D = ....................的平方再求和B += pSrcBitsLeft[(i+iter1->y)*_widthLeft + (j+iter1->x)];D += pSrcBitsLeft[(i+iter1->y)*_widthLeft + (j+iter1->x)] *pSrcBitsLeft[(i+iter1->y)*_widthLeft + (j+iter1->x)];}}Count4 = 0;//逐个小窗口计算Rfor (iter2 = FeaturePointRight.begin(); iter2 != FeaturePointRight.end(); ++iter2) {if (iter2->x < k || iter2->y < k || iter2->x > heightRight - k || iter2->y > widthRight - k) {cout << "Right point abort: " << Count4++ << endl;continue;} A = 0.0;C = 0.0;E = 0.0;j = iter2->x; //i,j为中心点i = iter2->y;j -= k; //i,j为顶点i -= k;int centerX = j + k;int centerY = i + k;//这里的centerX，centerY就是右片某个小窗口的中心//求A,C,Efor (m = 0; m < m_window3; ++m){for(n = 0; n < m_window3 ; ++n){//C = 左片小窗口所有像素灰度值求和//E = ........................平方的求和//A = 对应位置像素灰度直乘积之和C += pSrcBitsRight[(i+m)*_widthRight + (j+n)];E += pSrcBitsRight[(i+m)*_widthRight+(j+n)]*pSrcBitsRight[(i+m)*_widthRight+(j+n)];A += pSrcBitsRight[(i+m)*_widthRight+(j+n)] *pSrcBitsLeft[(m+iter1->y)*_widthLeft + (n+iter1->x)];}}double temp1 = D - B*B/m_window3/m_window3;double temp2 = E - C*C/m_window3/m_window3;if (temp1 != 0 && temp2 != 0){R = (A - B*C/m_window3/m_window3) / sqrt(temp1*temp2);}elseR = 0;if (R > MaxR){MaxR = R;lMaxWidth = centerX;lMaxHeight = centerY;cout << R << endl;}}if (MaxR >= m_CoRelated){MATCHEDPOINTS mp;mp.coefficient = MaxR;mp.x1 = iter1->x;mp.y1 = iter1->y;mp.x2 = lMaxWidth;mp.y2 = lMaxHeight;MatchedPoints.push_back(mp);}}//保存结果位图DrawCross(pSrcBitsLeft, pSrcBitsRight, MatchedPoints,_widthLeft, _widthRight, pixelCountLeft, pixelCountRight);SaveBand(widthLeft, heightLeft, _widthLeft*heightLeft, 8, pSrcBitsLeft, dstFileLeft);SaveBand(widthRight, heightRight, _widthRight*heightRight, 8, pSrcBitsRight, dstFileRight);//保存特征点列表SaveList(MatchedPoints, txtFile);//clean updelete[] pSrcBitsLeft;delete[] pSrcBitsRight;delete[] pTemplate;pSrcBitsLeft = NULL;pSrcBitsRight = NULL;pTemplate = NULL;if (m_flag1 == 1){ShellExecute(NULL, "open", txtFile, NULL, NULL, SW_SHOWNORMAL);}return nRetCode;}。