波段运算支持四则运算（加减乘除）

Qt、GDAL遥感影像显示

1、创建类

class GisBandMathProcess
	: public GisRasterProcess
{
	OGIS_DECLARE_MEMBERS(GisBandMathProcess);
public:

	GisBandMathProcess();

	~GisBandMathProcess();

	virtual QString processName() const;

	virtual OGis::ErrorStatus beginProcess();

	virtual OGis::ErrorStatus processData(ProcessData *pData) const;

	virtual OGis::ErrorStatus process(MatArray &desMat, const MatArray &srcMat, const MatArray &maskArray = MatArray());

private:

	void parseExperssion();

private:

	uint m_datasetIndex;
	QString m_exp, m_sPath;
	UIntArray m_bandIndexs;
	QStringList m_expStrings;

};

2、实现类

#include "StdAfx.h"

OGIS_CONS_DEFINE_MEMBERS_WORKER(GisBandMathProcess, GisRasterProcess)

GisBandMathProcess::GisBandMathProcess()
	: GisRasterProcess()
{

}

GisBandMathProcess::~GisBandMathProcess()
{

}

void GisBandMathProcess::parseExperssion()
{
	QString expStr(m_exp);
	QRegularExpression regExp("B|b[0-9]");
	QRegularExpressionMatch match;
	do {
		match = regExp.match(expStr);
		if(match.hasMatch()) {
			QString band = match.captured();
			QString tagName = QS("%1").arg(band);
			m_expStrings.append(tagName);
			expStr.remove(band, Qt::CaseInsensitive);
		}
		else {
			break;
		}
	} while(true);
}

QString GisBandMathProcess::processName() const
{
	return QS("波段运算");
}

OGis::ErrorStatus GisBandMathProcess::beginProcess()
{
	GisBandMathDialog dialog;
	if(QDialog::Rejected == dialog.exec()) {
		return OGis::eUserCancel;
	}
	m_exp = dialog.bandMathExp().trimmed().toLower();
	if(m_exp.isEmpty()) {
		return OGis::eInvalidString;
	}
	parseExperssion();
	GisVariableToBandsPairings varDialog;
	varDialog.setExpression(m_exp);
	varDialog.adjustSize();
	if(QDialog::Rejected == varDialog.exec()) {
		return OGis::eUserCancel;
	}
	m_sPath = varDialog.getSaveFileName();
	varDialog.getBandsPairing(m_datasetIndex, m_bandIndexs);
	return GisRasterProcess::beginProcess();
}

double GetMatValue(const cv::Mat &src, int row, int col)
{
	int type = src.type();
	switch(type) {
		case CV_8U: return (double)src.at<byte>(row, col);
		case CV_8S: return (double)src.at<char>(row, col);
		case CV_16U: return (double)src.at<ushort>(row, col);
		case CV_16S: return (double)src.at<short>(row, col);
		case CV_32S: return (double)src.at<int>(row, col);
		case CV_32F: return (double)src.at<float>(row, col);
		case CV_64F: return (double)src.at<double>(row, col);
		default:
			break;
	}
	return 0;
}

OGis::ErrorStatus GisBandMathProcess::processData(ProcessData *pData) const
{
	if(nullptr == pData) {
		return OGis::eNullObjectPointer;
	}
	RasterData *mpRasterData = (RasterData *)pData;
	mpRasterData->bDisplay = true;
	mpRasterData->datasetIndex = m_datasetIndex;
	mpRasterData->exportDataType = 7;
	mpRasterData->oFormat = OGis::kEnvi;
	mpRasterData->blockSize = {1024, 1024};
	mpRasterData->itemArray = m_bandIndexs;
	mpRasterData->outputFileName = m_sPath;
	mpRasterData->resultSize = {0, 0, 1};
	pData = mpRasterData;
	return OGis::eOk;
}

OGis::ErrorStatus GisBandMathProcess::process(MatArray &desMat, const MatArray &srcMat, const MatArray &maskArray)
{
	QJSEngine myEngine; int rows = srcMat[0].rows, cols = srcMat[0].cols;
	GisProgressBarPointer progressBar(QS("[%1]").arg(m_exp), rows, true);
	cv::Mat des = cv::Mat(rows, cols, CV_64FC1);
	for(int row = 0; row < rows; ++row) {
		progressBar.setProgressPos(row);
		for(int col = 0; col < cols; ++col) {
			QString exp = m_exp; bool bCal = true;
			for(int idx = 0; idx < srcMat.size(); ++idx) {
				cv::Mat mat = srcMat[idx]; double val = 0;
				if(maskArray.size() > 0) {
					cv::Mat mask = maskArray[idx];
					val = GetMatValue(mask, row, col);
					bCal = bCal && (val != 0);
				}
				val = GetMatValue(mat, row, col);
				exp.replace(m_expStrings[idx], QS("%1").arg(val));
			}
			if(bCal) {
				double val = myEngine.evaluate(exp).toNumber();
				des.at<double>(row, col) = val;
			}
			else {
				des.at<double>(row, col) = nan("");
			}
		}
	}
	desMat << des;
	return OGis::eOk;
}

GisBandMathDialog：波段运算表达式输出对话框；

GisVariableToBandsPairings：定义波段运算中（b1,b2,b3）变量对应的波段；

在 processData(ProcessData *pData) 接口中具体参数如下：

mpRasterData->bDisplay = true; //处理结果自动打开
mpRasterData->datasetIndex = m_datasetIndex; //数据集处理的索引
mpRasterData->exportDataType = 7; //输出数据的类型
mpRasterData->oFormat = OGis::kEnvi; //结果数据保存为ENVI格式
mpRasterData->blockSize = {1024, 1024}; //数据处理分块大小
mpRasterData->itemArray = m_bandIndexs; //数据处理的波段
mpRasterData->outputFileName = m_sPath; //结果数据保存的全路径
mpRasterData->resultSize = {0, 0, 1}; //结果数据的大小（0:为原始数据大小）

接口 process(MatArray &desMat, const MatArray &srcMat, const MatArray &maskArray) 是数据处理的核心：

利用Qt提供的QJSEngine类（支持四则运算），所以我们将波段表达式中的b1，b2，b3等参数进行像素值替换。例如：计算NDVI = (b1-b2)/(b1+b2)。b1是近红外波段，b2是红波段，我们将b1和b2替换成对应波段的像素值，即表达式为：(24-62)/(24+62)，然后通过double val = myEngine.evaluate(exp).toNumber(); 计算结果为：-0.4418，这个值也就是对应输出结果影像的像素值。