温馨提示：这篇文章已超过469天没有更新，请注意相关的内容是否还可用！
摘要：，，在QT框架中，利用GDAL（Geospatial Data Abstraction Library）实现多线程读取遥感图像至QImage。通过GDAL库，能够高效处理地理空间数据，结合QT的多线程功能，可并行读取多张遥感图像，显著提高数据处理效率。此过程涉及图像数据的解析、转换和加载，最终将遥感图像数据成功加载至QImage中，便于后续处理和分析。
GDAL 是一个很强大的可以读取很多格式 的带有GIS信息的栅格型图像。前阵子项目中需要读取遥感图像，并显示到QT界面，由于遥感图像一般很大，所以采取新开一个读图线程的方式来读取，防止界面假死。下面是代码共享，测试通过读取500MB的24000*24000像素GeoTiff图并在QT的QGraphicsView中显示。
（图片来源网络，侵删）
环境：VS2005+SP1, Qt 4.6.0, GDAL 1.6.2
文件：commontoolfunctions.h, gdalimagereaderthread.h, gdalimagereaderthread.cpp
（图片来源网络，侵删）
1. commontoolfunctions.h文件：
/*
* 工具函数：网上找来的根据两点的经纬度坐标计算两个点的大地距离。
*/
#ifndef COMMONTOOLFUNCTIONS_H
#define COMMONTOOLFUNCTIONS_H
#include
#define EARTH_RADIUS 6378.137
#define MY_PI 3.14159265358979
#define LB_CONVERT_TO_RAD(d) (d * MY_PI / 180.0)
/*
* 输入经纬度坐标，单位'度'。
* 返回值单位为公里。
* code from google maps. 近似值.
*/
double CalcDistanceByLB(double lat1, double lng1, double lat2, double lng2)
{
double radLat1 = LB_CONVERT_TO_RAD(lat1);
double radLat2 = LB_CONVERT_TO_RAD(lat2);
double a = radLat1 - radLat2;
double b = LB_CONVERT_TO_RAD(lng1) - LB_CONVERT_TO_RAD(lng2);
double s = 2 * asin(sqrt(pow(sin(a/2),2) + cos(radLat1)*cos(radLat2)*pow(sin(b/2),2)));
s = s * EARTH_RADIUS;
return s;
}
#endif
2. gdalimagereaderthread.h文件
/*
* QT线程类：使用GDAL读取图像文件并转为QImage的线程。
* 创建：MulinB
* 日期：2010-1-8
* 上次修改: 2010-06-20
* HUST-IPRAI-CVL
*/
#ifndef GDALIMAGEREADERTHREAD_H
#define GDALIMAGEREADERTHREAD_H
#include
#include //for GDAL
class QImage;
class RSI_MapCanvas;
class RSI_ImageReaderThread : public QThread
{
Q_OBJECT
public:
RSI_ImageReaderThread(RSI_MapCanvas* pCanvas);
~RSI_ImageReaderThread();
public:
void ReadFile(const QString& imageFileName); //使用该函数读取文件
protected:
void run(); //线程函数
private:
RSI_MapCanvas* m_pCanvas;
QString m_imageFileName;
private:
//GDAL读取图像文件时的临时变量
GDALRasterBand* m_pBand_GDAL_;
int m_nXSize_GDAL_;
int m_nYSize_GDAL_;
void* m_pDataBufStart_GDAL_;
int m_nBufXSize_GDAL_;
int m_nBufYSize_GDAL_;
GDALDataType m_eBufDataType_GDAL_;
int m_nBufPixSpace_GDAL_;
int m_nBufLineSpace_GDAL_;
private:
enum IMG_OPEN_DRIVER_T
{
RSI_IMG_DRIVER_UNDEF = 0, //未定义的打开方式
RSI_USE_OPENCV, //使用opencv
RSI_USE_GDAL, //使用GDAL
};
private:
//将图像宽度字节数32位对齐
inline int Get32bitAlignedWidthBytes(int nImgPixWidth, int nBytesPerPix);
//GDAL buffer 转 QImage
QImage* CvtGDALBufferToGrayQImage(uchar* &pafScanblock, int nImgSizeX, int nImgSizeY, GDALDataType pixType);
//将非8位的灰度图转换为8位灰度图
template
void CvtImgDataTo8bitGrayImg(T* pOrigData, int nImgSizeX, int nImgSizeXAligned32bit, int nImgSizeY, uchar* pNew8bitImgBuf, int nNewImgSizeXAligned32bit, T maxVal, T minVal);
};
#endif // GDALIMAGEREADERTHREAD_H 3. gdalimagereaderthread.cpp文件
/*
* QT线程类：使用GDAL读取图像文件并转为QImage的线程。
* 创建：MulinB
* 日期：2010-1-8
* 上次修改: 2010-06-20
* TODO: 1.图像分辨率转换计算貌似有问题...
2.代码比较乱，没时间整理...囧...
* HUST-IPRAI-CVL
*/
#include "GDALImageReaderThread.h"
#include "iplimagetoqimage.h" //for IplImage converting to QImage
#include "commontoolfunctions.h" //for calculate distance by LB
#include //for OpenCV
#include //for OpenCV
#include //for INT_MAX, INT_MIN, ...
#include //for FLT_MAX, FLT_MIN, ...
QMutex G_imageReadMutex;
QWaitCondition G_imageReadComplete;
RSI_ImageReaderThread::RSI_ImageReaderThread(RSI_MapCanvas* pCanvas)
: m_pCanvas(pCanvas)
{
}
RSI_ImageReaderThread::~RSI_ImageReaderThread()
{
}
//线程函数
void RSI_ImageReaderThread::run()
{
//这里只将GDAL最费时的RasterIO()函数作为新线程
/*
GDALRasterBand::RasterIO() parameters:
eRWFlag Either GF_Read to read a region of data, or GT_Write to write a region of data.
nXOff The pixel offset to the top left corner of the region of the band to be accessed. This would be zero to start from the left side.
nYOff The line offset to the top left corner of the region of the band to be accessed. This would be zero to start from the top.
nXSize The width of the region of the band to be accessed in pixels.
nYSize The height of the region of the band to be accessed in lines.
pData The buffer into which the data should be read, or from which it should be written. This buffer must contain at least nBufXSize * nBufYSize words of type eBufType. It is organized in left to right, top to bottom pixel order. Spacing is controlled by the nPixelSpace, and nLineSpace parameters.
nBufXSize the width of the buffer image into which the desired region is to be read, or from which it is to be written.
nBufYSize the height of the buffer image into which the desired region is to be read, or from which it is to be written.
eBufType the type of the pixel values in the pData data buffer. The pixel values will automatically be translated to/from the GDALRasterBand data type as needed.
nPixelSpace The byte offset from the start of one pixel value in pData to the start of the next pixel value within a scanline. If defaulted (0) the size of the datatype eBufType is used.
nLineSpace The byte offset from the start of one scanline in pData to the start of the next. If defaulted the size of the datatype eBufType * nBufXSize is used.
*/
G_imageReadMutex.lock();
m_pBand_GDAL_->RasterIO(GF_Read, 0, 0,
m_nXSize_GDAL_, m_nYSize_GDAL_,
m_pDataBufStart_GDAL_,
m_nBufXSize_GDAL_, m_nBufYSize_GDAL_,
m_eBufDataType_GDAL_,
m_nBufPixSpace_GDAL_,
m_nBufLineSpace_GDAL_);
//G_imageReadComplete.wakeAll();
G_imageReadMutex.unlock();
}
//读图函数：使用GDAL和OpenCV读取图像文件
void RSI_ImageReaderThread::ReadFile(const QString& imageFileName)
{
m_imageFileName = imageFileName;
if (m_imageFileName.isEmpty() || m_pCanvas == NULL)
return;
//选择打开图像的驱动方式: GDAL or OpenCV
const QString strOpenCV = tr("OpenCV");
const QString strGDAL = tr("GDAL");
IMG_OPEN_DRIVER_T eOpenDriverType;
QStringList items; //for input dialog
items height)
.arg(pIplImg->depth)
.arg(pIplImg->nChannels);
rsiFileInfo += tempStr;
//显示
uchar* imageDataBuffer = NULL;
QImage* pQImage = IplImageToQImage(pIplImg, &imageDataBuffer); //转换成QImage
//创建图层
RSI_ImageLayer* pImageLayer = new RSI_ImageLayer(*pQImage);
pImageLayer->setLayerName(m_imageFileName);
pImageLayer->setImageLayerInfo(rsiFileInfo);
pImageLayer->setBufferToBeFree(imageDataBuffer); //记录到m_bufferToBeFree中，待关闭图像时释放buffer内存
pImageLayer->setVisible(false); //默认不显示
pImageLayer->setFlag(QGraphicsItem::ItemIsMovable, false);
m_pCanvas->AddImageLayer(pImageLayer);
//释放临时内存
cvReleaseImage(&pIplImg);
delete pQImage;
goto _READ_IMAGE_FILE_OVER; //读取完毕直接返回
}
else if (eOpenDriverType == RSI_USE_GDAL)
{
//如果是遥感图像格式，则使用GDAL打开
GDALAllRegister();
GDALDataset* poDataset; //GDAL数据集
poDataset = (GDALDataset*)GDALOpen(m_imageFileName.toStdString().c_str(), GA_ReadOnly);
if( poDataset == NULL )
{
//QMessageBox::warning(NULL, tr("Failed"), tr("GDAL Open Image File Failed!"), tr("OK"));
return;
}
//读取遥感图像信息
QString rsiFileInfo; //遥感图像信息
QString tempStr;
tempStr = tr("GDAL Driver: %1/%2/n")
.arg(QString(poDataset->GetDriver()->GetDescription()))
.arg(QString(poDataset->GetDriver()->GetMetadataItem(GDAL_DMD_LONGNAME)));
rsiFileInfo += tempStr;
tempStr = tr("Size is %1 x %2 x %3(bands)/n")
.arg(poDataset->GetRasterXSize())
.arg(poDataset->GetRasterYSize())
.arg(poDataset->GetRasterCount());
rsiFileInfo += tempStr;
/*
地理参考坐标系统:
The returned string defines the projection coordinate system of the image in OpenGIS WKT format.
*/
double adfGeoTransform[6];
QString prjRefStr(poDataset->GetProjectionRef());
if (prjRefStr.isEmpty())
{
prjRefStr = tr("none"); //如果此字段为空，那么投影方式可能是GCPs
}
else //非空时才读取仿射变换信息
{
tempStr = tr("Projection(WKT format) is '%1'/n").arg(prjRefStr);
rsiFileInfo += tempStr;
//非空时才读取仿射变换信息
/*
仿射地理变换信息：
adfGeoTransform[0] // top left x
adfGeoTransform[1] // w-e pixel resolution
adfGeoTransform[2] // rotation, 0 if image is "north up"
adfGeoTransform[3] // top left y
adfGeoTransform[4] // rotation, 0 if image is "north up"
adfGeoTransform[5] // n-s pixel resolution
计算某一点的地理坐标可以如下计算：
Xgeo = GT(0) + Xpixel*GT(1) + Yline*GT(2)
Ygeo = GT(3) + Xpixel*GT(4) + Yline*GT(5)
*/
poDataset->GetGeoTransform( adfGeoTransform );
tempStr = tr("Origin Cords = (%1,%2)/n")
.arg(adfGeoTransform[0], 0, 'f', 6)
.arg(adfGeoTransform[3], 0, 'f', 6);
rsiFileInfo += tempStr;
tempStr = tr("Pixel Size(degrees) = (%1,%2)/n")
.arg(adfGeoTransform[1], 0, 'f', 6)
.arg(adfGeoTransform[5], 0, 'f', 6);
rsiFileInfo += tempStr;
}
//读取GCPs投影方式
QString prjGCPStr(poDataset->GetGCPProjection());
if (prjGCPStr.isEmpty())
{
prjGCPStr = tr("none");
}
else
{
tempStr = tr("GCPs Prj is '%1'/n").arg(prjGCPStr);
rsiFileInfo += tempStr;
int nGCPs = poDataset->GetGCPCount(); //获得GCP控制点的个数
const GDAL_GCP* pGCPs = poDataset->GetGCPs(); //获得GCP控制点
tempStr = tr("Number of GCPs: %1/n").arg(nGCPs);
rsiFileInfo += tempStr;
//由GCPs获得仿射变换参数
GDALGCPsToGeoTransform( nGCPs, pGCPs, adfGeoTransform, TRUE );
tempStr = tr("Origin Cords = (%1,%2)/n")
.arg(adfGeoTransform[0], 0, 'f', 6)
.arg(adfGeoTransform[3], 0, 'f', 6);
rsiFileInfo += tempStr;
tempStr = tr("Pixel Size(degrees) = (%1,%2)/n")
.arg(adfGeoTransform[1], 0, 'f', 6)
.arg(adfGeoTransform[5], 0, 'f', 6);
rsiFileInfo += tempStr;
}
//根据adfGeoTransform计算每个像素代表多少米距离
double oriCordL = adfGeoTransform[0];
double oriCordB = adfGeoTransform[3];
double tempL, tempB;
GDALApplyGeoTransform(adfGeoTransform, 100, 0, &tempL, &tempB); //计算(100,0)像素处的L,B坐标
double tempDist = CalcDistanceByLB(oriCordL, oriCordB, tempL, tempB);
double xMeterPerPix = tempDist * 1000.0 / 100.0; //单位为米
GDALApplyGeoTransform(adfGeoTransform, 0, 100, &tempL, &tempB); //计算(0,100)像素处的L,B坐标
tempDist = CalcDistanceByLB(oriCordL, oriCordB, tempL, tempB);
double yMeterPerPix = tempDist * 1000.0 / 100.0; //单位为米
tempStr = tr("Pixel Size(meters) = (%1,%2)/n")
.arg(xMeterPerPix, 0, 'f', 6)
.arg(yMeterPerPix, 0, 'f', 6);
rsiFileInfo += tempStr;
//读取图像大小
int nBandCount = poDataset->GetRasterCount();
int nImgSizeX = poDataset->GetRasterXSize();
int nImgSizeY = poDataset->GetRasterYSize();
//读取图像数据点类型
GDALRasterBand* preBand = poDataset->GetRasterBand(1); //预读取遥感的第一个波段
GDALDataType prePixType = preBand->GetRasterDataType();
int prePixSize = GDALGetDataTypeSize(prePixType) / 8; //GDALGetDataTypeSize得到的是bit
if (nBandCount >= 3 && prePixSize == 1) //多于3个波段的图像给用户选择是否合并前三个波段用RGB显示
{
//QMessageBox::StandardButton userAnswer = QMessageBox::Yes;
QMessageBox::StandardButton userAnswer = QMessageBox::question(m_pCanvas,
tr("Choose image show format"),
tr("%1 bands detected!/nWill you show it in RGB color format using the first three bands?").arg(nBandCount),
QMessageBox::Yes|QMessageBox::No);
if (userAnswer == QMessageBox::Yes)
{
//--------将前三个波段组合成RGB颜色显示, 使用QImage::Format_RGB32格式显示---------
uchar* pafScanblock = (uchar*)malloc(4*(nImgSizeX)*(nImgSizeY)); //Format_RGB32每个像素4个字节: RGBA
//assert(pafScanblock != NULL);
if (pafScanblock == NULL)
{
QMessageBox::warning(m_pCanvas, tr("Failed"), tr("Open Image File Failed! Out of memory!"), tr("OK"));
delete poDataset; //关闭数据集
goto _READ_IMAGE_FILE_OVER;
}
memset(pafScanblock, 0x00, 4*(nImgSizeX)*(nImgSizeY));
//读取图像
for (int nB=0; nBGetRasterBand(nB+1); //遥感的一个波段
assert(poBand != NULL);
GDALDataType pixType = poBand->GetRasterDataType();
int pixSize = GDALGetDataTypeSize(pixType) / 8; //GDALGetDataTypeSize得到的是bit
assert(pixSize == 1); //彩色图的每个通道的size应该是1个字节
//用RasterIO()读取图像数据
//poBand->RasterIO(GF_Read, 0, 0, nImgSizeX, nImgSizeY,
// pafScanblock+nB, nImgSizeX, nImgSizeY, pixType, 4, 0);
G_imageReadMutex.lock();
m_pBand_GDAL_ = poBand;
m_nXSize_GDAL_ = nImgSizeX;
m_nYSize_GDAL_ = nImgSizeY;
m_pDataBufStart_GDAL_ = pafScanblock+nB; //注意！
m_nBufXSize_GDAL_ = nImgSizeX;
m_nBufYSize_GDAL_ = nImgSizeY;
m_eBufDataType_GDAL_ = pixType;
m_nBufPixSpace_GDAL_ = 4; //注意！
m_nBufLineSpace_GDAL_ = 0;
start(); //开始新线程
G_imageReadMutex.unlock();
while (this->isRunning())
{
progressBarValue += 8;
progress.setValue(progressBarValue % 100);
QThread::msleep(1000); //注意：这里是让主线程sleep，而不是新线程
}
//G_imageReadComplete.wait(&G_imageReadMutex);
//progressBarValue += 20;
//progress.setValue(progressBarValue % 100);
//G_imageReadMutex.unlock();
}
//组合成彩色QImage
QImage* pQImage = NULL;
pQImage = new QImage(pafScanblock, nImgSizeX, nImgSizeY, QImage::Format_RGB32); //注意：QImage: each scanline of data in the image must also be 32-bit aligned
//assert(pQImage != NULL);
if (pQImage == NULL)
{
QMessageBox::warning(m_pCanvas, tr("Failed"), tr("Open Image File Failed! Out of memory!"), tr("OK"));
free(pafScanblock);
delete poDataset; //关闭数据集
goto _READ_IMAGE_FILE_OVER;
}
//创建图层
RSI_ImageLayer* pImageLayer = new RSI_ImageLayer(*pQImage);
if (pImageLayer == NULL)
{
QMessageBox::warning(m_pCanvas, tr("Failed"), tr("Open Image File Failed! Out of memory!"), tr("OK"));
free(pafScanblock);
delete pQImage;
delete poDataset; //关闭数据集
goto _READ_IMAGE_FILE_OVER;
}
pImageLayer->setLayerName(tr("3 bands composited color image") + QString("(") + m_imageFileName + QString(")"));
pImageLayer->setImageLayerInfo(rsiFileInfo);
pImageLayer->setGeoTransform(adfGeoTransform);
pImageLayer->setResolutionX(xMeterPerPix);
pImageLayer->setResolutionY(yMeterPerPix);
pImageLayer->setBufferToBeFree(pafScanblock); //记录到m_bufferToBeFree中，待关闭图像时释放buffer内存
pImageLayer->setVisible(false); //默认不显示
pImageLayer->setFlag(QGraphicsItem::ItemIsMovable, false);
m_pCanvas->AddImageLayer(pImageLayer);
//释放临时对象
delete pQImage;
delete poDataset; //关闭数据集
goto _READ_IMAGE_FILE_OVER; //读取完毕直接返回
}
}
//------------按照波段分别读取，每个波段显示成一个灰度图-------------
//读取图像
for (int nB=0; nBGetRasterBand(nB+1); //遥感的一个波段
assert(poBand != NULL);
GDALDataType pixType = poBand->GetRasterDataType();
int pixSize = GDALGetDataTypeSize(pixType) / 8; //GDALGetDataTypeSize得到的是bit
//注意：QImage: each scanline of data in the image must also be 32-bit aligned
int nImgSizeXAligned32bit = Get32bitAlignedWidthBytes(nImgSizeX, pixSize) / pixSize; //对齐后的像素宽度
uchar* pafScanblock = (uchar*)malloc(pixSize*(nImgSizeXAligned32bit)*(nImgSizeY));
//assert(pafScanblock != NULL);
if (pafScanblock == NULL)
{
QMessageBox::warning(m_pCanvas, tr("Failed"), tr("Open Image File Failed! Out of memory!"), tr("OK"));
delete poDataset; //关闭数据集
goto _READ_IMAGE_FILE_OVER;
}
//用RasterIO()读取图像数据
//poBand->RasterIO(GF_Read, 0, 0, nImgSizeX, nImgSizeY,
// pafScanblock, nImgSizeXAligned32bit, nImgSizeY, pixType, 0, 0);
G_imageReadMutex.lock();
m_pBand_GDAL_ = poBand;
m_nXSize_GDAL_ = nImgSizeX;
m_nYSize_GDAL_ = nImgSizeY;
m_pDataBufStart_GDAL_ = pafScanblock;
m_nBufXSize_GDAL_ = nImgSizeXAligned32bit;
m_nBufYSize_GDAL_ = nImgSizeY;
m_eBufDataType_GDAL_ = pixType;
m_nBufPixSpace_GDAL_ = 0;
m_nBufLineSpace_GDAL_ = 0;
start(); //开始新线程
G_imageReadMutex.unlock();
while (this->isRunning())
{
progressBarValue += 8;
progress.setValue(progressBarValue % 100);
QThread::msleep(1000); //注意：这里是让主线程sleep，而不是新线程
}
//G_imageReadComplete.wait(&G_imageReadMutex);
//progressBarValue += 20;
//progress.setValue(progressBarValue % 100);
//G_imageReadMutex.unlock();
//转化为QImage
QImage* pQImage = CvtGDALBufferToGrayQImage(pafScanblock, nImgSizeX, nImgSizeY, pixType);
if (pQImage == NULL)
{
QMessageBox::warning(m_pCanvas, tr("Failed"), tr("Open Image File Failed! Out of memory!"), tr("OK"));
free(pafScanblock);
delete poDataset; //关闭数据集
goto _READ_IMAGE_FILE_OVER;
}
//创建图层
RSI_ImageLayer* pImageLayer = new RSI_ImageLayer(*pQImage);
if (pImageLayer == NULL)
{
QMessageBox::warning(m_pCanvas, tr("Failed"), tr("Open Image File Failed! Out of memory!"), tr("OK"));
free(pafScanblock);
delete pQImage;
delete poDataset; //关闭数据集
goto _READ_IMAGE_FILE_OVER;
}
pImageLayer->setLayerName(tr("band") + QString("%1").arg(nB+1) + tr(": pixel type = ") + QString(GDALGetDataTypeName(pixType)) + QString("(") + QString("%1").arg(pixSize) + tr(" bytes)"));
pImageLayer->setImageLayerInfo(rsiFileInfo);
pImageLayer->setGeoTransform(adfGeoTransform);
pImageLayer->setResolutionX(xMeterPerPix);
pImageLayer->setResolutionY(yMeterPerPix);
pImageLayer->setBufferToBeFree(pafScanblock); //记录到m_bufferToBeFree中，待关闭图像时释放buffer内存
pImageLayer->setVisible(false); //默认不显示
pImageLayer->setFlag(QGraphicsItem::ItemIsMovable, false);
m_pCanvas->AddImageLayer(pImageLayer);
//释放临时对象
delete pQImage;
}
delete poDataset; //关闭数据集
goto _READ_IMAGE_FILE_OVER; //读取完毕直接返回
} //end if eOpenDriverType
_READ_IMAGE_FILE_OVER:
return;
}
//将图像宽度字节数32位对齐，返回对齐后的字节数
int RSI_ImageReaderThread::Get32bitAlignedWidthBytes(int nImgPixWidth, int nBytesPerPix)
{
return ((nImgPixWidth*nBytesPerPix/4) + (((nImgPixWidth*nBytesPerPix)%4)==0 ? 0 : 1)) * 4;
}
//GDAL buffer 转 QImage
QImage* RSI_ImageReaderThread::CvtGDALBufferToGrayQImage(uchar* &pafScanblock, int nImgSizeX, int nImgSizeY, GDALDataType pixType)
{
//准备颜色表
QVector vcolorTable;
for (int i = 0; i setColorTable(vcolorTable);
}
return pQImage;
}
//将非8位的灰度图转换为8位灰度图
template
void RSI_ImageReaderThread::CvtImgDataTo8bitGrayImg(T* pOrigData, int nImgSizeX, int nImgSizeXAligned32bit, int nImgSizeY, uchar* pNew8bitImgBuf, int nNewImgSizeXAligned32bit, T maxVal, T minVal)
{
int oldIdxPos = 0;
int newIdxPos = 0;
T minPixVal = maxVal;
T maxPixVal = minVal;
T curPixVal;
//找到最大值和最小值
for (int niY=0; niY 文章版权声明：除非注明，否则均为VPS857原创文章，转载或复制请以超链接形式并注明出处。
还没有评论，来说两句吧...