基于相关系数影像匹配实习报告(共15页)

1、精选优质文档-倾情为你奉上基于相关系数影像匹配实习报告1、 实习内容与目的 通过上机编程实现基于相关系数算法的影像匹配，从而更加熟练的掌握书本上关于影像匹配及相关系数算法的理论知识，并熟悉上机编程的操作。2、 实习原理影像匹配实质上是在两幅或是多幅影像之间识别同名点，其中基于相关系数的算法是实现影像匹配的基本算法之一。相关系数是标准化的协方差函数，协方差函数除以两信号的方差即得相关系数。计算相关系数的公式如下：g(x,y)与g(x,y)的相关系数为：其中：若(p0, q0)>(p, q)，(pp0,qq0)，则p0,q0为搜索区影像相对于目标区影像的位移参数。由离散灰度数据对相关系数的估

2、计为，考虑到计算工作量，相关系数的实用公式为3、 实习步骤及相关代码 在第一次点特征点实习的基础上我们组进行了第二次基于相关系数的影像匹配实习。所以程序实现的部分操作直接利用前次实习的中的程序。（1） 读入左右影像 FILE* pSrcFileLeft = NULL;pSrcFileLeft = fopen(srcFileLeft,"rb");if(pSrcFileLeft = NULL)printf("原始左片影像打开失败n");return FAILURE;FILE* pSrcFileRight = NULL;pSrcFileRight = fope

3、n(srcFileRight,"rb");if(pSrcFileRight = NULL)printf("原始右片影像打开失败n");return FAILURE;（2） 确定目标窗口大小 通过老师课堂上的讲授和我们在编程过程中的不断尝试，我们默认使用的窗口的尺寸为11*11。（3） 确定目标点的位置 由前次的实习中得到左相片中的特征点作为这次实习的目标点。（4） 预测右影像搜索范围 本程序中也在左相片中计算出了特征点，然后将左右两相片中的特征点进行计算相关系数。（5） 逐窗口计算相关系数并保存相应的储存单元 将左相片中的目标点与右相片中的特征点进行相关

4、系数的计算。（6） 比较相关系数，取最大相关系数作为匹配点 对于分别比较上一步骤中计算出的相关系数，选取出最大的一个相关系数作为匹配点。 程序操作如下：在程序主界面中点击Extend中选取ImageMatch选项，然后进入参数设置界面。设置读取和保存影像路径和窗口大小等参数后，点击Progress进行运算。4、 实习结果及分析 作为数据源的左右影像如些图所示：左相片 右相片 经过计算后得到的目标点左右图像分别如下所示：左相片 右相片 显示匹配点号： 通过老师得到的图片可以看出在左相片中得到的五个目标点中有四个在右相片中找到了同名点，由此可以看出本程序基本满足实习要求，并得到理想的结果。不足的地

5、方就是选取的特征点的数量不够多，这样就无法判断本程序的“误判”的情况如何。5、 实习体会 在学习到相关系数的时候，一看到计算相关系数的公式的时候我就有种被折服的感觉，那个公式实在是太复杂，它使我感觉我怎么也不能记住这个复杂的公式。不过经过这次实习后，我觉得这个公式又不是当初认为的那样无法记忆。我想这次实习我最大的收获就是对于之前在学习过程中遇到不懂知识以及难以记忆的知识都得到了加强，方便了我对他们的理解与记忆。 其次，在编程过程中更加锻炼自己动手编程的能力，使得我编程能力较之前有着一定的提高。每次编程实习都是一个积累经验的过程，然后通过一次次的实习后，我最终的能力就能够得到很好的体现。附件：v

6、oid 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<FEATURE

7、POINT>: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

8、 = 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.bfR

9、eserved2 = 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

10、 = sizeof(BITMAPFILEHEADER) + sizeof(BITMAPINFOHEADER) + byteCount;elseMessageBox(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;bmi

11、h.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 rgb256;int i(0);for

12、 (; i < 256; +i)rgbi.rgbBlue = i;rgbi.rgbGreen = i;rgbi.rgbRed = i;rgbi.rgbReserved = 0;/修改颜色表255对应的颜色rgb255.rgbBlue = 0;rgb255.rgbGreen = 0;rgb255.rgbRed = 255;f.Write(rgb, sizeof(RGBQUAD)*256);f.Write(pBits, byteCount);f.Close();elseMessageBox(NULL, "保存图像失败", "错误", MB_ICONER

13、ROR);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 (pBitsi = 255)pBitsi = 254;for (iter = FeaturePoint.begin(); iter != FeaturePoint.end(); +iter)i = it

14、er->y;j = iter->x;pBitsi*_width+j = 255;pBitsi*_width+j+1 = 255;pBitsi*_width+j-1 = 255;pBits(i+1)*_width+j = 255;pBits(i-1)*_width+j = 255;pBitsi*_width+j+2 = 255;pBitsi*_width+j-2 = 255;pBits(i+2)*_width+j = 255;pBits(i-2)*_width+j = 255;void DrawCross(BYTE* pBits1, BYTE* pBits2, vector<M

15、ATCHEDPOINTS> 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 (pBits1i = 255)pBits1i = 254;for (j = 0; j < byteCount2; +j)if (pBits2j = 255)pBits2j = 254;for (iter = MatchedPoi

16、nts.begin(); iter != MatchedPoints.end(); +iter)i = iter->y1;j = iter->x1;pBits1i*_width1+j = 255;pBits1i*_width1+j+1 = 255;pBits1i*_width1+j-1 = 255;pBits1(i+1)*_width1+j = 255;pBits1(i-1)*_width1+j = 255;pBits1i*_width1+j+2 = 255;pBits1i*_width1+j-2 = 255;pBits1(i+2)*_width1+j = 255;pBits1(i

17、-2)*_width1+j = 255;i = iter->y2;j = iter->x2;pBits2i*_width2+j = 255;/pBits2i*_width2+j+1 = 255;pBits2i*_width2+j-1 = 255;pBits2(i+1)*_width2+j = 255;pBits2(i-1)*_width2+j = 255;pBits2i*_width2+j+2 = 255;pBits2i*_width2+j-2 = 255;pBits2(i+2)*_width2+j = 255;pBits2(i-2)*_width2+j = 255;void Sa

18、veList(vector<FEATUREPOINT> FeaturePoint, CString txtFile)FILE* pTxtFile= fopen(txtFile, "wt");if (!pTxtFile)return;fprintf(pTxtFile, "#t特征点数: %dnn", FeaturePoint.size();fprintf(pTxtFile, "坐标Xt坐标Yt兴趣值n");vector<FEATUREPOINT>:iterator iter;for (iter = Feature

19、Point.begin(); iter != FeaturePoint.end(); +iter)fprintf(pTxtFile, "%dt%dt%dn", 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(pT

20、xtFile, "#同名点对数: %dnn",MatchedPoints.size();vector<MATCHEDPOINTS>:iterator iter;for (iter = MatchedPoints.begin(); iter != MatchedPoints.end(); +iter)fprintf(pTxtFile, "%dt%dt%dt%dt%fn",iter->x1, iter->y1, iter->x2, iter->y2, iter->coefficient);fclose(pTxtFil

21、e);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(

22、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_Rig

23、ht,window2_Right;CString CoRelated,window3,flag1;/启动参数srcFileLeft = argv1;srcFileRight = argv2;dstFileLeft = argv3;dstFileRight = argv4;txtFile = argv5;threshold_Left = argv6;/阈值window1_Left = argv7;window2_Left = argv8;threshold_Right = argv9;window1_Right = argv10;window2_Right = argv11;CoRelated

24、= argv12;window3 = argv13;flag1 = argv14;/判断是否在创建结果文件后显示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_Co

25、Related = (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)

26、= -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* pSrcFile

27、Right = 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(BITMAPINFOHEAD

28、ER), 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 filesDWORD pixelCountLeft = _widthLeft*heightLeft;fseek(pSrcFileLef

29、t, 256*sizeof(RGBQUAD), SEEK_CUR);BYTE* pSrcBitsLeft = new BYTEpixelCountLeft;/new 1fread(pSrcBitsLeft, sizeof(BYTE), pixelCountLeft, pSrcFileLeft);fclose(pSrcFileLeft);pSrcFileLeft = NULL;/读取右片数据BITMAPFILEHEADER bmfhRight;BITMAPINFOHEADER bmihRight;fread(&bmfhRight, sizeof(BITMAPFILEHEADER), 1,

30、 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

31、 filesDWORD pixelCountRight = _widthRight*heightRight;fseek(pSrcFileRight, 256*sizeof(RGBQUAD), SEEK_CUR);BYTE* pSrcBitsRight = new BYTEpixelCountRight;/new 2fread(pSrcBitsRight, sizeof(BYTE), pixelCountRight, pSrcFileRight);fclose(pSrcFileRight);pSrcFileRight = NULL;/开始计算左片特征点/计算左片的兴趣点存储于该vector中ve

32、ctor<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_

33、window1_R, m_window2_R,widthRight, heightRight, &Count2);cout << FeaturePointLeft.size() << endl;cout << FeaturePointRight.size() << endl;FeaturePointLeft = RemoveReplicative(FeaturePointLeft);FeaturePointRight = RemoveReplicative(FeaturePointRight);cout << "Fe

34、aturePointLeft.size() = " << FeaturePointLeft.size() << endl;cout << "FeaturePointRight.size() = " << FeaturePointRight.size() << endl;SaveList(FeaturePointLeft, "d:list1.txt");SaveList(FeaturePointRight, "d:list2.txt");/特征点计算完毕，存储于vec

35、tor中/遍历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 BYTEm_window3*m_window3;/new 5ZeroMemory(pTemplate, sizeof(BYTE)*m_window3*m_window3);BYTE* pTarget = new BYTEm_win

36、dow3*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 != FeaturePo

37、intLeft.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 abo

38、rt: " << 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) * pS

39、rcBitsLeft(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