c++实现影像构建多层金字塔提取特征点
发布时间
阅读量：
阅读量
c++实现影像构建多层金字塔提取特征点并还原到影像
文章目录
-
什么是影像金字塔？
-
如何构建影像金字塔？
-
通过影像金字塔分别提取特征点策略
-
二、使用步骤
-
- 方法一（基于ORB提取特征点）：
- 方法二（基于sift提取特征点）
- 方法三（基于sift提取特征点，这个更好，可以设置金字塔层数和尺度因子，但是目前数据量太多，内存溢出了，还需要优化，想要用的可以用）
-
总结
什么是影像金字塔？
影像金字塔由原始影像按一定规则生成的由细到粗不同分辨率的影像集。金字塔的底部是图像的高分辨率表示，也就是原始图像，而顶部是低分辨率的近似。最底层的分辨率最高，并且数据量最大，随着层数的增加，其分辨率逐渐降低，数据量也按比例减少。
如下图所示：【备注：此图从百度上下载的】
下图所示是通过一个1000*1000大小的单波段tif分别构建的0级、1级、2级、3级构建的四种级别金字塔显示。
如何构建影像金字塔？
可以通过pyrDown()和pryUp()两个函数分别进行下采样、升采样构建金字塔，根据个人需求进行使用，这里我们使用的是影像原图，因此只做下采样即可。
下采样：
void pyrDown(InputArray src, OutputArray dst, const Size& dstsize=Size());
上采样：
void pyrUp(InputArray src, OutputArray dst, const Size& dstsize=Size());
通过影像金字塔分别提取特征点策略
步骤一：
二、使用步骤
方法一（基于ORB提取特征点）：
这里我就一个主函数直接写了啊，可以直接用，根据需求改成sift、surf、harris等等均可；这里也给大家提供一个百分比截断拉伸算法（默认0.02、0.02即可，本人亲测效果最佳），之前网上的一位博主发的，找不到来源了，效果很好，大家拿去用吧。
代码如下（示例）：
main.cpp
/*tif创建金字塔*/
int main(int argc, const char* argv[])
{
#ifdef _DEBUG
printf("debug model");
const char* Image = "E:\ BaiduNetdiskDownload\ 1.tif";
string BaseMap = "E:\ BaiduNetdiskDownload\ 2.tif";
string outputpath = "E:\ BaiduNetdiskDownload\ 1.png";
#else
#endif
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO"); //支持中文路径
clock_t start = clock();
cout << "读取影像金字塔!" << endl;
string outputpath_parent = "";
outputpath_parent = parent(outputpath.c_str());
string outputpath_name = "";
outputpath_name = stem(outputpath.c_str());
待校正影像拉伸
string Image_8bit = "";
Image_8bit = outputpath_parent + "\ " + outputpath_name;
bool Linear_per2_s = -1;
Linear_per2_s = Linear_per2(Image, Image_8bit.c_str(), 0.02, 0.02);
cv::Mat img1 = cv::imread(Image_8bit, cv::IMREAD_GRAYSCALE);
cv::Mat img2 = cv::imread(BaseMap, cv::IMREAD_GRAYSCALE);
if (img1.empty() || img2.empty()) {
std::cerr << "Error loading images" << std::endl;
return -1;
}
// 构建金字塔 (示例中只构建了两层)
std::vector<cv::Mat> pyramid1, pyramid2;
pyramid1.push_back(img1);
pyramid2.push_back(img2);
for (int i = 0; i < 2; ++i) { // 创建两个级别的金字塔，这里根据需要构建金字塔级别数量进行修改
cv::Mat tmp1, tmp2;
pyrDown(pyramid1.back(), tmp1);
pyrDown(pyramid2.back(), tmp2);
pyramid1.push_back(tmp1);
pyramid2.push_back(tmp2);
}
// 提取特征点
cv::Ptr<cv::ORB> detector = cv::ORB::create();
std::vector<cv::KeyPoint> keypoints1, keypoints2;
cv::Mat descriptors1, descriptors2;
detector->detectAndCompute(pyramid1[1], cv::Mat(), keypoints1, descriptors1); // pyramid1[0]：第一层金字塔， pyramid1[1]：第二层金字塔，等等依次类推
detector->detectAndCompute(pyramid2[1], cv::Mat(), keypoints2, descriptors2);
// 特征点匹配
cv::BFMatcher matcher(cv::NORM_HAMMING);
std::vector<cv::DMatch> matches;
matcher.match(descriptors1, descriptors2, matches);
// 可视化匹配结果
cv::Mat imgMatches;
cv::drawMatches(pyramid1[1], keypoints1, pyramid2[1], keypoints2, matches, imgMatches);
cv::imwrite(outputpath, imgMatches);
clock_t end_match = clock();
double end1 = (double)(end_match - start) / CLOCKS_PER_SEC;
printf("处理花费时间: %f秒\n", end1);
return 0;
}
Linear_per2.cpp
/*tif创建金字塔*/
/*
百分比截断(16位降8位)
const char *inFilename：待拉伸影像
const char *dstFilename：拉伸后保存影像
double min_v：拉伸阈值_数据范围的后min_v的像素映射值（如0.02即后2%）
double max_v：拉伸阈值_数据范围的前min_v的像素映射值（如0.02即前2%）
备注：当min_v及max_v为均为0时，为最大最小值拉伸
*/
int Linear_per2(const char* inFilename, const char* dstFilename, double min_v, double max_v)
{
GDALAllRegister();
//为了支持中文路径
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");
int src_height = 0;
int src_width = 0;
string namefile = "";
namefile = name(inFilename);
GDALDataset* poIn = (GDALDataset*)GDALOpen(inFilename, GA_ReadOnly); //打开影像
//获取影像大小
src_width = poIn->GetRasterXSize();
src_height = poIn->GetRasterYSize();
//获取影像波段数
int InBands = poIn->GetRasterCount();
if (poIn == NULL)
{
printf("无法打开影像！\n");
return -1;
}
//影像遍历
long long int wh = (long long int)src_width * (long long int)src_height;
unsigned short* srcData_all = NULL;
srcData_all = (unsigned short*)malloc(sizeof(unsigned short) * wh);
memset(srcData_all, 0, sizeof(unsigned short) * wh);
poIn->GetRasterBand(1)->RasterIO(GF_Read, 0, 0, src_width, src_height, srcData_all, src_width, src_height, GDT_UInt16, 0, 0);
long long int a_b_1 = 0;
for(long long int i =0; i<wh; i++)
{
if (srcData_all[i]==0)
{
a_b_1 ++;
}
}
free(srcData_all);
srcData_all = NULL;
int a_b_2 = 0;
a_b_2 = (a_b_1 * 100) / wh;
cout << a_b_2 <<endl;
if (a_b_2 > 98 )
{
printf("待校正影像有效区域值小于99%！不进行拉伸处理！\n");
GDALClose((GDALDataset*)poIn);
return -1;
}
//获取影像格式
GDALDataType eDataType = poIn->GetRasterBand(1)->GetRasterDataType();
//定义存储影像的空间参考数组
double adfInGeoTransform[6] = { 0 };
const char* pszWKT = NULL;
//获取影像空间参考
poIn->GetGeoTransform(adfInGeoTransform);
pszWKT = poIn->GetProjectionRef();
//创建文件
GDALDriver* poDriver = (GDALDriver*)GDALGetDriverByName("GTiff");
if (InBands == 4)
{
InBands = 3;
}
GDALDataset* poOutputDS = poDriver->Create(dstFilename, src_width, src_height, InBands, GDT_Byte, NULL);
if (poOutputDS == NULL)
{
printf("创建输出影像失败！\n");
return -1;
}
//设置拉伸后图像的空间参考以及地理坐标
poOutputDS->SetGeoTransform(adfInGeoTransform);
poOutputDS->SetProjection(pszWKT);
//读取影像
//printf("%s\n 16位影像降到8位影像处理...\n", namefile.c_str());
for (int iBand = 0; iBand < InBands; iBand++)
{
//cout << "正在处理第 " << iBand + 1 << " 波段" << endl;
//读取影像
uint16_t* srcData = (uint16_t*)malloc(sizeof(uint16_t) * src_width * src_height * 1);
memset(srcData, 0, sizeof(uint16_t) * 1 * src_width * src_height);
int src_max = 0, src_min = 65500;
//读取多光谱影像到缓存
poIn->GetRasterBand(iBand + 1)->RasterIO(GF_Read, 0, 0, src_width, src_height, srcData + 0 * src_width * src_height, src_width, src_height, GDT_UInt16, 0, 0);
//}
//统计最大值
for (int src_row = 0; src_row < src_height; src_row++)
{
for (int src_col = 0; src_col < src_width; src_col++)
{
uint16_t src_temVal = *(srcData + src_row * src_width + src_col);
if (src_temVal > src_max)
src_max = src_temVal;
if (src_temVal < src_min)
src_min = src_temVal;
}
}
double* numb_pix = (double*)malloc(sizeof(double) * (src_max + 1)); //存像素值直方图，即每个像素值的个数
memset(numb_pix, 0, sizeof(double) * (src_max + 1));
//------- 统计像素值直方图 ------------//
for (int src_row = 0; src_row < src_height; src_row++)
{
for (int src_col = 0; src_col < src_width; src_col++)
{
uint16_t src_temVal = *(srcData + src_row * src_width + src_col);
*(numb_pix + src_temVal) += 1;
}
}
double* frequency_val = (double*)malloc(sizeof(double) * (src_max + 1)); //像素值出现的频率
memset(frequency_val, 0.0, sizeof(double) * (src_max + 1));
double frequency_value0;
for (int val_i = 0; val_i <= src_max; val_i++)
{
*(frequency_val + val_i) = *(numb_pix + val_i) / double(src_width * src_height);
if (val_i == 0)
{
frequency_value0 = *(frequency_val + val_i);
}
}
//cout<< sum<<endl;
double* accumlt_frequency_val = (double*)malloc(sizeof(double) * (src_max + 1)); //像素出现的累计频率
memset(accumlt_frequency_val, 0.0, sizeof(double) * (src_max + 1));
for (int val_i = 0; val_i <= src_max; val_i++)
{
for (int val_j = 0; val_j < val_i; val_j++)
{
*(accumlt_frequency_val + val_i) += *(frequency_val + val_j);
}
}
//统计像素两端截断值
int minVal = 0, maxVal = 0;
for (int val_i = 1; val_i < src_max; val_i++)
{
double acc_fre_temVal0 = *(frequency_val + 0);
double acc_fre_temVal = *(accumlt_frequency_val + val_i);
if ((acc_fre_temVal - acc_fre_temVal0) > min_v)
{
minVal = val_i;
break;
}
}
for (int val_i = src_max - 1; val_i > 0; val_i--)
{
double acc_fre_temVal0 = *(accumlt_frequency_val + src_max);
double acc_fre_temVal = *(accumlt_frequency_val + val_i);
if (acc_fre_temVal < (acc_fre_temVal0 - max_v))
{
maxVal = val_i;
break;
}
}
free(accumlt_frequency_val);
accumlt_frequency_val = NULL;
free(frequency_val);
frequency_val = NULL;
for (int src_row = 0; src_row < src_height; src_row++)
{
uint8_t* dstData = (uint8_t*)malloc(sizeof(uint8_t) * src_width);
memset(dstData, 0, sizeof(uint8_t) * src_width);
for (int src_col = 0; src_col < src_width; src_col++)
{
uint16_t src_temVal = *(srcData + src_row * src_width + src_col);
double stre_temVal = (src_temVal - minVal) / double(maxVal - minVal);
if (src_temVal < minVal)
{
if (src_temVal == 0)
{
*(dstData + src_col) = 0;
}
else
{
*(dstData + src_col) = 1;//原来不是0的如果直接拉伸0将变为nodata
}
//*(dstData + src_col) = (src_temVal) *(20.0/double(minVal)) ;
//*(dstData + src_col) = 0;
}
else if (src_temVal > maxVal)
{
//stre_temVal = (src_temVal - src_min) / double(src_max - src_min);
*(dstData + src_col) = 255;
}
else
{
//*(dstData + src_col) = pow(stre_temVal,0.7) * 250;
//*(dstData + src_col) = (src_temVal - minVal)/(double)(maxVal-minVal) * 255.0;
int tmp = (src_temVal - minVal) / (double)(maxVal - minVal) * 255.0;
if (tmp <= 0)
{
if (src_temVal == 0)
{
*(dstData + src_col) = 0;
}
else
{
*(dstData + src_col) = 1;//原来不是0的如果直接拉伸0将变为nodata
//cout<<src_temVal<<","<<*(dstData + src_col)<<endl;
}
}
else if (tmp >= 255)
{
*(dstData + src_col) = 255;
}
else
{
*(dstData + src_col) = (int)tmp;
}
}
}
poOutputDS->GetRasterBand(iBand + 1)->RasterIO(GF_Write, 0, src_row, src_width, 1, dstData, src_width, 1, GDT_Byte, 0, 0);
free(dstData);
dstData = NULL;
}
//设置无效数据nodata为1
//poOutputDS->GetRasterBand(iBand + 1)->SetNoDataValue(1);
free(numb_pix);
numb_pix = NULL;
free(srcData);
srcData = NULL;
}
GDALClose((GDALDataset*)poIn);
GDALClose((GDALDataset*)poOutputDS);
return 0;
}
方法二（基于sift提取特征点）
这里我把这一部分定义成一个函数，直接使用即可；
代码如下（示例）：
bool CreatePyramid2(const char* bzFile, const char* pzFile, vector<KeyPoint>& keypoints1, vector<KeyPoint>& keypoints2, Mat& descriptors1, Mat& descriptors2)
{
cv::Mat img1 = cv::imread(bzFile, cv::IMREAD_GRAYSCALE);
cv::Mat img2 = cv::imread(pzFile, cv::IMREAD_GRAYSCALE);
cout<< bzFile <<endl;
cout << pzFile << endl;
if (img1.empty() || img2.empty()) {
std::cerr << "Error loading images" << std::endl;
return false;
}
cout << "这里没错误！" << endl;
// 构建金字塔 (示例中只构建了两层)
std::vector<cv::Mat> pyramid1, pyramid2;
pyramid1.push_back(img1);
pyramid2.push_back(img2);
for (int i = 0; i < 2; ++i) { // 创建两个级别的金字塔，这里根据需要构建金字塔级别数量进行修改
cv::Mat tmp1, tmp2;
pyrDown(pyramid1.back(), tmp1);
pyrDown(pyramid2.back(), tmp2);
pyramid1.push_back(tmp1);
pyramid2.push_back(tmp2);
}
// 提取特征点
Ptr<Feature2D> detector = cv::SIFT::create();
//detector->detect(img_1, keypoints_1);
//cv::Ptr<cv::sift> detector = cv::ORB::create();
//std::vector<cv::KeyPoint> keypoints1, keypoints2;
//cv::Mat descriptors1, descriptors2;
detector->detectAndCompute(pyramid1[0], cv::Mat(), keypoints1, descriptors1); // // pyramid1[0]：第一层金字塔， pyramid1[1]：第二层金字塔，等等依次类推
detector->detectAndCompute(pyramid2[0], cv::Mat(), keypoints2, descriptors2);
/*detector->detect(pyramid1[1], keypoints1);
detector->detect(pyramid2[1], keypoints2);*/
// 特征点匹配
/*cv::BFMatcher matcher(cv::NORM_HAMMING);
std::vector<cv::DMatch> matches;
matcher.match(descriptors1, descriptors2, matches);*/
return true;
}
方法三（基于sift提取特征点，这个更好，可以设置金字塔层数和尺度因子，但是目前数据量太多，内存溢出了，还需要优化，想要用的可以用）
这里我把这一部分定义成一个函数，直接使用即可；
代码如下（示例）：
/*该方法需要优化，一开始提点合并导致数组越界，需要改为每个金字塔的特征点分别提取和匹配筛选后在还原到原图输出*/
bool CreatePyramid(const char* bzFile, vector<KeyPoint>& keypoint_s, Mat& descriptor_Mat)
{
// 加载TIF图像
Mat img = imread(bzFile, IMREAD_GRAYSCALE);
if (img.empty()) {
cout << "图像加载失败" << endl;
return false;
}
// 创建SIFT特征检测器
Ptr<SIFT> sift = SIFT::create();
// 定义金字塔层数和尺度因子
int pyramidLevels = 4; // 金字塔的层数
double scaleFactor = 1.5; // 每层之间的尺度因子
// 创建图像金字塔
//vector<Mat> pyramid; //构建金字塔数组，可以通过构建多个金字塔分别待校正影像和基准影像匹配
Mat currentImg = img;
/*vector<KeyPoint> keypoint_s*/;
vector<Mat> descriptor_s;
for (int level = 0; level < pyramidLevels; ++level) {
// 在当前层上提取特征点
vector<KeyPoint> keypoints;
Mat descriptors;
sift->detectAndCompute(currentImg, noArray(), keypoints, descriptors);
cout << "在层级 " << level << " 上找到的特征点数量: " << keypoints.size() << endl;
// 将当前图像加入到金字塔中
//pyramid.push_back(currentImg);
// 准备下一层图像
Mat resizedImg;
resize(currentImg, resizedImg, Size(), 1.0 / scaleFactor, 1.0 / scaleFactor);
currentImg = resizedImg;
if (level != 0)
{
float scale = pow(2, level); // pyramid_level是控制点所在的金字塔层级
for (auto& kp : keypoints) {
kp.pt.x *= scale;
kp.pt.y *= scale;
}
}
for (int i = 0; i < keypoints.size(); i++)
{
keypoint_s.push_back(keypoints[i]);
descriptor_s.push_back(descriptors);
}
keypoints.clear();
keypoints.shrink_to_fit();
resizedImg.release();
descriptors.release();
}
//多个相同的尺寸和数据类型的Mat合并到一个Mat
// 步骤2: 计算合并后的矩阵需要多少行
long long int totalRows = 0;
for (const auto& descriptor : descriptor_s) {
totalRows += descriptor.rows;
}
// 假设所有描述子矩阵都有相同的列数和类型
if (descriptor_s.empty())
{
return false; // 确保descriptors不为空
}
int cols = descriptor_s[0].cols;
int type = descriptor_s[0].type();
// 步骤3: 创建合并后的描述子矩阵
cv::Mat mergedDescriptors(totalRows, cols, type);
// 步骤4: 复制每个矩阵到合并矩阵中
long long int currentRow = 0;
for (const auto& descriptor : descriptor_s) {
descriptor.copyTo(mergedDescriptors.rowRange(currentRow, currentRow + descriptor.rows));
currentRow += descriptor.rows;
}
descriptor_Mat = mergedDescriptors.clone(); // 使用clone()来复制数据
currentImg.release();
for (int i = 0; i < descriptor_s.size(); i++)
{
descriptor_s[i].release();
}
descriptor_s.clear();
descriptor_s.shrink_to_fit();
/*string tifoutput_a = "";
for (int i=0; i< pyramid.size(); i++)
{
tifoutput_a = tifoutput + "\ " + to_string(i + 1) + ".tif";
cv::imwrite(tifoutput_a, pyramid[i]);
}*/
return true;
}
总结
金字塔构建和匹配这一块，如果直接百度查询，发现很多都是零零散散的，对这一块不是很熟悉的话，不是很好写，东一块西一块，思路很简单，但是实现总是比较棘手，在这里我也是以明确的思路，结合网上的一些文献和GPT4.0进行查询逐步实现。GPT4.0确实好用，但也要结合实际需求，很多难度系数大的给的都是错误答案，在使用需要拆解很细，可以得到想要的答案。思路是重点！！！！
全部评论 (0)
还没有任何评论哟~