SIFT特征提取C++

#include "SiftDetector.h" using std::vector; #define INTERVALS 2 #define SIGMA_ANTIALIAS 0.5 #define SIGMA_PREBLUR 1.0 #define EDGE_THRESHOLD 10.0 #define CONTRAST_THRESHOLD 0.03 #define FEATURE_NUM_BINS 36 #define MAX_KERNEL_SIZE 20 #define M_PI 3.14159265358979 #define KEYPOINT_WINDOW_SIZE 16 #define DESCRIPTOR_NUM_BINS 8 #define DESCRIPTOR_THRESHOLD 0.2 void SaveFloatingPointImage(const char *filename, IplImage *img) { IplImage *dup = cvCreateImage(cvGetSize(img), 8, 1); cvCvtScale(img, dup, 255.0); cvSaveImage(filename, dup); cvReleaseImage(&dup); } SiftDetector::SiftDetector(IplImage* img, int octaves) { m_srcImage = cvCloneImage(img); m_numOctaves = octaves; InitializeImage(); } SiftDetector::SiftDetector(const char* filename, int octaves) { m_srcImage = cvLoadImage(filename); m_numOctaves = octaves; InitializeImage(); } void SiftDetector::InitializeImage() { int i=0; // step 1 创建一个包含不同高斯模糊图像的2D数组 m_scaleSpace = new IplImage**[m_numOctaves]; for (i=0; i<m_numOctaves; i++) m_scaleSpace[i] = new IplImage*[INTERVALS+3]; // step 2 创建一个包含不同的DoG图像的2D数组 m_DoGImageList = new IplImage**[m_numOctaves]; for (i=0; i<m_numOctaves; i++) m_DoGImageList[i] = new IplImage*[INTERVALS+2]; // step 1 创建一个包含图像的模糊量的2D数组 m_Sigma = new double*[m_numOctaves]; for (i=0; i<m_numOctaves; i++) m_Sigma[i] = new double[INTERVALS+3]; } SiftDetector::~SiftDetector() { int i,j; for (i=0; i<m_numOctaves; i++) { // 释放某组尺度空间中的图像 for (j=0; j<INTERVALS+3; j++) cvReleaseImage(&m_scaleSpace[i][j]); for (j=0; j<INTERVALS+2; j++) cvReleaseImage(&m_DoGImageList[i][j]); // 释放数组内存 delete[] m_scaleSpace[i]; delete[] m_DoGImageList[i]; delete[] m_Sigma[i]; } for (i=0; i<m_numKeypoints; i++) { cvReleaseImage(&m_magnitude[i]); cvReleaseImage(&m_orientation[i]); } // 释放2D数组内存 delete[] m_scaleSpace; delete[] m_DoGImageList; delete[] m_Sigma; delete[] m_magnitude; delete[] m_orientation; } void SiftDetector::ShowKeypoints() { IplImage* img = cvCloneImage(m_srcImage); for (int i=0; i<m_numKeypoints; i++) { Keypoint& kp = m_extrema[i]; double scale = pow(2.0, kp.octave-1); int kr = kp.kr * scale; int kc = kp.kc * scale; double length = 10 * kp.maxMag; double orien = kp.maxOri; cvLine(img, cvPoint(kc, kr), cvPoint(kc, kr), CV_RGB(255, 255, 255), 3); cvLine(img, cvPoint(kc, kr), cvPoint(kc+length*cos(orien), kr+length*sin(orien)), CV_RGB(255, 255, 255), 1); } cvNamedWindow("Keypoints"); cvShowImage("Keypoints", img); } void SiftDetector::DoSiftDetect() { BuildScaleSpace(); DetectExtrema(); AssignOrientations(); ExtractKeypointDescriptors(); } void SiftDetector::BuildScaleSpace() { printf("Generating scale space...\n"); IplImage* imgGray = cvCreateImage(cvGetSize(m_srcImage), 32, 1); IplImage* imgTemp = cvCreateImage(cvGetSize(m_srcImage), 8, 1); // 如果图像是彩色的，把它转化为黑白的 if (m_srcImage->nChannels == 3) cvCvtColor(m_srcImage, imgTemp, CV_BGR2GRAY); else cvCopy(m_srcImage, imgTemp); // 生成浮点图像，将[0, 255]转变为[0.0, 1.0] for (int c=0; c<imgTemp->width; c++) for (int r=0; r<imgTemp->height; r++) cvSetReal2D(imgGray, r, c, cvGetReal2D(imgTemp, r, c)/255.0); // Lowe给出用0.5的模糊量来模糊图像，并放大2倍尺寸能够增加稳定特征点的个数 cvSmooth(imgGray, imgGray, CV_GAUSSIAN, 0, 0, SIGMA_ANTIALIAS); m_scaleSpace[0][0] = cvCreateImage(cvSize(imgGray->width*2, imgGray->height*2), 32, 1); cvPyrUp(imgGray, m_scaleSpace[0][0]); // 模糊第一幅图像 cvSmooth(m_scaleSpace[0][0], m_scaleSpace[0][0], CV_GAUSSIAN, 0, 0, SIGMA_PREBLUR); double initSigma = sqrt(2.0f); m_Sigma[0][0] = initSigma*0.5; for (int i=0; i<m_numOctaves; i++) { double sigma = initSigma; CvSize currentSize = cvGetSize(m_scaleSpace[i][0]); //char filename[200]; for (int j=1; j<INTERVALS+3; j++) { m_scaleSpace[i][j] = cvCreateImage(currentSize, 32, 1); // step 1 模糊当前图像以得到下一幅图像，更新模糊量sigma cvSmooth(m_scaleSpace[i][j-1], m_scaleSpace[i][j], CV_GAUSSIAN, 0, 0, sigma); sigma = sqrt(2.0f)*sigma; m_Sigma[i][j] = sqrt(2.0f)*m_Sigma[i][j-1]; // step 2 计算Difference of Gaussian图像 m_DoGImageList[i][j-1] = cvCreateImage(currentSize, 32, 1); cvSub(m_scaleSpace[i][j-1], m_scaleSpace[i][j], m_DoGImageList[i][j-1]); // 保存生成的图像for debug /*sprintf(filename, "images\\Gaussian\\g_octave_%d_scale_%d.jpg", i, j); SaveFloatingPointImage(filename, m_scaleSpace[i][j]); sprintf(filename, "images\\DoG\\dog_octave_%d_scale_%d.jpg", i, j); SaveFloatingPointImage(filename, m_DoGImageList[i][j-1]);*/ } if (i<m_numOctaves-1) { // 减少大小 currentSize.width/=2; currentSize.height/=2; // step 1 分配内存，创建下一组尺度空间的第一幅图像 m_scaleSpace[i+1][0] = cvCreateImage(currentSize, 32, 1); cvPyrDown(m_scaleSpace[i][0], m_scaleSpace[i+1][0]); m_Sigma[i+1][0] = m_Sigma[i][INTERVALS]; /*sprintf(filename, "images\\Gaussian\\g_octave_%d_scale_%d.jpg", i+1, 0); SaveFloatingPointImage(filename, m_scaleSpace[i+1][0]);*/ } } } void SiftDetector::DetectExtrema() { printf("Detecting extrema..\n"); int num = 0; int numRemoved = 0; double edge_threshold = (EDGE_THRESHOLD+1)*(EDGE_THRESHOLD+1)/EDGE_THRESHOLD; for (int i=0; i<m_numOctaves; i++) { CvSize currentSize = cvGetSize(m_DoGImageList[i][0]); for (int j=1; j<INTERVALS+1; j++) { IplImage* middle = m_DoGImageList[i][j]; IplImage* up = m_DoGImageList[i][j+1]; IplImage* down = m_DoGImageList[i][j-1]; for (int c=1; c<middle->width-1; c++) for (int r=1; r<middle->height-1; r++) { bool isExtrema = false; double currentPixel = cvGetReal2D(middle, r, c); if ( // step 3 检测是否是极大值 currentPixel > cvGetReal2D(middle, r-1, c-1) && currentPixel > cvGetReal2D(middle, r-1, c ) && currentPixel > cvGetReal2D(middle, r-1, c+1) && currentPixel > cvGetReal2D(middle, r, c-1) && currentPixel > cvGetReal2D(middle, r, c+1) && currentPixel > cvGetReal2D(middle, r+1, c-1) && currentPixel > cvGetReal2D(middle, r+1, c ) && currentPixel > cvGetReal2D(middle, r+1, c+1) && currentPixel > cvGetReal2D(up, r-1, c-1) && currentPixel > cvGetReal2D(up, r-1, c ) && currentPixel > cvGetReal2D(up, r-1, c+1) && currentPixel > cvGetReal2D(up, r, c-1) && currentPixel > cvGetReal2D(up, r, c+1) && currentPixel > cvGetReal2D(up, r+1, c-1) && currentPixel > cvGetReal2D(up, r+1, c ) && currentPixel > cvGetReal2D(up, r+1, c+1) && currentPixel > cvGetReal2D(down, r-1, c-1) && currentPixel > cvGetReal2D(down, r-1, c ) && currentPixel > cvGetReal2D(down, r-1, c+1) && currentPixel > cvGetReal2D(down, r, c-1) && currentPixel > cvGetReal2D(down, r, c+1) && currentPixel > cvGetReal2D(down, r+1, c-1) && currentPixel > cvGetReal2D(down, r+1, c ) && currentPixel > cvGetReal2D(down, r+1, c+1) || // step 3 检测是否是极小值 currentPixel < cvGetReal2D(middle, r-1, c-1) && currentPixel < cvGetReal2D(middle, r-1, c ) && currentPixel < cvGetReal2D(middle, r-1, c+1) && currentPixel < cvGetReal2D(middle, r, c-1) && currentPixel < cvGetReal2D(middle, r, c+1) && currentPixel < cvGetReal2D(middle, r+1, c-1) && currentPixel < cvGetReal2D(middle, r+1, c ) && currentPixel < cvGetReal2D(middle, r+1, c+1) && currentPixel < cvGetReal2D(up, r-1, c-1) && currentPixel < cvGetReal2D(up, r-1, c ) && currentPixel < cvGetReal2D(up, r-1, c+1) && currentPixel < cvGetReal2D(up, r, c-1) &&