diff --git a/NiirsLinux.sln b/NiirsLinux.sln
new file mode 100644
index 0000000..9e865ae
--- /dev/null
+++ b/NiirsLinux.sln
@@ -0,0 +1,22 @@
+
+Microsoft Visual Studio Solution File, Format Version 12.00
+# Visual Studio Version 16
+VisualStudioVersion = 16.0.33027.164
+MinimumVisualStudioVersion = 10.0.40219.1
+Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "NiirsLinux", "NiirsLinux\NiirsLinux.vcxproj", "{22B7E044-B017-466B-985E-D3FDF6435C2E}"
+EndProject
+Global
+ GlobalSection(SolutionConfigurationPlatforms) = preSolution
+ Release|x64 = Release|x64
+ EndGlobalSection
+ GlobalSection(ProjectConfigurationPlatforms) = postSolution
+ {22B7E044-B017-466B-985E-D3FDF6435C2E}.Release|x64.ActiveCfg = Release|x64
+ {22B7E044-B017-466B-985E-D3FDF6435C2E}.Release|x64.Build.0 = Release|x64
+ EndGlobalSection
+ GlobalSection(SolutionProperties) = preSolution
+ HideSolutionNode = FALSE
+ EndGlobalSection
+ GlobalSection(ExtensibilityGlobals) = postSolution
+ SolutionGuid = {4D31A63B-64E0-4A37-A0E6-6AD47E2C0284}
+ EndGlobalSection
+EndGlobal
diff --git a/NiirsLinux/NiirsLinux.vcxproj b/NiirsLinux/NiirsLinux.vcxproj
new file mode 100644
index 0000000..02c57ea
--- /dev/null
+++ b/NiirsLinux/NiirsLinux.vcxproj
@@ -0,0 +1,85 @@
+
+
+
+
+ Release
+ x64
+
+
+
+ {22B7E044-B017-466B-985E-D3FDF6435C2E}
+ QtVS_v304
+ 10.0.22000.0
+ $(MSBuildProjectDirectory)\QtMsBuild
+
+
+
+ Application
+ v142
+
+
+
+
+
+
+ Qt5.11.2
+ core;xml;gui;widgets;printsupport
+ release
+
+
+
+
+
+
+
+
+
+
+
+
+ D:\qgis\osgeo4w\include;D:\qgis\osgeo4w\apps\Qt5\include;D:\OpenCV\opencv455\include;D:\OpenCV\opencv455\include\opencv2;D:\qgis\osgeo4w\apps\Qt5\include\QtXml;$(IncludePath)
+ D:\qgis\osgeo4w\lib;D:\qgis\osgeo4w\apps\Qt5\lib;D:\OpenCV\opencv455\x64\vc16\lib;$(LibraryPath)
+
+
+
+ opencv_core455.lib;opencv_imgcodecs455.lib;opencv_imgproc455.lib;gdal_i.lib;%(AdditionalDependencies)
+
+
+ true
+
+
+
+
+ true
+ true
+ EditAndContinue
+ Disabled
+
+
+ Windows
+ true
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
\ No newline at end of file
diff --git a/NiirsLinux/NiirsLinux.vcxproj.filters b/NiirsLinux/NiirsLinux.vcxproj.filters
new file mode 100644
index 0000000..2dabaa6
--- /dev/null
+++ b/NiirsLinux/NiirsLinux.vcxproj.filters
@@ -0,0 +1,53 @@
+
+
+
+
+ {4FC737F1-C7A5-4376-A066-2A32D752A2FF}
+ qml;cpp;c;cc;cxx;def;odl;idl;hpj;bat;asm;asmx
+
+
+ {93995380-89BD-4b04-88EB-625FBE52EBFB}
+ h;hh;hpp;hxx;hm;inl;inc;xsd
+
+
+ {67DA6AB6-F800-4c08-8B7A-83BB121AAD01}
+ qrc;rc;ico;cur;bmp;dlg;rc2;rct;bin;rgs;gif;jpg;jpeg;jpe;resx;tiff;tif;png;wav;mfcribbon-ms
+
+
+ {99349809-55BA-4b9d-BF79-8FDBB0286EB3}
+ ui
+
+
+ {639EADAA-A684-42e4-A9AD-28FC9BCB8F7C}
+ ts
+
+
+
+
+ Source Files
+
+
+ Source Files
+
+
+
+
+ Form Files
+
+
+
+
+ Resource Files
+
+
+
+
+ Header Files
+
+
+
+
+ Form Files
+
+
+
\ No newline at end of file
diff --git a/NiirsLinux/main.cpp b/NiirsLinux/main.cpp
new file mode 100644
index 0000000..d7fc870
--- /dev/null
+++ b/NiirsLinux/main.cpp
@@ -0,0 +1,28 @@
+#include "niirs.h"
+#include
+#include
+#include
+
+int main(int argc, char *argv[])
+{
+ QString xmlInterface = "";
+
+ for (int i = 0; i < argc; i++)
+ {
+ QString temp = QString::fromLocal8Bit(argv[i]);
+ if (temp == "--xml" && i + 1 < argc)
+ {
+ xmlInterface = argv[i + 1];
+ break;
+ }
+ }
+
+ QApplication a(argc, argv);
+ QStyle* style = QStyleFactory::create("Fusion");
+ a.setStyle(style);
+
+ NIIRS w;
+ w.setXmlInterface(xmlInterface);
+ w.show();
+ return a.exec();
+}
diff --git a/NiirsLinux/niirs.cpp b/NiirsLinux/niirs.cpp
new file mode 100644
index 0000000..a84ae4d
--- /dev/null
+++ b/NiirsLinux/niirs.cpp
@@ -0,0 +1,3229 @@
+#include "niirs.h"
+
+NIIRS::NIIRS(QWidget* parent)
+ : QMainWindow(parent)
+{
+ ui.setupUi(this);
+ connect(ui.radioVisible, &QRadioButton::clicked, this, &NIIRS::onRadioBtn_Visi);
+ connect(ui.radioMidInfrared, &QRadioButton::clicked, this, &NIIRS::onRadioBtn_Infra);
+ connect(ui.radioLonInfrared, &QRadioButton::clicked, this, &NIIRS::onRadioBtn_Infra);
+
+ ui.textFormula->setText(QString::fromLocal8Bit("可见光数据计算公式:\n") +
+ "NIIRS = 10.251 - a*lgGSD + b*lgRER - 0.656*H - 0.334*(G/SNR) | RER>=0.9:a=3.32,b=1.559;RER<0.9:a=3.16,b=2.817");
+
+ connect(ui.pbtCalculate, &QPushButton::clicked, this, &NIIRS::calculateNiirsScore);
+ connect(ui.pbtObjectivePara, &QPushButton::clicked, this, &NIIRS::calculateObjectivePara);
+ //connect(ui.pbtObjectivePara, &QPushButton::clicked, this, &NIIRS::calculateEvaluationImage);
+ //connect(ui.pbtRER, &QPushButton::clicked, this, &NIIRS::calculateRER);
+
+
+ ui.lineEdit_GSDX->setText("1.0");
+ ui.lineEdit_GSDY->setText("1.0");
+ //ui.lineEdit_RERX->setText("0.986");
+ //ui.lineEdit_RERY->setText("0.965");
+ //ui.lineEdit_MTFC_H->setText("0.9");
+
+
+ ui.lineEdit_RER->setText("0.986");
+ ui.lineEdit_H->setText("0.9");
+
+ ui.lineEdit_MTFC_G->setText("1");
+ ui.lineEdit_SNR->setText("20");
+ ui.lineEdit_PIQE->setText("0.986");
+
+ QString inputImgPath = mXmlImgPath;
+ GDALAllRegister();
+ CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");
+ string tiffEdgePath = inputImgPath.toStdString();
+ GDALDataset* poDataset;
+ poDataset = (GDALDataset*)GDALOpen(tiffEdgePath.c_str(), GA_ReadOnly);
+ if (NULL != poDataset)
+ {
+ for (int i = 0; i < poDataset->GetRasterCount(); i++)
+ ui.comboBoxBand->addItem(QString("band:%1").arg(i + 1));
+ }
+ GDALClose(poDataset);
+}
+
+QStringList NIIRS::getAllFiles(QString path, QString fileType)
+{
+ // 判断路径是否存在
+ QDir dir(path);
+ if (!dir.exists())
+ {
+ QMessageBox::warning(this, "warning", "no such path:\n" + path);
+ return QStringList();
+ }
+
+ dir.setFilter(QDir::Files | QDir::NoSymLinks);
+ QFileInfoList infoList = dir.entryInfoList();
+
+ int fileCount = infoList.count();
+ if (fileCount <= 0)
+ {
+ QMessageBox::warning(this, "warning", "No files in folder:\n" + path);
+ return QStringList();
+ }
+
+ QStringList files;
+ for (int i = 0; i < fileCount; i++)
+ {
+ QFileInfo fileInfo = infoList.at(i);
+ QString suffix = fileInfo.suffix();
+ if (QString::compare(suffix, fileType, Qt::CaseInsensitive) == 0)
+ {
+ QString absoluteFilePath = fileInfo.absoluteFilePath();
+ files.append(absoluteFilePath);
+ }
+ }
+ return files;
+}
+
+void NIIRS::onRadioBtn_Visi()
+{
+ ui.textFormula->setText(QString::fromLocal8Bit("可见光数据计算公式:\n") +
+ "NIIRS = 10.251 - a*lgGSD + b*lgRER - 0.656*H - 0.334*(G/SNR) | RER>=0.9:a=3.32,b=1.559;RER<0.9:a=3.16,b=2.817");
+}
+void NIIRS::onRadioBtn_Infra()
+{
+ if (ui.radioMidInfrared->isChecked())
+ ui.textFormula->setText(QString::fromLocal8Bit("中红外数据计算公式:\n")
+ + "NIIRS = 10.751 - a*lgGSD + b*lgRER - 0.656*H - 0.334*(G/SNR) | RER>=0.9:a=3.32,b=1.559;RER<0.9:a=3.16,b=2.817");
+ else if (ui.radioLonInfrared->isChecked())
+ ui.textFormula->setText(QString::fromLocal8Bit("长波红外数据计算公式:\n")
+ + "NIIRS = 10.751 - a*lgGSD + b*lgRER - 0.656*H - 0.334*(G/SNR) | RER>=0.9:a=3.32,b=1.559;RER<0.9:a=3.16,b=2.817");
+}
+
+void NIIRS::on_pbtOpenImg_clicked()
+{
+ QString openImg = QFileDialog::getOpenFileName(this, "", "", "*.tif *.tiff");
+ if (!openImg.isEmpty())
+ ui.lineEdit_img->setText(openImg);
+ else
+ ui.lineEdit_INFO->setText(QString::fromLocal8Bit("未打开数据"));
+ //QString openImg = QFileDialog::getExistingDirectory(this, "", "");
+ //if (!openImg.isEmpty())
+ // ui.lineEdit_img->setText(openImg);
+ //else
+ // ui.lineEdit_INFO->setText(QString::fromLocal8Bit("未打开路径"));
+}
+
+//void NIIRS::on_pbtOpenEdgeX_clicked()
+//{
+// QString edgeX = QFileDialog::getOpenFileName(this, "", "", "*.tif *.tiff");
+// if (!edgeX.isEmpty())
+// ui.lineEdit_EdgeX->setText(edgeX);
+// else
+// ui.lineEdit_INFO->setText(QString::fromLocal8Bit("未打开数据"));
+// //QString edgeX= QFileDialog::getOpenFileName(this, "", "", "*.jpg");
+// //if (!edgeX.isEmpty())
+// // ui.lineEdit_EdgeX->setText(edgeX);
+// //else
+// // ui.lineEdit_INFO->setText(QString::fromLocal8Bit("未打开数据"));
+//}
+
+//void NIIRS::on_pbtOpenEdgeY_clicked()
+//{
+// QString edgeY = QFileDialog::getOpenFileName(this, "", "", "*.tif *.tiff");
+// if (!edgeY.isEmpty())
+// ui.lineEdit_EdgeY->setText(edgeY);
+// else
+// ui.lineEdit_INFO->setText(QString::fromLocal8Bit("未打开数据"));
+// //QString edgeY = QFileDialog::getOpenFileName(this, "", "", "*.jpg");
+// //if (!edgeY.isEmpty())
+// // ui.lineEdit_EdgeY->setText(edgeY);
+// //else
+// // ui.lineEdit_INFO->setText(QString::fromLocal8Bit("未打开数据"));
+//}
+
+int NIIRS::getImageInfo(const char* srcImgPath, struct ImgInfo* imgInfo)
+{
+ GDALAllRegister();
+ CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO"); //中文路径
+ GDALDataset* hSrcDS;
+ hSrcDS = (GDALDataset*)GDALOpen(srcImgPath, GA_ReadOnly);
+
+ if (NULL == hSrcDS)
+ {
+ GDALClose(hSrcDS);
+ return -1;
+ }
+
+ GDALRasterBand* rasterBand;
+ rasterBand = hSrcDS->GetRasterBand(1);
+
+ imgInfo->dataType = rasterBand->GetRasterDataType();
+ imgInfo->bandCount = hSrcDS->GetRasterCount();
+ imgInfo->xSize = hSrcDS->GetRasterXSize();
+ imgInfo->ySize = hSrcDS->GetRasterYSize();
+ imgInfo->xOffset = 0;
+ imgInfo->yOffset = 0;
+
+ strcpy(imgInfo->projWkt, hSrcDS->GetProjectionRef());
+ hSrcDS->GetGeoTransform(imgInfo->geoTransform);
+ imgInfo->noData = rasterBand->GetNoDataValue();
+
+ GDALClose(hSrcDS);
+ return 0;
+}
+
+int NIIRS::readTiffFile(const char* srcImgPath, ImgInfo imgInfo, void* imgArray, int bandIndex)
+{
+ GDALAllRegister();
+ CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO"); //中文路径
+ GDALDataset* poSrcDS = (GDALDataset*)GDALOpen(srcImgPath, GA_ReadOnly);
+
+ //返回错误信息
+ int errInfo = 0;
+
+ if (NULL == poSrcDS)
+ {
+ GDALClose(poSrcDS);
+ return -1;
+ }
+
+ //输入波段索引大于数据波段数
+ if (bandIndex > imgInfo.bandCount)
+ {
+ GDALClose(poSrcDS);
+ return -1;
+ }
+ //读取全部波段
+ if (0 == bandIndex)
+ {
+ errInfo = poSrcDS->RasterIO(GF_Read, imgInfo.xOffset, imgInfo.yOffset, imgInfo.xSize, imgInfo.ySize,
+ imgArray, imgInfo.xSize, imgInfo.ySize, imgInfo.dataType, imgInfo.bandCount, 0, 0, 0, 0);
+ }
+ //读取指定波段
+ else
+ {
+ GDALRasterBand* poBand = poSrcDS->GetRasterBand(bandIndex);
+ errInfo = poBand->RasterIO(GF_Read, imgInfo.xOffset, imgInfo.yOffset, imgInfo.xSize, imgInfo.ySize,
+ imgArray, imgInfo.xSize, imgInfo.ySize, imgInfo.dataType, 0, 0);
+ }
+ GDALClose((GDALDatasetH)poSrcDS);
+ return errInfo;
+}
+
+void NIIRS::initTableWidget(int rowNum, int colNum)
+{
+ ui.tableWidget->clear();
+ ui.tableWidget->setColumnCount(colNum);
+ ui.tableWidget->setRowCount(rowNum);
+ QStringList header;
+ for (int i = 0; i < colNum; i++)
+ header << "band " + QString::number(i + 1);
+ //header << "band 1" << "band 2" << "band 3" << "band 4";
+ ui.tableWidget->setHorizontalHeaderLabels(header);
+ header.clear();
+ header << QString::fromLocal8Bit("亮度") << QString::fromLocal8Bit("标准差") << QString::fromLocal8Bit("对比度") << QString::fromLocal8Bit("熵")
+ << QString::fromLocal8Bit("信噪比") << QString::fromLocal8Bit("陡度") << QString::fromLocal8Bit("清晰度") << QString::fromLocal8Bit("灰度最大值")
+ << QString::fromLocal8Bit("灰度最小值") << QString::fromLocal8Bit("PIQE")
+ << QString::fromLocal8Bit("0°二阶矩") << QString::fromLocal8Bit("90°二阶矩")
+ << QString::fromLocal8Bit("45°二阶矩") << QString::fromLocal8Bit("135°二阶矩")
+ << QString::fromLocal8Bit("RER");
+ ui.tableWidget->setVerticalHeaderLabels(header);
+
+ ui.tableWidget->setVerticalScrollMode(QAbstractItemView::ScrollPerPixel);
+ ui.tableWidget->setHorizontalScrollMode(QAbstractItemView::ScrollPerPixel);
+
+}
+
+void NIIRS::setTableWidgetValue(double val, int row, int col)
+{
+ ui.tableWidget->setItem(row, col, new QTableWidgetItem(QString::number(val, 10, 6)));
+}
+
+void NIIRS::setXmlInterface(QString strXmlInterface)
+{
+ mXmlInterface = strXmlInterface;
+
+ //qDebug() << "setXmlInterface:" << mXmlInterface;
+ std::cout << "setXmlInterface:" << mXmlInterface.toStdString();
+
+ //解析xml
+ if (mXmlInterface != "")
+ {
+ QFile fileXML(mXmlInterface);
+ if (fileXML.open(QIODevice::ReadOnly))
+ {
+ QDomDocument dom;
+ QString errorStr;
+ int errorLine = 0;
+ int errorCol = 0;
+ if (!dom.setContent(&fileXML, true, &errorStr, &errorLine, &errorCol))
+ {
+ ui.lineEdit_INFO->setText(errorStr + ", line: " + QString::number(errorLine) + ", col: " + QString::number(errorCol));
+ fileXML.close();
+ return;
+ }
+ else
+ {
+ QDomElement root = dom.documentElement();
+ QDomNode current = root;
+ QStack stack;
+ //遍历XML节点
+ while (1)
+ {
+ if (!current.isNull())
+ {
+ //查找并读取XML中参数
+ QDomElement element = current.toElement();
+ if (element.tagName() == "InputImg")
+ mXmlImgPath = element.firstChild().toCharacterData().data();
+ else if (element.tagName() == "TopLeftX")
+ mXmlTopLeftX = element.firstChild().toCharacterData().data();
+ else if (element.tagName() == "TopLeftY")
+ mXmlTopLeftY = element.firstChild().toCharacterData().data();
+ else if (element.tagName() == "BottomRightX")
+ mXmlBottomRightX = element.firstChild().toCharacterData().data();
+ else if (element.tagName() == "BottomRightY")
+ mXmlBottomRightY = element.firstChild().toCharacterData().data();
+ else if (element.tagName() == "ProductDir")
+ mXmlProductDir = element.firstChild().toCharacterData().data();
+ else if (element.tagName() == "ProductResultXML")
+ mResultXmlDir = element.firstChild().toCharacterData().data();
+
+ stack.push(current);
+ current = current.firstChild();
+ }
+ else if (!stack.count())
+ break;
+ else
+ current = stack.pop().nextSibling();
+ }
+ }
+ fileXML.close();
+ }
+ }
+ if (mXmlImgPath != "")
+ ui.lineEdit_img->setText(mXmlImgPath);
+
+ qDebug() << mXmlImgPath;
+ qDebug() << mXmlTopLeftX;
+ qDebug() << mXmlTopLeftY;
+ qDebug() << mXmlBottomRightX;
+ qDebug() << mXmlBottomRightY;
+ qDebug() << mXmlProductDir;
+
+ ui.lineEditLTx->setText(mXmlTopLeftX);
+ ui.lineEditLTy->setText(mXmlTopLeftY);
+ ui.lineEditRBx->setText(mXmlBottomRightX);
+ ui.lineEditRBy->setText(mXmlBottomRightY);
+
+}
+
+void NIIRS::calculateObjectivePara()
+{
+ ui.lineEdit_INFO->clear();
+ ui.tableWidget->clear();
+ ui.tableWidget->setColumnCount(0);
+ ui.tableWidget->setRowCount(0);
+
+ GDALAllRegister();
+ CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO"); //中文路径
+ //读取tiff
+ QString evaluateImg = ui.lineEdit_img->text();
+ if (evaluateImg.isEmpty())
+ {
+ //QMessageBox::warning(this, "warning", "The image to be evaluated is not opened");
+ ui.lineEdit_INFO->setText(QString::fromLocal8Bit("未打开待评估的图像"));
+ return;
+ }
+
+ string tiffPath = evaluateImg.toStdString();
+ GDALDataset* poDataset;
+ poDataset = (GDALDataset*)GDALOpen(tiffPath.c_str(), GA_ReadOnly);
+ if (NULL == poDataset)
+ {
+ GDALClose(poDataset);
+ return;
+ }
+ GDALRasterBand* pRasterBand;
+ GDALRasterBand* pClipBand;
+ vector pBuf;
+
+ //存储各波段参数
+ vectorvecGrayMax;
+ vectorvecGrayMin;
+
+ vectorvecLight;
+ vectorvecStdDev;
+ vectorvecContrastRadio;
+ vectorvecSNR;
+ vectorvecGradient;
+ vectorvecDefinition;
+ vectorvecPIQE;
+ vectorvecEntropy;
+ vector>vecASM;
+ vectortemp;
+
+ //poDataset->GetMetadata();
+
+ GDALRasterBand* pTestBand = poDataset->GetRasterBand(1);
+ GDALDataType dataType = pTestBand->GetRasterDataType();
+
+ if (dataType != GDT_UInt16)// && dataType != GDT_Byte
+ {
+ ui.lineEdit_INFO->setText(QString::fromLocal8Bit("输入影像数据类型无法用于计算客观参数"));
+ return;
+ }
+ //qDebug()<<"Dataset data type:GDT_UInt16";
+ int bandCount = poDataset->GetRasterCount();
+ int iwidth = pTestBand->GetXSize();
+ int iheight = pTestBand->GetYSize();
+
+ // 根据裁切范围确定裁切后的图像宽高
+// int topLeftX = mXmlTopLeftX.toInt();
+// int topLeftY = mXmlTopLeftY.toInt();
+// int bottomRightX = mXmlBottomRightX.toInt();
+// int bottomRightY = mXmlBottomRightY.toInt();
+
+ int topLeftX = ui.lineEditLTx->text().toInt();
+ int topLeftY = ui.lineEditLTy->text().toInt();
+ int bottomRightX = ui.lineEditRBx->text().toInt();
+ int bottomRightY = ui.lineEditRBy->text().toInt();
+
+ int clipWidth = bottomRightX - topLeftX;
+ int clipHeight = bottomRightY - topLeftY;
+
+ bool isCliped = false;
+
+ //逐波段计算参数
+ for (int i = 0; i < bandCount; i++)
+ {
+ unsigned short* bandBuf = new unsigned short[iwidth * iheight];
+ unsigned short* bandBufClip = new unsigned short[clipWidth * clipHeight];
+
+ vecASM.push_back(temp);
+ pRasterBand = poDataset->GetRasterBand(i + 1);
+
+ //在当前文件夹创建单波段裁剪图
+ QString tempTiffPath = QCoreApplication::applicationDirPath();
+ tempTiffPath = tempTiffPath + "/test.tiff";
+
+ //计算1.亮度、2.标准差，5.信噪比=Mean/StdDev
+ double bandMin = 0;
+ double bandMax = 0;
+ double bandMean = 0;
+ double bandStdDev = 0;
+
+ //clipWidth = 0;
+ //判断是否裁剪
+ if (clipWidth > 0 && clipHeight > 0)
+ {
+ isCliped = true;
+
+ //解析波段
+ if (pRasterBand->RasterIO(GF_Read, topLeftX, topLeftY, clipWidth, clipHeight, bandBufClip, clipWidth, clipHeight, dataType, 0, 0) != 0)
+ {
+ ui.lineEdit_INFO->setText(QString::fromLocal8Bit("解析裁剪数据出现错误, 波段:") + QString::number(i + 1));
+ continue;
+ }
+ pBuf.push_back(bandBufClip);
+
+ //计算1.亮度，2.标准差，5.信噪比
+ GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTIFF");
+ GDALDataset* clipDataset = poDriver->Create(tempTiffPath.toStdString().c_str(), clipWidth, clipHeight, 1, dataType, 0);
+ pClipBand = clipDataset->GetRasterBand(1);
+ pClipBand->RasterIO(GF_Write, 0, 0, clipWidth, clipHeight, bandBufClip, clipWidth, clipHeight, dataType, 0, 0);
+
+ pClipBand->ComputeStatistics(false, &bandMin, &bandMax, &bandMean, &bandStdDev, NULL, NULL);
+ vecLight.push_back(bandMean);
+ vecStdDev.push_back(bandStdDev);
+ vecSNR.push_back(bandMean / bandStdDev);
+
+ vecGrayMax.push_back(bandMax);
+ vecGrayMin.push_back(bandMin);
+
+ //4.计算图像熵
+ //std::cout<<"calculateTiffEntropy bandMax: "<RasterIO(GF_Read, 0, 0, iwidth, iheight, bandBuf, iwidth, iheight, dataType, 0, 0) != 0)
+ {
+ ui.lineEdit_INFO->setText(QString::fromLocal8Bit("输入数据解析错误, 波段:") + QString::number(i + 1));
+ continue;
+ }
+ pBuf.push_back(bandBuf);
+
+ //计算1.亮度，2.标准差，5.信噪比
+
+ pRasterBand->ComputeStatistics(false, &bandMin, &bandMax, &bandMean, &bandStdDev, NULL, NULL);
+ vecLight.push_back(bandMean);
+ vecStdDev.push_back(bandStdDev);
+ vecSNR.push_back(bandMean / bandStdDev);
+
+ vecGrayMax.push_back(bandMax);
+ vecGrayMin.push_back(bandMin);
+
+ //4.计算图像熵
+ std::cout << "calculateTiffEntropy bandMax: " << bandMax;
+ double tiffEntropy = calculateTiffEntropy(pBuf[i], iwidth * iheight, bandMax);
+ vecEntropy.push_back(tiffEntropy);
+
+ }
+
+
+ //数组转到mat，判断是否是裁剪过的
+ cv::Mat imageBand;
+ int imgW, imgH;
+ if (isCliped)
+ {
+ imgW = clipWidth;
+ imgH = clipHeight;
+ }
+ else
+ {
+ imgW = iwidth;
+ imgH = iheight;
+ }
+ if (dataType == GDT_Byte)
+ imageBand = cv::Mat(imgH, imgW, CV_8UC1, pBuf[i]);
+ else if (dataType == GDT_UInt16)
+ imageBand = cv::Mat(imgH, imgW, CV_16UC1, pBuf[i]);
+ else
+ {
+ ui.lineEdit_INFO->setText(QString::fromLocal8Bit("输入影像数据类型无法用于计算客观参数"));
+ continue;
+ }
+ //else if(dataType == GDT_UInt32)
+ // imageBand = cv::Mat(iheigth, iwidth, CV_32SC1, pBuf[i]);
+
+
+ //3.计算对比度
+ double contrastRadio = calculateContrastRatio(imageBand);
+ vecContrastRadio.push_back(contrastRadio);
+
+ //6.计算陡度
+ double sobel = gradient_Sobel(imageBand);
+ vecGradient.push_back(sobel);
+
+ //7.计算laplation
+ double laplation = definition_Laplacian(imageBand);
+ vecDefinition.push_back(laplation);
+
+ //计算角二阶矩、图像熵
+ cv::Mat dst_0degree, dst_90degree, dst_45degree, dst_135degree;
+ get_GLCM_0deg(imageBand, dst_0degree);
+ get_GLCM_90deg(imageBand, dst_90degree);
+ get_GLCM_45deg(imageBand, dst_45degree);
+ get_GLCM_135deg(imageBand, dst_135degree);
+ //qDebug() << "--calculate asm band" << i + 1;
+ //vecASM.push_back(temp);
+ //qDebug() << "-----0 degree-----";
+ vecASM[i].push_back(calculateASM(dst_0degree));
+ //qDebug() << "-----90 degree-----";
+ vecASM[i].push_back(calculateASM(dst_90degree));
+ //qDebug() << "-----45 degree-----";
+ vecASM[i].push_back(calculateASM(dst_45degree));
+ //qDebug() << "-----135 degree-----";
+ vecASM[i].push_back(calculateASM(dst_135degree));
+
+ //9.计算PIQE
+ double piqe = calculatePIQEScore(imageBand);
+ vecPIQE.push_back(piqe);
+
+ delete[] bandBufClip;
+ delete[] bandBuf;
+ }
+
+ initTableWidget(14, bandCount);
+ /*
+ * 亮度
+ * 标准差
+ * 对比度
+ * 熵
+ * 信噪比
+ * 陡度
+ * 清晰度
+ * PIQE
+ * 0°二阶矩,90°二阶矩,45°二阶矩,135°二阶矩
+ * --RER--
+ */
+ int currentBandCount = ui.comboBoxBand->currentIndex();
+ ui.lineEdit_SNR->setText(QString::number(vecSNR[currentBandCount], 10, 4));
+ ui.lineEdit_PIQE->setText(QString::number(vecPIQE[currentBandCount], 10, 4));
+
+ imgLight = vecLight[currentBandCount];
+ imgStdDev = vecStdDev[currentBandCount];
+ imgContrastRadio = vecContrastRadio[currentBandCount];
+ imgEntropy = vecEntropy[currentBandCount];
+ imgSNR = vecSNR[currentBandCount];
+ imgGradient = vecGradient[currentBandCount];
+ imgDefinition = vecDefinition[currentBandCount];
+ imgASM = (vecASM[0][currentBandCount] + vecASM[currentBandCount][1] + vecASM[currentBandCount][2] + vecASM[currentBandCount][3]) / 4;
+ imgPIQE = vecPIQE[currentBandCount];
+
+ for (int i = 0; i < bandCount; i++)
+ {
+ setTableWidgetValue(vecLight[i], 0, i);
+ setTableWidgetValue(vecStdDev[i], 1, i);
+ setTableWidgetValue(vecContrastRadio[i], 2, i);
+ setTableWidgetValue(vecEntropy[i], 3, i);
+ setTableWidgetValue(vecSNR[i], 4, i);
+ setTableWidgetValue(vecGradient[i], 5, i);
+ setTableWidgetValue(vecDefinition[i], 6, i);
+ setTableWidgetValue(vecGrayMax[i], 7, i);
+ setTableWidgetValue(vecGrayMin[i], 8, i);
+ setTableWidgetValue(vecPIQE[i], 9, i);
+ setTableWidgetValue(vecASM[i][0], 10, i);
+ setTableWidgetValue(vecASM[i][1], 11, i);
+ setTableWidgetValue(vecASM[i][2], 12, i);
+ setTableWidgetValue(vecASM[i][3], 13, i);
+ }
+
+ //qDebug() << "----------------------------------------------------------\n";
+ //释放内存
+ GDALClose(poDataset);
+
+ QString outTxt;
+ //在输出文件夹生成txt
+ if (mXmlProductDir != "")
+ {
+ QDir productDir(mXmlProductDir);
+ if (!productDir.exists())
+ {
+ ui.lineEdit_INFO->setText(QString::fromLocal8Bit("product.txt not exists"));
+ return;
+ }
+ outTxt = productDir.absoluteFilePath("product.txt");
+ //qDebug() << "product.txt: " << outTxt;
+ QFile txtFile(outTxt);
+ if (!txtFile.open(QIODevice::WriteOnly | QIODevice::Text | QIODevice::Truncate))//Truncate:覆盖文本///Append:在文本最后追加内容
+ return;
+ QDateTime time = QDateTime::currentDateTime();
+ QString currentTime = time.toString("yyyyMMddhhmmss");
+ QTextStream stream(&txtFile);
+ stream << QString::fromLocal8Bit("日期时间: %1\n").arg(currentTime);
+ stream << QString::fromLocal8Bit("文件路径: %1\n").arg(evaluateImg);
+ if (isCliped) {
+ stream << QString::fromLocal8Bit("图像宽度: %1, ").arg(clipWidth);
+ stream << QString::fromLocal8Bit("高度: %1\n").arg(clipHeight);
+ }
+ else {
+ stream << QString::fromLocal8Bit("图像宽度: %1, ").arg(iwidth);
+ stream << QString::fromLocal8Bit("高度: %1\n").arg(iheight);
+ }
+ for (int i = 0; i < bandCount; i++)
+ {
+ stream << QString::fromLocal8Bit("波段: %1\n").arg(i + 1);
+
+ stream << QString::fromLocal8Bit(" 图像灰度最大值: %1\n").arg(vecGrayMax[i]);
+ stream << QString::fromLocal8Bit(" 图像灰度最小值: %1\n").arg(vecGrayMin[i]);
+ stream << QString::fromLocal8Bit(" 图像亮 度: %1\n").arg(vecLight[i]);
+ stream << QString::fromLocal8Bit(" 图像标准差: %1\n").arg(vecStdDev[i]);
+ stream << QString::fromLocal8Bit(" 图像对比度: %1\n").arg(vecContrastRadio[i]);
+ stream << QString::fromLocal8Bit(" 图像熵: %1\n").arg(vecEntropy[i]);
+ stream << QString::fromLocal8Bit(" 图像信噪比: %1\n").arg(vecSNR[i]);
+ stream << QString::fromLocal8Bit(" 图像陡 度: %1\n").arg(vecGradient[i]);
+ stream << QString::fromLocal8Bit(" 图像清晰度: %1\n").arg(vecDefinition[i]);
+ stream << QString::fromLocal8Bit(" 图像PIQE: %1\n").arg(vecPIQE[i]);
+ stream << QString::fromLocal8Bit(" 图像0°二阶矩: %1\n").arg(vecASM[i][0]);
+ stream << QString::fromLocal8Bit(" 图像90°二阶矩: %1\n").arg(vecASM[i][1]);
+ stream << QString::fromLocal8Bit(" 图像45°二阶矩: %1\n").arg(vecASM[i][2]);
+ stream << QString::fromLocal8Bit(" 图像135°二阶矩: %1\n").arg(vecASM[i][3]);
+ }
+ stream << "\n";
+ txtFile.close();
+ ui.lineEdit_INFO->setText("result files save at:" + outTxt);
+ }
+ else
+ {
+ ui.lineEdit_INFO->setText(QString::fromLocal8Bit("product.txt输出文件夹为空或不存在"));
+ }
+}
+
+//void NIIRS::calculateEvaluationImage()
+//{
+// QString evaluateImgPath = ui.lineEdit_img->text();
+// // QStringList evaluateFiles = getAllFiles(evaluateImgPath, "BMP");
+// // if (evaluateFiles.size() == 0)
+// // {
+// // ui.lineEdit_INFO->setText("No BMP files in folder:" + evaluateImgPath);
+// // return;
+// // }
+// // QString textFileName = "D:\\A_testdata\\niirs_test\\tid2008\\a_distorted_images.txt";
+// // QFile file(textFileName);
+// // if (!file.open(QFile::Append | QFile::Text))
+// // {
+// // QMessageBox::warning(this, "warning", "Please select the correct file");
+// // return;
+// // }
+// // QTextStream textStream(&file);
+// // textStream << "imgName,imgLight,imgStdDev,imgContrastRadio,imgEntropy,imgSNR,imgGradient,imgDefinition,imgASM,imgPIQE,\n";
+// // for (int i = 0; i < evaluateFiles.size(); i++)
+// // {
+// // cv::Mat evaluateImg = cv::imread(evaluateFiles[i].toStdString());
+// cv::Mat evaluateImg = cv::imread(evaluateImgPath.toStdString());
+// cv::Mat imgGray;
+// if (evaluateImg.channels() == 3)
+// cv::cvtColor(evaluateImg, imgGray, CV_BGR2GRAY);
+// else
+// imgGray = evaluateImg;
+//
+// cv::Scalar meanScalar, stdDevScalar;
+// cv::meanStdDev(imgGray, meanScalar, stdDevScalar);
+// QVectorparas;
+// double imgLight = meanScalar.val[0];
+// double imgStdDev = stdDevScalar.val[0];
+// double imgContrastRadio = calculateContrastRatio(imgGray);
+// double imgEntropy = calculateEntropy(imgGray);
+// double imgSNR = calculateSNR(imgGray);
+// double imgSobel = gradient_Sobel(imgGray);
+// double imgLaplacian = definition_Laplacian(imgGray);
+// paras.append(imgLight);
+// paras.append(imgStdDev);
+// paras.append(imgContrastRadio);
+// paras.append(imgEntropy);
+// paras.append(imgSNR);
+// paras.append(imgSobel);
+// paras.append(imgLaplacian);
+//
+// cv::Mat dst_0, dst_90, dst_45, dst_135;
+// get_GLCM_0deg(imgGray, dst_0);
+// get_GLCM_90deg(imgGray, dst_90);
+// get_GLCM_45deg(imgGray, dst_45);
+// get_GLCM_135deg(imgGray, dst_135);
+// double degree_0 = calculateASM(dst_0);
+// double degree_90 = calculateASM(dst_90);
+// double degree_45 = calculateASM(dst_45);
+// double degree_135 = calculateASM(dst_135);
+// double imgASM = (degree_0 + degree_90 + degree_45 + degree_135) / 4.0;
+// double imgPIQE = calculatePIQEScore(imgGray);
+// paras.append(imgASM);
+// paras.append(imgPIQE);
+//
+// // QFileInfo fileInfo(evaluateFiles[i]);
+// // QString imgName = fileInfo.fileName();
+// // textStream << imgName<<",";
+// // for (int j = 0; j < paras.size(); j++)
+// // textStream << QString::number(paras[j], 10, 6)<<",";
+// // textStream << "\n";
+// // }
+// // file.close();
+// // qDebug() << "calculateEvaluationImage done";
+//}
+
+double NIIRS::calculateNiirsScore()
+{
+ calculateRER();
+
+ QString textGSDX = ui.lineEdit_GSDX->text();
+ QString textGSDY = ui.lineEdit_GSDY->text();
+ double GSDX, GSDY;
+ QString textRER = ui.lineEdit_RER->text();
+ QString text_H = ui.lineEdit_H->text();
+ double RER, MTFC_H;
+ QString textMTFC_G = ui.lineEdit_MTFC_G->text();
+ double MTFC_G;
+ QString textSNR = ui.lineEdit_SNR->text();
+ QString textPIQE = ui.lineEdit_PIQE->text();
+ double SNR, PIQE;
+ //1.检查GSD输入是否正确
+ if (textGSDX.toDouble() <= 0.0 || textGSDY.toDouble() <= 0.0)
+ {
+ ui.lineEdit_INFO->setText(QString::fromLocal8Bit("检查GSD输入"));
+ return 0.0;
+ }
+ else
+ {
+ GSDX = textGSDX.toDouble();
+ GSDY = textGSDY.toDouble();
+ ui.lineEdit_INFO->clear();
+ }
+ //2.检查RER输入是否正确
+ if (textRER.toDouble() == 0.0)
+ {
+ ui.lineEdit_INFO->setText(QString::fromLocal8Bit("请输入RER数据"));
+ return 0.0;
+ }
+ else
+ {
+ RER = textRER.toDouble();
+
+ ui.lineEdit_INFO->clear();
+ }
+ //3.检查MTFC_H输入是否正确
+ if (text_H.toDouble() < 0.0 || text_H.toDouble() > 1.90)
+ {
+ ui.lineEdit_INFO->setText(QString::fromLocal8Bit("检查几何过冲H输入"));
+ return 0.0;
+ }
+ else
+ {
+ MTFC_H = text_H.toDouble();
+ ui.lineEdit_INFO->clear();
+ }
+ //4.检查MTFC_G输入是否正确
+ if (textMTFC_G.toDouble() < 1.0 || textMTFC_G.toDouble() > 19.0)
+ {
+ ui.lineEdit_INFO->setText(QString::fromLocal8Bit("检查噪声增益G输入"));
+ return 0.0;
+ }
+ else
+ {
+ MTFC_G = textMTFC_G.toDouble();
+ ui.lineEdit_INFO->clear();
+ }
+ //5.检查SNR输入是否正确
+ if (textSNR.toDouble() < 0.0 || textSNR.toDouble() > 130.0)
+ {
+ ui.lineEdit_INFO->setText(QString::fromLocal8Bit("检查信噪比SNR输入"));
+ return 0.0;
+ }
+ else
+ {
+ SNR = textSNR.toDouble();
+ ui.lineEdit_INFO->clear();
+ }
+ //6.检查PIQE输入是否正确
+ if (textPIQE.toDouble() == 0.0)
+ {
+ ui.lineEdit_INFO->setText(QString::fromLocal8Bit("检查PIQE"));
+ return 0.0;
+ }
+ else
+ {
+ PIQE = textPIQE.toDouble();
+ ui.lineEdit_INFO->clear();
+ }
+
+ double NIIRS = 0.0;
+ double formulaPara = 0.0;
+ double GSD_inch;
+ if (ui.radioVisible->isChecked())
+ formulaPara = 10.251;
+ else if (ui.radioMidInfrared->isChecked() || ui.radioLonInfrared->isChecked())
+ formulaPara = 10.751;
+ //需要将米转换到英寸
+ GSD_inch = sqrtf(GSDX * GSDY) * 39.3700787402;
+
+ //RER>=0.9:a=3.32,b=1.559;RER<0.9:a=3.16,b=2.817
+ double para_a = RER >= 0.9 ? 3.32 : 3.16;
+ double para_b = RER >= 0.9 ? 1.559 : 2.817;
+ //qDebug() << "--a:" << a << ", b:" << b << ",para:" << para;
+ //可见10.251 - a*lgGSD + b*lgRER - 0.656*H - 0.334*(G/SNR)
+ //红外10.751 - a*lgGSD + b*lgRER - 0.656*H - 0.334*(G/SNR)
+ NIIRS = formulaPara - (para_a * log10(GSD_inch)) + (para_b * log10(RER)) - (0.656 * MTFC_H) - (0.334 * (MTFC_G / SNR));
+
+ ui.lineEdit_NIIRS->setText(QString::number(NIIRS, 10, 4));
+
+
+
+
+ QString outXMLPath = mResultXmlDir;
+ std::cout << outXMLPath.toStdString() << std::endl;
+ if (outXMLPath == "")
+ QMessageBox::information(this, tr("information"), tr("The Product Dir is Empty"));
+ else
+ {
+ //outXMLPath = outXMLPath + "/result.xml";
+ //ui.lineEdit_INFO->setText(outXMLPath);
+ std::cout << outXMLPath.toStdString() << std::endl;
+ QFile file(outXMLPath);
+ if (!file.open(QIODevice::WriteOnly | QIODevice::Truncate))
+ {//之前的内容被清空
+ qDebug() << QStringLiteral("打开") << outXMLPath << QStringLiteral("文件失败！");
+ return 0.0;
+ }
+ QDomDocument doc;
+ QDomProcessingInstruction processInstruction = doc.createProcessingInstruction("xml", "version=\"1.0\" encoding=\"UTF-8\"");
+ doc.appendChild(processInstruction);
+
+ //根节点
+ QDomElement root = doc.createElement("Root");
+ doc.appendChild(root);
+
+ {
+ //节点1 Metadata
+ QDomElement paraMetadata = doc.createElement("Metadata");
+ root.appendChild(paraMetadata);
+
+ QDomElement paraMeta1 = doc.createElement("Path");
+ QDomText paraMeta1Text = doc.createTextNode(ui.lineEdit_img->text());
+ paraMeta1.appendChild(paraMeta1Text);
+ paraMetadata.appendChild(paraMeta1);
+
+ QDomElement paraMeta2 = doc.createElement("Time");
+ QDomText paraMeta2Text = doc.createTextNode("2023.02.14");
+ paraMeta2.appendChild(paraMeta2Text);
+ paraMetadata.appendChild(paraMeta2);
+
+ QDomElement paraMetaLTx = doc.createElement("LeftTopX");
+ QDomText paraMetaLTxText = doc.createTextNode(mXmlTopLeftX);
+ paraMetaLTx.appendChild(paraMetaLTxText);
+ paraMetadata.appendChild(paraMetaLTx);
+
+ QDomElement paraMetaLTy = doc.createElement("LeftTopY");
+ QDomText paraMetaLTyText = doc.createTextNode(mXmlTopLeftY);
+ paraMetaLTy.appendChild(paraMetaLTyText);
+ paraMetadata.appendChild(paraMetaLTy);
+
+ QDomElement paraMetaRBx = doc.createElement("RightBottomX");
+ QDomText paraMetaRBxText = doc.createTextNode(mXmlBottomRightX);
+ paraMetaRBx.appendChild(paraMetaRBxText);
+ paraMetadata.appendChild(paraMetaRBx);
+
+ QDomElement paraMetaRBy = doc.createElement("RightBottomY");
+ QDomText paraMetaRByText = doc.createTextNode(mXmlBottomRightY);
+ paraMetaRBy.appendChild(paraMetaRByText);
+ paraMetadata.appendChild(paraMetaRBy);
+ }
+ {
+ //节点2 NIIRS参数及结果
+ QDomElement paraResultNiirs = doc.createElement("NIIRS");
+ root.appendChild(paraResultNiirs);
+
+ QDomElement paraNiirsPara = doc.createElement("Parameters");
+ paraResultNiirs.appendChild(paraNiirsPara);
+
+ QDomElement paraGSDX = doc.createElement("GSDX");
+ QDomText paraGSDXText = doc.createTextNode(ui.lineEdit_GSDX->text());
+ paraGSDX.appendChild(paraGSDXText);
+ paraNiirsPara.appendChild(paraGSDX);
+
+ QDomElement paraGSDY = doc.createElement("GSDY");
+ QDomText paraGSDYText = doc.createTextNode(ui.lineEdit_GSDY->text());
+ paraGSDY.appendChild(paraGSDYText);
+ paraNiirsPara.appendChild(paraGSDY);
+
+ QDomElement paraRER = doc.createElement("RER");
+ QDomText paraRERTest = doc.createTextNode(ui.lineEdit_RER->text());
+ paraRER.appendChild(paraRERTest);
+ paraNiirsPara.appendChild(paraRER);
+
+ QDomElement paraH = doc.createElement("H");
+ QDomText paraHText = doc.createTextNode(ui.lineEdit_H->text());
+ paraH.appendChild(paraHText);
+ paraNiirsPara.appendChild(paraH);
+
+ QDomElement paraG = doc.createElement("NoiseGain");
+ QDomText paraGText = doc.createTextNode(ui.lineEdit_MTFC_G->text());
+ paraG.appendChild(paraGText);
+ paraNiirsPara.appendChild(paraG);
+
+ QDomElement paraSNR = doc.createElement("SNR");
+ QDomText paraSNRText = doc.createTextNode(ui.lineEdit_SNR->text());
+ paraSNR.appendChild(paraSNRText);
+ paraNiirsPara.appendChild(paraSNR);
+
+ QDomElement paraNiirsPIQE = doc.createElement("PIQE");
+ QDomText paraNiirsPIQEText = doc.createTextNode(ui.lineEdit_PIQE->text());
+ paraNiirsPIQE.appendChild(paraNiirsPIQEText);
+ paraNiirsPara.appendChild(paraNiirsPIQE);
+
+
+
+ QDomElement paraNiirsScore = doc.createElement("Result");
+ paraResultNiirs.appendChild(paraNiirsScore);
+
+ QDomElement paraTiffType = doc.createElement("Type");
+ QDomText paraTiffTypeText;
+ if (ui.radioVisible->isChecked())
+ paraTiffTypeText = doc.createTextNode("Visible");
+ else if (ui.radioMidInfrared->isChecked())
+ paraTiffTypeText = doc.createTextNode("MiddleInfrared");
+ else if (ui.radioLonInfrared->isChecked())
+ paraTiffTypeText = doc.createTextNode("LongWaveInfrared");
+ paraTiffType.appendChild(paraTiffTypeText);
+ paraNiirsScore.appendChild(paraTiffType);
+
+ QDomElement paraScore = doc.createElement("NiirsScore");
+ QDomText paraScoreText = doc.createTextNode(ui.lineEdit_NIIRS->text());
+ paraScore.appendChild(paraScoreText);
+ paraNiirsScore.appendChild(paraScore);
+ }
+ {
+ //节点3 十个客观参数
+ QDomElement paraObjectiveParameter = doc.createElement("ObjectiveParameter");
+ root.appendChild(paraObjectiveParameter);
+
+ QDomElement paraLight = doc.createElement("ImageLight");
+ QDomText paraLightText = doc.createTextNode(QString::number(imgLight, 10, 6));
+ paraLight.appendChild(paraLightText);
+ paraObjectiveParameter.appendChild(paraLight);
+
+ QDomElement paraStdDev = doc.createElement("ImageStdDev");
+ QDomText paraStdDevText = doc.createTextNode(QString::number(imgStdDev, 10, 6));
+ paraStdDev.appendChild(paraStdDevText);
+ paraObjectiveParameter.appendChild(paraStdDev);
+
+ QDomElement paraContrastRadio = doc.createElement("ImageContrastRadio");
+ QDomText paraContrastRadioText = doc.createTextNode(QString::number(imgContrastRadio, 10, 6));
+ paraContrastRadio.appendChild(paraContrastRadioText);
+ paraObjectiveParameter.appendChild(paraContrastRadio);
+
+ QDomElement paraEntropy = doc.createElement("ImageEntropy");
+ QDomText paraEntropyText = doc.createTextNode(QString::number(imgEntropy, 10, 6));
+ paraEntropy.appendChild(paraEntropyText);
+ paraObjectiveParameter.appendChild(paraEntropy);
+
+ QDomElement paraSNR = doc.createElement("ImageSNR");
+ QDomText paraSNRText = doc.createTextNode(QString::number(imgSNR, 10, 6));
+ paraSNR.appendChild(paraSNRText);
+ paraObjectiveParameter.appendChild(paraSNR);
+
+ QDomElement paraRER = doc.createElement("ImageRER");
+ QDomText paraRERText = doc.createTextNode(QString::number(imgRER, 10, 6));
+ paraRER.appendChild(paraRERText);
+ paraObjectiveParameter.appendChild(paraRER);
+
+ QDomElement paraGradient = doc.createElement("ImageGradient");
+ QDomText paraGradientText = doc.createTextNode(QString::number(imgGradient, 10, 6));
+ paraGradient.appendChild(paraGradientText);
+ paraObjectiveParameter.appendChild(paraGradient);
+
+ QDomElement paraDefinition = doc.createElement("ImageDefinition");
+ QDomText paraDefinitionText = doc.createTextNode(QString::number(imgDefinition, 10, 6));
+ paraDefinition.appendChild(paraDefinitionText);
+ paraObjectiveParameter.appendChild(paraDefinition);
+
+ QDomElement paraASM = doc.createElement("ImageASM");
+ QDomText paraASMText = doc.createTextNode(QString::number(imgASM, 10, 6));
+ paraASM.appendChild(paraASMText);
+ paraObjectiveParameter.appendChild(paraASM);
+
+ QDomElement paraPIQE = doc.createElement("ImagePIQE");
+ QDomText paraPIQEText = doc.createTextNode(QString::number(imgPIQE, 10, 6));
+ paraPIQE.appendChild(paraPIQEText);
+ paraObjectiveParameter.appendChild(paraPIQE);
+ }
+ QTextStream outFile(&file);
+ doc.save(outFile, 4);//缩进4格
+ file.close();
+
+ ui.lineEdit_INFO->setText("result files save at:" + outXMLPath);
+ }
+
+ return NIIRS;
+}
+
+double NIIRS::calculateRER()
+{
+ ui.lineEdit_INFO->clear();
+
+ GDALAllRegister();
+ CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");
+
+ //QString inEdgeImg = mXmlImgPath;
+ QString inEdgeImg = ui.lineEdit_img->text();
+ //ui.lineEdit_INFO->setText(inEdgeImg);
+ //double imgRERValue = 0.0, imgHValue = 0.0;
+ cv::Mat imgTarget;
+
+ if (inEdgeImg == "")
+ {
+ ui.lineEdit_INFO->setText(QString::fromLocal8Bit("无刃边数据"));
+ return 0.0;
+ }
+
+ string tiffEdgePath = inEdgeImg.toStdString();
+ GDALDataset* poDataset;
+ poDataset = (GDALDataset*)GDALOpen(tiffEdgePath.c_str(), GA_ReadOnly);
+ if (poDataset == nullptr)
+ {
+ GDALClose(poDataset);
+ ui.lineEdit_INFO->setText(QString::fromLocal8Bit("刃边图像无法解析"));
+ return 0.0;
+ }
+
+ //GDALRasterBand* rasterBand = poDataset->GetRasterBand(1);
+ int currentBand = ui.comboBoxBand->currentIndex();
+ GDALRasterBand* rasterBand = poDataset->GetRasterBand(currentBand + 1);
+ GDALDataType dataType = rasterBand->GetRasterDataType();
+ int inImgWidth = rasterBand->GetXSize();
+ int inImgHeight = rasterBand->GetYSize();
+
+ //int topLeftX = mXmlTopLeftX.toInt();
+ //int topLeftY = mXmlTopLeftY.toInt();
+ //int bottomRightX = mXmlBottomRightX.toInt();
+ //int bottomRightY = mXmlBottomRightY.toInt();
+
+ int topLeftX = ui.lineEditLTx->text().toInt();
+ int topLeftY = ui.lineEditLTy->text().toInt();
+ int bottomRightX = ui.lineEditRBx->text().toInt();
+ int bottomRightY = ui.lineEditRBy->text().toInt();
+
+ int clipWidth = bottomRightX - topLeftX;
+ int clipHeight = bottomRightY - topLeftY;
+ if (clipWidth == 0 && clipHeight == 0)
+ {
+ clipWidth = inImgWidth;
+ clipHeight = inImgHeight;
+ }
+
+ //读取x方向刃边图
+ if (dataType == GDT_Byte)
+ {
+ unsigned char* bandBuf = new unsigned char[clipWidth * clipHeight];
+ if (rasterBand->RasterIO(GF_Read, topLeftX, topLeftY, clipWidth, clipHeight, bandBuf, clipWidth, clipHeight, dataType, 0, 0) != 0)
+ {
+ ui.lineEdit_INFO->setText(QString::fromLocal8Bit("输入数据解析错误"));
+ return -1.0;
+ }
+ imgTarget = cv::Mat(clipHeight, clipWidth, CV_8UC1, bandBuf);
+ }
+ else if (dataType == GDT_UInt16)
+ {
+ unsigned short* bandBuf = new unsigned short[clipWidth * clipHeight];
+
+ if (rasterBand->RasterIO(GF_Read, topLeftX, topLeftY, clipWidth, clipHeight, bandBuf, clipWidth, clipHeight, dataType, 0, 0) != 0)
+ {
+ ui.lineEdit_INFO->setText(QString::fromLocal8Bit("输入数据解析错误"));
+ return -1.0;
+ }
+ cv::Mat img = cv::Mat(clipHeight, clipWidth, CV_16UC1, bandBuf);
+ double min, max;
+ cv::minMaxLoc(img, &min, &max);
+ img.convertTo(imgTarget, CV_8UC1, 255.0 / max, 0);
+ }
+ else
+ {
+ ui.lineEdit_INFO->setText(QString::fromLocal8Bit("输入影像数据类型无法用于RER计算"));
+ return -1.0;
+ }
+ GDALClose(poDataset);
+
+ //根据两点，计算斜率k，判断为x方向或y方向
+ //double slopeK = 0.0;
+ //QPairedgeVals;
+ //if (slopeK >= -1 && slopeK <= 1)
+ // edgeVals = calculateERx(imgTarget);
+ //else
+ // edgeVals = calculateERy(imgTarget);
+ //imgRERValue = edgeVals.first;
+ //imgHValue = edgeVals.second;
+ //ui.lineEdit_RER->setText(QString::number(imgRERValue, 10, 4));
+ //ui.lineEdit_H->setText(QString::number(imgHValue, 10, 4));
+ cv::Mat imgGray, imgGrayX4;
+ if (imgTarget.channels() == 3)
+ cv::cvtColor(imgTarget, imgGray, CV_BGR2GRAY);
+ else
+ imgGray = imgTarget;
+
+ cv::pyrUp(imgGray, imgGrayX4);
+ cv::pyrUp(imgGrayX4, imgGrayX4);//重采样为原来4倍，用于1/4像元大小的超采样
+ cv::Mat imgThres;
+ cv::threshold(imgGrayX4, imgThres, 130, 255, cv::THRESH_BINARY);
+ QString thresholdPath = mXmlProductDir + "/threshold.bmp";
+ cv::imwrite(thresholdPath.toStdString().c_str(), imgThres);
+
+ cv::Mat imgCanny;
+ cv::Canny(imgThres, imgCanny, 100, 200);//canny边缘提取
+ QString cannyPath = mXmlProductDir + "/canny.bmp";
+ cv::imwrite(cannyPath.toStdString().c_str(), imgCanny);
+
+ vectorR_lines;
+ //cv::HoughLinesP(imgCanny, R_lines, 1, CV_PI / 180, 100, 25, 100);//提取直线
+ cv::HoughLinesP(imgCanny, R_lines, 1, CV_PI / 180, 50, 50, 10);//提取直线
+
+ QMapmap_rer_h;//key RER, value H
+
+ if (R_lines.size() == 0)
+ {
+ //qDebug() << "find no line in edge x";
+ QMessageBox::information(this, "information", "find no edge in image\n" + inEdgeImg);
+ //ui.lineEdit_INFO->setText(QString::fromLocal8Bit("无法从输入数据中提取刃边"));
+
+ return 0.0;
+ }
+ double maxRER = 0.0;
+
+ omp_set_num_threads(8);
+#pragma omp parallel for
+ for (int z = 0; z < R_lines.size(); z++)
+ {
+ cv::Vec4i p = R_lines[z];
+
+ cv::Point pointA = cv::Point(p[0], p[1]);
+ cv::Point pointB = cv::Point(p[2], p[3]);
+
+ QPairrer_h = calculate_rer_h(imgGrayX4, pointA, pointB);
+ if (rer_h.first > maxRER)
+ {
+ mRerEdgePointA = pointA;
+ mRerEdgePointB = pointB;
+ }
+ map_rer_h.insert(rer_h.first, rer_h.second);
+ }
+
+ QString edgeOutPath = mXmlProductDir + "/edges.bmp";
+ QString outTxt;
+ //在输出文件夹生成txt
+ QDir productDir(mXmlProductDir);
+ outTxt = productDir.absoluteFilePath("edges.txt");
+ QFile txtFile(outTxt);
+ txtFile.open(QIODevice::WriteOnly | QIODevice::Text | QIODevice::Truncate);//Truncate:覆盖文本///Append:在文本最后追加内容
+ QTextStream stream(&txtFile);
+ stream << QString("max rer edge: LTx:%1, LTy:%2; RBx:%3, RBy:%4\n")
+ .arg(mRerEdgePointA.x).arg(mRerEdgePointA.y).arg(mRerEdgePointB.x).arg(mRerEdgePointB.y);
+ cv::cvtColor(imgGrayX4, imgGrayX4, CV_GRAY2BGR);
+ for (int z = 0; z < R_lines.size(); z++)
+ {
+ cv::Vec4i p = R_lines[z];
+
+ cv::Point pointA = cv::Point(p[0], p[1]);
+ cv::Point pointB = cv::Point(p[2], p[3]);
+
+ cv::line(imgGrayX4, pointA, pointB, cv::Scalar(255, 255, 0));
+ stream << QString("edge:%1, LTx:%2, LTy:%3; RBx:%4, RBy:%5\n").arg(z).arg(p[0]).arg(p[1]).arg(p[2]).arg(p[3]);
+ }
+ stream << "\n";
+ txtFile.close();
+
+ cv::imwrite(edgeOutPath.toStdString().c_str(), imgGrayX4);
+
+ maxRER = map_rer_h.lastKey();
+ double maxH = map_rer_h.value(maxRER);
+
+ ui.lineEdit_RER->setText(QString::number(maxRER, 10, 4));
+ ui.lineEdit_H->setText(QString::number(maxH, 10, 4));
+
+ imgRER = maxRER;
+
+ return maxRER;
+
+
+ /////////////////////////////////////////////////////////////////////////////////////////////
+
+
+ //QString pathTargetX = mEdgeXPath;
+ //QString pathTargetY = mEdgeYPath;
+ ////cv::Mat imgTargetX, imgTargetY;
+ //if (pathTargetX == "" || pathTargetY == "")
+ //{
+ // ui.lineEdit_INFO->setText(QString::fromLocal8Bit("未打开刃边数据"));
+ // return 0.0;
+ //}
+
+ ////GDALAllRegister();
+ ////CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");
+ //string tiffXPath = pathTargetX.toStdString();
+ //string tiffYPath = pathTargetY.toStdString();
+ //GDALDataset* poDatasetX;
+ //GDALDataset* poDatasetY;
+ //poDatasetX = (GDALDataset*)GDALOpen(tiffXPath.c_str(), GA_ReadOnly);
+ //poDatasetY = (GDALDataset*)GDALOpen(tiffYPath.c_str(), GA_ReadOnly);
+ //if (NULL == poDatasetX)
+ //{
+ // GDALClose(poDatasetX);
+ // GDALClose(poDatasetY);
+ // ui.lineEdit_INFO->setText(QString::fromLocal8Bit("X方向刃边图像无法解析"));
+ // return -1.0;
+ //}
+ //if (NULL == poDatasetY)
+ //{
+ // GDALClose(poDatasetX);
+ // GDALClose(poDatasetY);
+ // ui.lineEdit_INFO->setText(QString::fromLocal8Bit("Y方向刃边图像无法解析"));
+ // return -1.0;
+ //}
+ //GDALRasterBand* rasterXBand = poDatasetX->GetRasterBand(1);
+ //GDALRasterBand* rasterYBand = poDatasetY->GetRasterBand(1);
+ //GDALDataType dataTypeX = rasterXBand->GetRasterDataType();
+ //GDALDataType dataTypeY = rasterYBand->GetRasterDataType();
+ //int xwidth = rasterXBand->GetXSize();
+ //int xheight = rasterXBand->GetYSize();
+ //int ywidth = rasterYBand->GetXSize();
+ //int yheight = rasterYBand->GetYSize();
+
+ //cv::Mat imgTargetX, imgTargetY;
+ ////读取x方向刃边图
+ //if (dataTypeX == GDT_Byte)
+ //{
+ // unsigned char* bandBuf = new unsigned char[xwidth * xheight];
+ // if (rasterXBand->RasterIO(GF_Read, 0, 0, xwidth, xheight, bandBuf, xwidth, xheight, dataTypeX, 0, 0) != 0)
+ // {
+ // ui.lineEdit_INFO->setText(QString::fromLocal8Bit("输入数据解析错误"));
+ // return -1.0;
+ // }
+ // imgTargetX = cv::Mat(xheight, xwidth, CV_8UC1, bandBuf);
+ //}
+ //else if (dataTypeX == GDT_UInt16)
+ //{
+ // unsigned short* bandBuf = new unsigned short[xwidth * xheight];
+ // if (rasterXBand->RasterIO(GF_Read, 0, 0, xwidth, xheight, bandBuf, xwidth, xheight, dataTypeX, 0, 0) != 0)
+ // {
+ // ui.lineEdit_INFO->setText(QString::fromLocal8Bit("输入数据解析错误"));
+ // return -1.0;
+ // }
+ // cv::Mat img = cv::Mat(xheight, xwidth, CV_16UC1, bandBuf);
+ // double min, max;
+ // cv::minMaxLoc(img, &min, &max);
+ // img.convertTo(imgTargetX, CV_8UC1, 255.0 / max, 0);
+ // //cv::imwrite("D:\\A_testdata\\niirs_test\\imgTargetX.png", imgTargetX);
+ //}
+ //else
+ //{
+ // ui.lineEdit_INFO->setText(QString::fromLocal8Bit("输入影像数据类型无法用于RER计算"));
+ // return -1.0;
+ //}
+ ////读取y方向刃边图
+ //if (dataTypeY == GDT_Byte)
+ //{
+ // unsigned char* bandBuf = new unsigned char[ywidth * yheight];
+ // if (rasterYBand->RasterIO(GF_Read, 0, 0, ywidth, yheight, bandBuf, ywidth, yheight, dataTypeY, 0, 0) != 0)
+ // {
+ // ui.lineEdit_INFO->setText(QString::fromLocal8Bit("输入数据解析错误"));
+ // return -1.0;
+ // }
+ // imgTargetY = cv::Mat(yheight, ywidth, CV_8UC1, bandBuf);
+ //}
+ //else if (dataTypeY == GDT_UInt16)
+ //{
+ // unsigned short* bandBuf = new unsigned short[ywidth * yheight];
+ // if (rasterYBand->RasterIO(GF_Read, 0, 0, ywidth, yheight, bandBuf, ywidth, yheight, dataTypeY, 0, 0) != 0)
+ // {
+ // ui.lineEdit_INFO->setText(QString::fromLocal8Bit("输入数据解析错误"));
+ // return -1.0;
+ // }
+ // cv::Mat img = cv::Mat(yheight, ywidth, CV_16UC1, bandBuf);
+ // double min, max;
+ // cv::minMaxLoc(img, &min, &max);
+ // img.convertTo(imgTargetY, CV_8UC1, 255.0 / max, 0);
+ // //cv::imwrite("D:\\A_testdata\\niirs_test\\imgTargetY.png", imgTargetY);
+ //}
+ //else
+ //{
+ // ui.lineEdit_INFO->setText(QString::fromLocal8Bit("输入影像数据类型无法用于RER计算"));
+ // return -1.0;
+ //}
+
+ //double imgRER, imgH;
+ //QPairedgeXVals = calculateERx(imgTargetX);
+ //imgRER = edgeXVals.first;
+ //imgH = edgeXVals.second;
+ //ui.lineEdit_RER->setText(QString::number(imgRER, 10, 4));
+ //ui.lineEdit_H->setText(QString::number(imgH, 10, 4));
+
+ //qDebug() << "--6.img RER: " << imgRER;
+ //qDebug() << "----img H:" << imgH;
+
+ //GDALClose(poDatasetX);
+ //GDALClose(poDatasetY);
+
+ //return imgRER;
+}
+
+double NIIRS::calculateMeanStdDev(cv::Mat img)
+{
+ GDALAllRegister();
+ CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO"); //中文路径
+ //gdal读取待评价影像，逐波段计算平均值计算亮度
+ QString evaluateImg = ui.lineEdit_img->text();
+ if (evaluateImg.isEmpty())
+ {
+ QMessageBox::warning(this, "warning", "The image to be evaluated is not opened");
+ return -1;
+ }
+ string tiffPath = evaluateImg.toStdString();
+ GDALDataset* hSrcDS;
+ hSrcDS = (GDALDataset*)GDALOpen(tiffPath.c_str(), GA_ReadOnly);
+ if (NULL == hSrcDS)
+ {
+ GDALClose(hSrcDS);
+ return -1;
+ }
+ GDALRasterBand* rasterBand;
+ qDebug() << "--Img Light&StdDev";
+ for (int i = 0; i < hSrcDS->GetRasterCount(); i++)
+ {
+ rasterBand = hSrcDS->GetRasterBand(i + 1);
+ double bandMin = 0;
+ double bandMax = 0;
+ double bandMean = 0;
+ double bandStdDev = 0;
+ rasterBand->ComputeStatistics(false, &bandMin, &bandMax, &bandMean, &bandStdDev, NULL, NULL);
+ qDebug() << "---band:" << i + 1 << endl << "----Light:" << bandMean << endl << "----StdDev:" << bandStdDev;
+ }
+ delete hSrcDS;
+ return 0.0;
+}
+
+//float NIIRS::calculateStdDev(cv::Mat img)
+//{
+// //转为灰度图
+// cv::Mat imgGray;
+// if (img.channels() == 3)
+// cv::cvtColor(img, imgGray, CV_BGR2GRAY);
+// else
+// imgGray = img;
+//
+// cv::Scalar meanScalar, stdDevScalar;
+// cv::meanStdDev(imgGray, meanScalar, stdDevScalar);
+//
+// float imgStdDev = stdDevScalar.val[0];
+// return imgStdDev;
+//}
+
+double NIIRS::calculateTiffEntropy(unsigned short* bandBuf, int bufSize, int maxVal)
+{
+ //获取最大值，赋值temp，用于计算每级亮度累加值，使用vector
+ //1.计算每个亮度累加值
+ //遍历数组，每个像元的亮度使对应vector位置的数量+1
+ QVectortemp(maxVal + 1);
+ std::cout << "max:" << maxVal + 1 << std::endl;
+ for (int i = 0; i < bufSize; i++)
+ {
+ auto pos = bandBuf[i];
+ int tempAddr = pos / 1;
+ temp[tempAddr] = temp[tempAddr] + 1;
+ }
+
+ //2.计算每个亮度的概率
+ for (int i = 0; i < temp.size(); i++)
+ temp[i] = temp[i] / double(bufSize);
+ //3.计算图像信息熵
+ float result = 0;
+ for (int i = 0; i < temp.size(); i++)
+ {
+ if (temp[i] == 0.0)
+ result = result;
+ else
+ result = result - temp[i] * (std::log(temp[i]) / std::log(2.0));//公式
+ }
+
+ return result;
+}
+
+//double NIIRS::calculateEntropy(cv::Mat img)
+//{
+// double temp[256] = { 0.0 };
+//
+// //计算每个像素的累积值
+// for (int m = 0; m < img.rows; m++)
+// {
+// const uchar* t = img.ptr(m);
+// for (int n = 0; n < img.cols; n++)
+// {
+// int i = t[n];
+// temp[i] = temp[i] + 1;
+// }
+// }
+// //计算每个像素的概率
+// for (int i = 0; i < 256; i++)
+// {
+// temp[i] = temp[i] / (img.rows * img.cols);
+// }
+// float result = 0;
+// //计算图像信息熵
+// for (int i = 0; i < 256; i++)
+// {
+// if (temp[i] == 0.0)
+// result = result;
+// else
+// result = result - temp[i] * (log(temp[i]) / log(2.0));
+// }
+//
+// return result;
+//}
+
+double NIIRS::calculateContrastRatio(cv::Mat img)
+{
+ cv::Mat imgGray;
+ if (img.channels() == 3)
+ cv::cvtColor(img, imgGray, CV_BGR2GRAY);
+ else
+ imgGray = img;
+
+ cv::Mat imgExtend;
+ cv::copyMakeBorder(imgGray, imgExtend, 1, 1, 1, 1, cv::BORDER_REPLICATE);//以边缘值四周各扩充1像元
+ imgExtend = imgExtend / 1.0;
+
+ int row = imgExtend.rows;
+ int col = imgExtend.cols;
+ double bb = 0.0;
+ for (int i = 1; i < row - 1; i++)
+ {
+ for (int j = 1; j < col - 1; j++)
+ {
+ double center = imgExtend.at(i, j);
+ bb += pow(center - imgExtend.at(i, j - 1), 2) + pow(center - imgExtend.at(i, j + 1), 2)
+ + pow(center - imgExtend.at(i - 1, j), 2) + pow(center - imgExtend.at(i + 1, j), 2);
+ }
+ }
+ double cr = bb / (4 * (row - 2) * (col - 2) + 3 * (2 * (row - 2) + 2 * (col - 2)) + 2 * 4);
+ return cr;
+}
+
+//double NIIRS::calculateSNR(cv::Mat img)
+//{
+// cv::Mat imgGray;
+// if (img.channels() == 3)
+// cv::cvtColor(img, imgGray, CV_BGR2GRAY);
+// else
+// imgGray = img;
+//
+// cv::Scalar meanScalar, stdDevScalar;
+// cv::meanStdDev(imgGray, meanScalar, stdDevScalar);
+// double snr = meanScalar.val[0] / stdDevScalar.val[0];
+// return snr;
+//}
+
+//float NIIRS::SFRCalculation(cv::Mat& ROI, double gamma)
+//{
+// if (ROI.empty())
+// {
+// std::cerr << "Open the ROI image error" << std::endl;
+// return 0;
+// }
+// int height = ROI.rows, width = ROI.cols;
+// //进行伽马解码
+// de_Gamma(ROI, gamma);
+// int i, j;
+// double slope = 0, intercept = 0;
+// //中心质心偏移Center centroid offset
+// double CCoffset = 0;
+// std::vector y_shifts(height);
+// //计算图像中心的质心和质心之间的偏移
+// std::vector Cen_Shifts = CentroidFind(ROI, y_shifts, &CCoffset);
+// if (Cen_Shifts.empty())
+// return 0;
+// //斜边拟合，线性回归
+// SLR(Cen_Shifts, y_shifts, &intercept, &slope);
+// //将数据行数截断为最大斜率周期，该周期将具有整数个完整相位旋转
+// ReduceRows(slope, &height);
+// //将CCoffset更新为图像的原始中点与计算的参考中点之间的偏移
+// CCoffset = CCoffset + 0.5 + intercept - width / 2;
+// //进行过采样
+// int SamplingLen = width * 4;
+// std::vector OverSamplingData = OverSampling(ROI, slope, CCoffset, height, width, &SamplingLen);
+// //利用汉明窗对数据两侧的纹波信号进行滤波
+// OverSamplingData = HammingWindows(OverSamplingData, SamplingLen);
+// //离散傅里叶变换
+// DFT(OverSamplingData, SamplingLen);
+// width = int(SamplingLen / 4);
+// double maxData = 0;
+// for (i = 0; i < SamplingLen; ++i)
+// {
+// if (OverSamplingData[i] != 0)
+// {
+// maxData = OverSamplingData[i];
+// break;
+// }
+// }
+// for (int i = 0; i < SamplingLen; ++i)
+// OverSamplingData[i] /= maxData;
+// vector >vec_mtf;
+// for (i = 0; i <= width; ++i)
+// {
+// double frequency = (double)i / width;
+// pairmtf_single;
+// mtf_single.first = frequency;
+// mtf_single.second = OverSamplingData[i];
+// vec_mtf.push_back(mtf_single);
+// //mtf_file << frequency << "," << OverSamplingData[i] << "\n";
+// }
+//
+//
+// return 0.0;
+//}
+
+//void NIIRS::DFT(std::vector& data, int size)
+//{
+// int i, j;
+// std::complex* arr = new std::complex[size];
+// for (i = 0; i < size; ++i)
+// {
+// arr[i] = std::complex(data[i], 0);
+// }
+//
+// for (i = 0; i < size / 2.0; ++i)
+// {
+// std::complex temp = 0;
+// for (j = 0; j < size; ++j)
+// {
+// double w = 2 * 3.1415 * i * j / size;
+// std::complex deg(cos(w), -sin(w));
+// temp += arr[j] * deg;
+// }
+// data[i] = sqrt(temp.real() * temp.real() + temp.imag() * temp.imag());
+// }
+//}
+
+//std::vector NIIRS::HammingWindows(std::vector& deSampling, int SamplingLen)
+//{
+// int i, j;
+// std::vector tempData(SamplingLen);
+//
+// //We want to shift the peak data to the center of line spread function data
+// //Because we will do the hamming window later, this will keep the important data away from filtering
+// //In case there are two peaks, we use two variable to record the peak data position
+// int L_location = -1, R_location = -1;
+//
+// double SamplingMax = 0;
+// for (i = 0; i < SamplingLen; ++i)
+// {
+// if (fabs(deSampling[i]) > fabs(SamplingMax))
+// SamplingMax = deSampling[i];
+// }
+//
+// for (i = 0; i < SamplingLen; ++i)
+// {
+// if (deSampling[i] == SamplingMax)
+// {
+// if (L_location < 0) { L_location = i; }
+// R_location = i;
+// }
+// }
+//
+// //the shift amount
+// int PeakOffset = (R_location + L_location) / 2 - SamplingLen / 2;
+//
+// if (PeakOffset)
+// {
+// for (i = 0; i < SamplingLen; ++i)
+// {
+// int newIndex = i - PeakOffset;
+// if (newIndex >= 0 && newIndex < SamplingLen) { tempData[newIndex] = deSampling[i]; }
+// }
+// }
+// else
+// {
+// for (i = 0; i < SamplingLen; ++i) { tempData[i] = deSampling[i]; }
+// }
+//
+// //do the hamming window filtering
+// for (int i = 0; i < SamplingLen; ++i)
+// {
+// tempData[i] = tempData[i] * (0.54 - 0.46 * cos(2 * CV_PI * i / (SamplingLen - 1)));
+// }
+//
+// return tempData;
+//}
+
+//std::vector NIIRS::OverSampling(cv::Mat& Src, double slope, double CCoffset, int height, int width, int* SamplingLen)
+//{
+// std::vector RowShifts(height, 0);
+//
+// int i, j, k;
+// int halfY = height >> 1;
+//
+// //calculate the pixel shift of each row to align the centroid of each row as close as possible
+// for (i = 0; i < height; ++i) { RowShifts[i] = (double)(i - halfY) / slope + CCoffset; }
+//
+// //DataMap -> to record the index of pixel after shift
+// //Datas -> to record the original pixel value
+// std::vector DataMap(height * width, 0);
+// std::vector Datas(height * width, 0);
+//
+// for (i = 0, k = 0; i < height; ++i)
+// {
+// int baseindex = width * i;
+// for (j = 0; j < width; ++j)
+// {
+// DataMap[baseindex + j] = j - RowShifts[i];
+// Datas[baseindex + j] = int(Src.at(i, j));
+// }
+// }
+//
+// std::vector SamplingBar(*SamplingLen, 0);
+// std::vector MappingCount(*SamplingLen, 0);
+//
+// //Start to mapping the original data to 4x sampling data and record the count of each pixel in 4x sampling data
+// for (i = 0; i < height * width; ++i)
+// {
+// int MappingIndex = static_cast(4 * DataMap[i]);
+// if (MappingIndex >= 0 && MappingIndex < *SamplingLen)
+// {
+// SamplingBar[MappingIndex] = SamplingBar[MappingIndex] + Datas[i];
+// MappingCount[MappingIndex]++;
+// }
+// }
+//
+// //average the pixel value in 4x sampling data, if the pixel value in pixel is zero, copy the value of close pixel
+// for (i = 0; i < *SamplingLen; ++i) {
+// j = 0;
+// k = 1;
+// if (MappingCount[i] == 0) {
+//
+// if (i == 0) {
+// while (!j)
+// {
+// if (MappingCount[i + k] != 0)
+// {
+// SamplingBar[i] = SamplingBar[i + k] / ((double)MappingCount[i + k]);
+// j = 1;
+// }
+// else ++k;
+// }
+// }
+// else {
+// while (!j && ((i - k) >= 0))
+// {
+// if (MappingCount[i - k] != 0)
+// {
+// SamplingBar[i] = SamplingBar[i - k]; /* Don't divide by counts since it already happened in previous iteration */
+// j = 1;
+// }
+// else ++k;
+// }
+// if ((i - k) < 0)
+// {
+// k = 1;
+// while (!j)
+// {
+// if (MappingCount[i + k] != 0)
+// {
+// SamplingBar[i] = SamplingBar[i + k] / ((double)MappingCount[i + k]);
+// j = 1;
+// }
+// else ++k;
+// }
+// }
+// }
+// }
+// else
+// SamplingBar[i] = (SamplingBar[i]) / ((double)MappingCount[i]);
+// }
+//
+// // reduce the length of sampling data
+// // because the datas at the edge are only matters, we truncate the data close to the two side of sampling data
+// // which has very small contribution to the result
+//
+// //truncating the data smaller than 10% and bigger than 90% of original length
+// int originalSamplingLen = *SamplingLen;
+// *SamplingLen = *SamplingLen * 0.8;
+//
+// std::vector deSampling(*SamplingLen, 0);
+// //derivative sampling data(which is ESF) to get the line spread function
+// for (i = originalSamplingLen * 0.1, j = 1; i < originalSamplingLen, j < *SamplingLen; ++i, ++j)
+// {
+// deSampling[j] = SamplingBar[i + 1] - SamplingBar[i];
+// }
+//
+// return deSampling;
+//}
+
+//void NIIRS::ReduceRows(double slope, int* ImgHeight)
+//{
+// double tempSlope = fabs(slope);
+// int cycs = (*ImgHeight) / tempSlope;
+// if (tempSlope * cycs < *ImgHeight) { *ImgHeight = tempSlope * cycs; }
+//}
+
+//void NIIRS::SLR(std::vector& Cen_Shifts, std::vector& y_shifts, double* a, double* b)
+//{
+// //a -> intercept, b->slope of slanted edge
+//
+// int y_size = y_shifts.size();
+// //using simple linear regression to solve equation and get slope and intercept
+// double xsquare = 0, xavg = 0, yavg = 0;
+// int i;
+// for (i = 0; i < y_size; ++i)
+// {
+// yavg += y_shifts[i];
+// xavg += Cen_Shifts[i];
+// }
+// xavg /= (double)y_size;
+// yavg /= (double)y_size;
+//
+// //simple linear regession
+// for (i = 0; i < y_size; ++i)
+// {
+// double temp = (Cen_Shifts[i] - xavg);
+// *b += temp * (y_shifts[i] - yavg);
+// xsquare += temp * temp;
+// }
+//
+// *b /= xsquare;
+// *a = yavg - (*b) * xavg;
+//}
+
+//std::vector NIIRS::CentroidFind(cv::Mat& Src, std::vector& y_shifts, double* CCoffset)
+//{
+// std::vector Cen_Shifts(Src.rows);
+// int i, j, height = Src.rows, width = Src.cols;
+//
+// cv::Mat tempSrc(Src.size(), CV_8UC1);
+//
+// //Do the bilaterFilter on template ROI image to make sure we can find the slanted edge more acurrate
+// cv::bilateralFilter(Src, tempSrc, 3, 200, 200);
+//
+// // calculate the centroid of every row in ROI
+// // centroid formuld: Total moments/Total amount
+// for (i = 0; i < height; ++i)
+// {
+// double molecule = 0, denominator = 0, temp = 0;
+// uchar* tempSrcPtr = tempSrc.ptr(i);
+// for (j = 0; j < width - 1; ++j)
+// {
+// temp = (double)tempSrcPtr[j + 1] - (double)tempSrcPtr[j];
+// molecule += temp * j;
+// denominator += temp;
+// }
+// Cen_Shifts[i] = molecule / denominator;
+// }
+//
+// //Eliminate the noise far from slant edge(+/- 10 pixels)
+// tempSrc = Src.clone();
+// cv::GaussianBlur(Src, Src, cv::Size(3, 3), 0);
+// for (i = 0; i < height; ++i)
+// {
+// uchar* SrcPtr = Src.ptr(i);
+// uchar* tempPtr = tempSrc.ptr(i);
+// for (j = int(Cen_Shifts[i]) - 10; j < int(Cen_Shifts[i]) + 10; ++j)
+// {
+// SrcPtr[j] = tempPtr[j];
+// }
+// }
+//
+// // check whether the edge in image is too close to the image corners
+// if (Cen_Shifts[0] < 2.0 || width - Cen_Shifts[0] < 2.0)
+// {
+// std::cerr << "The edge in ROI is too close to the image corners" << std::endl;
+// return {};
+// }
+//
+// if (Cen_Shifts[height - 1] < 2 || width - Cen_Shifts[height - 1] < 2)
+// {
+// std::cerr << "The edge in ROI is too close to the image corners" << std::endl;
+// return {};
+// }
+//
+// int half_ySize = height / 2;
+// int CC = Cen_Shifts[half_ySize]; //Centroid of image center
+// *CCoffset = CC;
+// for (i = 0; i < height; ++i)
+// {
+// Cen_Shifts[i] -= CC; //Calculate the Shifts between Centroids and Centroid of image center
+// y_shifts[i] = i - half_ySize; //Calculate the shifts between height of image center and each row
+// }
+//
+// return Cen_Shifts;
+//}
+
+//void NIIRS::de_Gamma(cv::Mat& Src, double gamma)
+//{
+// if (Src.channels() != 1) { return; }
+//
+// for (int i = 0; i < Src.rows; ++i)
+// {
+// uchar* SrcP = Src.ptr(i);
+// for (int j = 0; j < Src.cols; ++j)
+// {
+// SrcP[j] = 255 * (pow((double)SrcP[j] / 255, 1 / gamma));
+// }
+// }
+//}
+
+QPair NIIRS::calculate_rer_h(cv::Mat img, cv::Point pointA, cv::Point pointB)
+{
+ QPairRER_H = QPair();
+ cv::Point pointLeft, pointRight;
+ pointLeft = pointA;
+ pointRight = pointB;
+ if (pointLeft.x > pointRight.x)
+ {
+ cv::Point temp = pointLeft;
+ pointLeft = pointRight;
+ pointRight = temp;
+ }
+
+ //pair中:vector 存放所有相同距离时像元值，之后求平均
+ //pair中:int 表示离刃边距离
+ QMap>qmapESF;
+ vectorvecTemp;
+ vecTemp.push_back(127);
+ for (int i = -20; i < 22; i++)
+ qmapESF.insert(i, vecTemp);
+
+
+ cv::Mat mask;
+ cv::Point pt1, pt2;
+ double k = (double(pointRight.y) - pointLeft.y) / (double(pointRight.x) - pointLeft.x);
+ if (k < -1 || k > 1)
+ {
+ cv::Point pointBottom, pointTop;
+ pointBottom = pointA;
+ pointTop = pointB;
+
+ if (pointBottom.y > pointTop.y)
+ {
+ cv::Point temp = pointBottom;
+ pointBottom = pointTop;
+ pointTop = temp;
+ }
+
+ int rangeL = 0, rangeT = 0, rangeR = 0, rangeB = 0;
+ int cutImgTopx = 0, cutImgTopy = 0, cutImgBotx = 0, cutImgBoty = 0;
+ if (pointBottom.x == pointTop.x)
+ {
+ rangeL = pointBottom.x - 20;
+ rangeR = pointBottom.x + 20 + 4;
+ rangeT = pointTop.y;
+ rangeB = pointTop.y;
+
+ cutImgTopx = 20;
+ cutImgTopy = 0;
+ cutImgBotx = pointTop.x - pointBottom.x + 20;
+ cutImgBoty = pointTop.y - pointBottom.y;
+ }
+ else if (k > 0)
+ {
+ double k2 = 1 / k;
+ int deltay = 20, deltax;
+ while (25 > pow((k * deltay), 2) + pow(double(deltay), 2))
+ deltay = deltay + 20;
+ deltax = ceil(deltay / k2);
+
+ rangeL = pointBottom.x - deltax;
+ rangeR = pointTop.x + deltax;
+ rangeT = pointTop.y + deltay + 4;
+ rangeB = pointBottom.y - deltay;
+
+ cutImgTopx = pointBottom.x - rangeL;
+ cutImgTopy = pointBottom.y - rangeB;
+ cutImgBotx = pointTop.x - rangeL;
+ cutImgBoty = pointTop.y - rangeB;
+ }
+ else//k<0
+ {
+ double k2 = 1 / k;
+ int deltay = 20, deltax;
+ while (25 > pow((k * deltay), 2) + pow(double(deltay), 2))
+ deltay = deltay + 20;
+ deltax = ceil(deltay / -k2);
+ //qDebug() << deltay << deltay / k2 << "->" << deltax;
+ rangeL = pointTop.x - deltax;
+ rangeR = pointBottom.x + deltax;
+ rangeB = pointBottom.y - deltay;
+ rangeT = pointTop.y + deltay + 4;
+
+ cutImgBotx = pointTop.x - rangeL;
+ cutImgBoty = pointTop.y - rangeB;
+ cutImgTopx = pointBottom.x - rangeL;
+ cutImgTopy = pointBottom.y - rangeB;
+ }
+
+ if (rangeL - rangeR == 0)
+ return RER_H;
+ if (rangeL < 0)
+ rangeL = 0;
+ if (rangeR > img.cols)
+ rangeR = img.cols;
+ if (rangeT > img.rows)
+ rangeT = img.rows;
+ if (rangeB < 0)
+ rangeB = 0;
+ //qDebug() << QString("Y rangeB:%1,rangeT:%2; rangeL:%3,rangeR:%4").arg(rangeB).arg(rangeT).arg(rangeL).arg(rangeR);
+ //QString outBmpPath = outPath + "/" + filename;
+ cv::Mat miniEdgeX = img(cv::Range(rangeB, rangeT), cv::Range(rangeL, rangeR));
+ cv::Mat evalueImg;
+ miniEdgeX.copyTo(evalueImg);
+
+ pt1 = cv::Point(cutImgTopx, cutImgTopy);
+ pt2 = cv::Point(cutImgBotx, cutImgBoty);
+ cv::Mat temp(evalueImg.rows, evalueImg.cols, evalueImg.type(), cv::Scalar(0));
+ mask = temp;
+ int rows = evalueImg.rows;
+ int cols = evalueImg.cols;
+ for (int row = 0; row < rows; row++)
+ {
+ for (int col = 0; col < cols; col++)
+ {
+ //当前采样点
+ cv::Point samplingPoint(col, row);
+ //用于计算当前点到刃边距离
+ arraysampling_to_edge = { samplingPoint, pt1, pt2 };
+ //求当前点到刃边的垂足，获取垂足的col，用于比较确认dis的正负
+ cv::Point2d samplingPointFoot = calculateFootPoint(sampling_to_edge);
+ if (samplingPointFoot.y < pt1.y || samplingPointFoot.y > pt2.y)
+ continue;
+ //采样点距离刃边的距离
+ int dis = pointDistance(samplingPointFoot, samplingPoint);
+ if (dis > 21)
+ continue;
+ mask.at(row, col) = evalueImg.at(row, col);
+
+ if (samplingPointFoot.x > col)
+ dis = 0 - dis;
+
+ // int value = evalueImg.at(row, col);
+ //存入vecESF中
+ //遍历vecESF，寻找相同距离的值
+ // if (qmapESF.contains(dis))
+ // {
+ vectorpixValues = qmapESF.value(dis);
+ pixValues.push_back(evalueImg.at(row, col));
+ qmapESF.insert(dis, pixValues);
+ // }
+ // else
+ // {
+ // vectorpixValues;
+ // pixValues.push_back(evalueImg.at(row, col));
+ // qmapESF.insert(dis, pixValues);
+ // }
+ }
+ }
+ }
+ else
+ {
+ int rangeL = 0, rangeT = 0, rangeR = 0, rangeB = 0;
+ int cutImgLTx = 0, cutImgLTy = 0, cutImgRBx = 0, cutImgRBy = 0;
+ if (k == 0)
+ {
+ rangeL = pointLeft.x;
+ rangeR = pointRight.x;
+ rangeT = pointRight.y + 20 + 4;
+ rangeB = pointLeft.y - 20;
+
+ cutImgLTx = pointLeft.x - pointLeft.x;
+ cutImgLTy = pointLeft.y + 20 - pointLeft.y;
+ cutImgRBx = pointRight.x - pointLeft.x;
+ cutImgRBy = pointRight.y + 20 - pointLeft.y;
+ }
+ else if (k > 0)
+ {
+ double k2 = 1 / k;
+ int deltaY = 20, deltaX;
+ while (25 > pow((k * deltaY), 2) + pow(double(deltaY), 2))
+ deltaY = deltaY + 20;
+ deltaX = ceil(deltaY / k2);
+
+ rangeL = pointLeft.x - deltaX;
+ rangeR = pointRight.x + deltaX;
+ rangeT = pointRight.y + deltaY + 4;
+ rangeB = pointLeft.y - deltaY;
+
+ cutImgLTx = pointLeft.x - rangeL;
+ cutImgLTy = pointLeft.y - rangeB;
+ cutImgRBx = pointRight.x - rangeL;
+ cutImgRBy = pointRight.y - rangeB;
+ }
+ else//k<0
+ {
+ double k2 = 1 / k;
+ int deltay = 20, deltax;
+ while (25 > pow((k * deltay), 2) + pow(double(deltay), 2))
+ deltay = deltay + 20;
+ deltax = ceil(deltay / -k2);
+ rangeL = pointLeft.x - deltax;
+ rangeR = pointRight.x + deltax;
+ rangeB = pointRight.y - deltay;
+ rangeT = pointLeft.y + deltay + 4;
+
+ cutImgLTx = pointLeft.x - rangeL;
+ cutImgLTy = pointLeft.y - rangeB;
+ cutImgRBx = pointRight.x - rangeL;
+ cutImgRBy = pointRight.y - rangeB;
+ }
+ if (rangeL - rangeR == 0)
+ return RER_H;
+ if (rangeL < 0)
+ rangeL = 0;
+ if (rangeB < 0)
+ rangeB = 0;
+ if (rangeT > img.rows)
+ rangeT = img.rows;
+ if (rangeR > img.cols)
+ rangeR = img.cols;
+ //qDebug() << QString("rangeB:%1,rangeT:%2; rangeL:%3,rangeR:%4").arg(rangeB).arg(rangeT).arg(rangeL).arg(rangeR);
+
+ cv::Mat miniEdgeX = img(cv::Range(rangeB, rangeT), cv::Range(rangeL, rangeR));
+ cv::Mat evalueImg;
+ miniEdgeX.copyTo(evalueImg);
+
+ pt1 = cv::Point(cutImgLTx, cutImgLTy);
+ pt2 = cv::Point(cutImgRBx, cutImgRBy);
+ cv::Mat temp(evalueImg.rows, evalueImg.cols, evalueImg.type(), cv::Scalar(0));
+ mask = temp;
+
+ int rows = evalueImg.rows;
+ int cols = evalueImg.cols;
+ for (int row = 0; row < rows; row++)
+ {
+ for (int col = 0; col < cols; col++)
+ {
+ //当前采样点
+ cv::Point samplingPoint(col, row);
+ //用于计算当前点到刃边距离
+ arraysampling_to_edge = { samplingPoint, pt1, pt2 };
+ //求当前点到刃边的垂足，获取垂足的col，用于比较确认dis的正负
+ cv::Point2d samplingPointFoot = calculateFootPoint(sampling_to_edge);
+ if (samplingPointFoot.x < pt1.x || samplingPointFoot.x > pt2.x)
+ continue;
+ //采样点距离刃边的距离
+ int dis = pointDistance(samplingPointFoot, samplingPoint);
+ if (dis > 21)
+ continue;
+ mask.at(row, col) = evalueImg.at(row, col);
+
+ if (samplingPointFoot.y > row)
+ dis = 0 - dis;
+
+ // int value = evalueImg.at(row, col);
+ //存入vecESF中
+ //遍历vecESF，寻找相同距离的值
+ // if (qmapESF.contains(dis))
+ // {
+ vectorpixValues = qmapESF.value(dis);
+ pixValues.push_back(evalueImg.at(row, col));
+ qmapESF.insert(dis, pixValues);
+ // }
+ // else
+ // {
+ // vectorpixValues;
+ // pixValues.push_back(evalueImg.at(row, col));
+ // qmapESF.insert(dis, pixValues);
+ // }
+ }
+ }
+
+ }
+ //cv::line(mask, pt1, pt2, cv::Scalar(0, 0, 0));
+ //vecESF中存的值求平均得ESF
+ QMapESF;
+ QMap>::const_iterator it = qmapESF.constBegin();
+ while (it != qmapESF.constEnd())
+ {
+ if (it.value().size() == 1)
+ ESF.insert(it.key(), it.value().at(0));
+ else
+ {
+ vectorsamdisValue = it.value();
+ int sum = 0;
+ int mean = 0;
+ for (int valuesize = 0; valuesize < samdisValue.size(); valuesize++)
+ sum = sum + samdisValue[valuesize];
+ int sameSize = samdisValue.size();
+ mean = sum / samdisValue.size();
+ ESF.insert(it.key(), mean);
+ }
+ ++it;
+ }
+ QMap::const_iterator it2 = ESF.constBegin();
+ int min = it2.value();
+ int max = it2.value();
+ while (it2 != ESF.constEnd())
+ {
+ if (min > it2.value())
+ min = it2.value();
+ if (max < it2.value())
+ max = it2.value();
+ ++it2;
+ }
+ int edgeValuePdot5 = ESF.value(2);
+ int edgeValueNdot5 = ESF.value(-2);
+ double ER_p05 = (double(edgeValuePdot5) - min) / (double(max) - min);
+ double ER_n05 = (double(edgeValueNdot5) - min) / (double(max) - min);
+ double resultRER = ER_p05 - ER_n05;
+ int edgeValue125 = ESF.value(5);
+ if (resultRER < 0)
+ {
+ edgeValue125 = ESF.value(-5);
+ resultRER = abs(resultRER);
+ }
+ double resultH = (double(edgeValue125) - min) / (double(max) - min);
+
+ RER_H.first = resultRER;
+ RER_H.second = resultH;
+
+ if (resultRER > 0.55)
+ {
+ cv::line(mask, pt1, pt2, cv::Scalar(0, 0, 0));
+ qDebug() << "resultRER:" << resultRER << "resultH:" << resultH;
+ }
+
+
+ return RER_H;
+}
+
+cv::Point2d NIIRS::calculateFootPoint(array& ps)
+{
+ cv::Point2d p(0.0, 0.0); // 垂足
+ if (ps.at(1).x == ps.at(2).x) // 线与x轴垂直
+ {
+ p.x = ps.at(1).x;
+ p.y = ps.at(0).y;
+ }
+ else if (ps.at(1).y == ps.at(2).y) // 线与y轴垂直
+ {
+ p.x = ps.at(0).x;
+ p.y = ps.at(1).y;
+ }
+ else // 线与x轴，y轴都不垂直
+ {
+ double a1 = double(ps.at(1).y) - double(ps.at(2).y);
+ double b1 = double(ps.at(2).x) - double(ps.at(1).x);
+ double c1 = double(ps.at(1).x) * double(ps.at(2).y) - double(ps.at(2).x) * float(ps.at(1).y);
+ //pFoot.x = (B * B * pA.x - A * B * pA.y - A * C) / (A * A + B * B);
+ //pFoot.y = (A * A * pA.y - A * B * pA.x - B * C) / (A * A + B * B);
+ double footx = (b1 * b1 * double(ps.at(0).x) - a1 * b1 * double(ps.at(0).y) - a1 * c1) / (a1 * a1 + b1 * b1);
+ double footy = (a1 * a1 * double(ps.at(0).y) - a1 * b1 * double(ps.at(0).x) - b1 * c1) / (a1 * a1 + b1 * b1);
+ //p.x = (b1 * b1 * ps.at(0).x - a1 * b1 * ps.at(0).y - a1 * c1) / (a1 * a1 + b1 * b1);
+ //p.y = (a1 * a1 * ps.at(0).y - a1 * b1 * ps.at(0).x - b1 * c1) / (a1 * a1 + b1 * b1);
+ p.x = footx;
+ p.y = footy;
+ }
+ return p;
+}
+
+QPair NIIRS::calculateERx(cv::Mat img)
+{
+ QPairER_H;
+
+ cv::Mat imgGray, imgGrayX4;
+ if (img.channels() == 3)
+ cv::cvtColor(img, imgGray, CV_BGR2GRAY);
+ else
+ imgGray = img;
+ //重采样为原来4倍，用于1/4像元大小的超采样
+ cv::pyrUp(imgGray, imgGrayX4);
+ cv::pyrUp(imgGrayX4, imgGrayX4);
+
+ ////高斯滤波
+ //cv::Mat imgGaussBlur;
+ //cv::GaussianBlur(imgGray, imgGaussBlur, cv::Size(3, 3), 3, 3);
+ //cv::imshow("imgGaussBlur", imgGaussBlur);
+
+ cv::Mat imgThres;
+ cv::threshold(imgGray, imgThres, 127, 255, cv::THRESH_BINARY);
+ //cv::imwrite("D:\\A_testdata\\niirs_test\\imgThresx.png", imgThres);
+ //cv::imshow("imgThres", imgThres);
+ //cv::waitKey(1);
+
+ //canny边缘提取
+ cv::Mat imgCanny;
+ cv::Canny(imgThres, imgCanny, 100, 200);
+ //cv::imshow("imgCanny", imgCanny);
+ //cv::waitKey(1);
+
+ vectorR_lines;
+ cv::HoughLinesP(imgCanny, R_lines, 1, CV_PI / 180, 10, 25, 100);
+ //获取所有直线的端点坐标，寻找两两之间距离最大的作为直线端点pt1、pt2
+ //for (size_t i = 0; i < R_lines.size(); i++)
+ //{
+ // cv::Vec4i l = R_lines[i];
+ // cv::line(imgGrayX4, cv::Point(l[0] * 4, l[1] * 4), cv::Point(l[2] * 4, l[3] * 4), cv::Scalar(255 / (i+1)), 3);
+ //}
+ //cv::imshow("imgGrayX4", imgGrayX4);
+ //return 0.0;
+
+ vectortop_points;
+ vectorbot_points;
+ if (R_lines.size() == 0)
+ {
+ qDebug() << "find no line in edge x";
+ ui.lineEdit_INFO->setText(QString::fromLocal8Bit("无法从输入数据中提取刃边"));
+ return ER_H;
+ }
+ for (size_t z = 0; z < R_lines.size(); z++)
+ {
+ cv::Vec4i p = R_lines[z];
+ //cout << "line x0: " << cv::Point(p[0], p[1]) <<", x1: "<< cv::Point(p[2], p[3]) << std::endl;
+ top_points.push_back(cv::Point(p[0], p[1]));
+ bot_points.push_back(cv::Point(p[2], p[3]));
+ }
+ //用于存放点之间距离
+ cv::Mat img_dis(top_points.size(), bot_points.size(), CV_32FC1, cv::Scalar(0));
+ for (int row = 0; row < top_points.size(); row++)
+ {
+ cv::Vec* data_Ptr = img_dis.ptr>(row);
+ for (int col = 0; col < bot_points.size(); col++)
+ {
+ if (row == col)
+ data_Ptr[col] = 0.0;
+ else
+ data_Ptr[col] = pointDistance(top_points[row], bot_points[col]);
+ }
+ }
+ double minv, maxv;
+ cv::Point pt_min, pt_max;
+ cv::minMaxLoc(img_dis, &minv, &maxv, &pt_min, &pt_max);
+ //std::cout << "img_dis = " << std::endl << img_dis << std::endl;
+ //std::cout << "maxv = " << maxv << std::endl;
+ //std::cout << "idx_max = " << pt_max << std::endl;
+
+ //刃边端点
+ cv::Point pt1 = cv::Point(top_points[pt_max.y].x * 4, top_points[pt_max.y].y * 4);
+ cv::Point pt2 = cv::Point(bot_points[pt_max.x].x * 4, bot_points[pt_max.x].y * 4);
+ //std::cout << "Edge End Point: " << pt1 << ", " << pt2 << std::endl;
+
+ cv::Mat masksrc(imgCanny.rows * 4, imgCanny.cols * 4, imgCanny.type(), cv::Scalar(0));
+ int rows = masksrc.rows;
+ int cols = masksrc.cols;
+ //cout << "row:" << rows << ", col:" << cols << std::endl;
+
+ //cv::line(masksrc, pt1, pt2, cv::Scalar(255), 1, 8);
+ //cv::imshow("xxxedge", masksrc);
+ //cv::waitKey(1);
+ //cv::imwrite("D:\\A_testdata\\niirs_test\\imgLinex.png", masksrc);
+
+ //pair中:vector 存放所有相同距离时像元值，之后求平均
+ //pair中:int 表示离刃边距离
+ vector>>vecESF;
+
+ for (int row = 0; row < rows; row++)
+ {
+ for (int col = 0; col < cols; col++)
+ {
+ //当前采样点
+ cv::Point samplingPoint(col, row);
+ //用于计算当前点到刃边距离
+ arraysampling_to_edge = { samplingPoint, pt1, pt2 };
+ //采样点距离刃边的距离
+ int dis = 0;
+ //cout << "mask at row:" << row << ",col:" << col;
+ //求当前点到刃边的垂足，获取垂足的col，用于比较确认dis的正负
+ cv::Point samplingPointFoot = calculate_foot_point(sampling_to_edge);
+ //cout << ",edge foot at row:" << samplingPointFoot.y << ",col:" << samplingPointFoot.x;
+ if (samplingPointFoot.x < col)
+ {
+ //dis = calculate_distance(sampling_to_edge);
+ dis = pointDistance(samplingPointFoot, samplingPoint);
+ }
+ else
+ {
+ //dis = 0 - calculate_distance(sampling_to_edge);
+ dis = 0 - pointDistance(samplingPointFoot, samplingPoint);
+ }
+ int value = imgGrayX4.at(row, col);
+ //cout << ",distance to edge:" << dis << ",value:"
+ // << value << std::endl;
+ pair>pix_dis_value;
+ //存入vecESF中
+ if (vecESF.size() == 0)//vecESF大小为零，直接放入
+ {
+ pix_dis_value.first = dis;
+ pix_dis_value.second.push_back(imgGrayX4.at(row, col));
+ vecESF.push_back(pix_dis_value);
+ }
+ else//vecESF大小不为零，查找有无距离相同的点存入
+ {
+ //遍历vecESF，寻找相同距离的值
+ for (int vecsize = 0; vecsize < vecESF.size(); vecsize++)
+ {
+ if (vecESF[vecsize].first == dis)//有，存入
+ {
+ vecESF[vecsize].second.push_back(imgGrayX4.at(row, col));
+ break;
+ }
+ if (vecsize == vecESF.size() - 1)//遍历至最后找不到，新建存入
+ {
+ pix_dis_value.first = dis;
+ pix_dis_value.second.push_back(imgGrayX4.at(row, col));
+ vecESF.push_back(pix_dis_value);
+ break;//上一行vecESF.size+1,break防止无限循环
+ }
+ }
+ }
+ }
+ }
+ //cout << "vecESF.size():" << vecESF.size() << std::endl;
+ //vecESF中存的值求平均得ESF
+ vector>ESF;
+ //求ESF
+ for (int vecsize = 0; vecsize < vecESF.size(); vecsize++)
+ {
+ if (vecESF[vecsize].second.size() == 1)
+ {
+ pairsamedisESF(vecESF[vecsize].first, vecESF[vecsize].second[0]);
+ ESF.push_back(samedisESF);
+ }
+ else
+ {
+ vectorsamdisValue = vecESF[vecsize].second;
+ int sum = 0;
+ int mean = 0;
+ for (int valuesize = 0; valuesize < samdisValue.size(); valuesize++)
+ {
+ sum = sum + samdisValue[valuesize];
+ }
+ int sameSize = samdisValue.size();
+ mean = sum / samdisValue.size();
+ pairsamedisESF(vecESF[vecsize].first, mean);
+ ESF.push_back(samedisESF);
+ }
+ }
+ //遍历ESF，距离大于0，求平均，距离小于0，求平均
+ pairDN_max_min;
+ vectorpositive_side, negative_side;
+ // 距离刃边+0.5 距离刃边-0.5 距离刃边+0.5
+ double edgeValue_p05, edgeValue_n05, edgeValue_p1dot25;
+ for (int i = 0; i < ESF.size(); i++)
+ {
+ if (ESF[i].first > 0)
+ positive_side.push_back(ESF[i].second);
+ else if (ESF[i].first < 0)
+ negative_side.push_back(ESF[i].second);
+ if (ESF[i].first == 2)
+ edgeValue_p05 = ESF[i].second;
+ if (ESF[i].first == -2)
+ edgeValue_n05 = ESF[i].second;
+ if (ESF[i].first == 5 && ESF[i].second > 0)
+ edgeValue_p1dot25 = abs(ESF[i].second);
+ }
+ double sumpos = accumulate(begin(positive_side), end(positive_side), 0.0);
+ double sumneg = accumulate(begin(negative_side), end(negative_side), 0.0);
+ double mean_positive = sumpos / positive_side.size();
+ double mean_negative = sumneg / negative_side.size();
+
+ DN_max_min.first = mean_positive > mean_negative ? mean_positive : mean_negative;
+ DN_max_min.second = mean_positive < mean_negative ? mean_positive : mean_negative;
+ //cout << "lager side:" << DN_max_min.first << ",small side:" << DN_max_min.second << std::endl;
+ //edgeValue_n05
+ //edgeValue_p05
+
+ double ER_p05 = (edgeValue_p05 - DN_max_min.second) / (DN_max_min.first - DN_max_min.second);
+ double ER_n05 = (edgeValue_n05 - DN_max_min.second) / (DN_max_min.first - DN_max_min.second);
+ double Hx = double(edgeValue_p1dot25 - DN_max_min.second) / double(DN_max_min.first - DN_max_min.second);
+ double ER = ER_p05 - ER_n05;
+ ER_H.first = ER;
+ ER_H.second = Hx;
+ return ER_H;
+
+ /*-----------------------------------------------------------------------------------------------------------*/
+
+ ////ESF求斜率求LSF
+ ////在ESF最后扩展一个用于求斜率
+ //pairtempend(ESF[ESF.size()].first + 1, ESF[ESF.size()].second);
+ //ESF.push_back(tempend);
+ //for (int i = 0; i < ESF.size() - 1; i++)//遍历至倒数第二位
+ //{
+ // pairtemp;
+ // temp.first = ESF[i].first;
+ // temp.second = (ESF[i].second - ESF[i + 1].second) / (ESF[i].first - ESF[i + 1].first);
+ // LSF.push_back(temp);
+ //}
+ //vectorvec_LSF[2];
+ //for (int vecsize = 0; vecsize < LSF.size(); vecsize++)
+ //{
+ // vec_LSF[0].push_back(LSF[vecsize].first);
+ // vec_LSF[1].push_back(LSF[vecsize].second);
+ //}
+ ////一行
+ ////cv::Mat mat_LSF = cv::Mat::zeros(1, vec_LSF[1].size(), CV_32F);// = cv::Mat(vec_LSF)
+ ////for (int i = 0; i < vec_LSF[1].size(); i++)
+ //// mat_LSF.at(0, i) = float(vec_LSF[1][i]);
+ ////cv::imwrite("mat_LSF.png", mat_LSF);
+ ////std::cout <<"mat_LSF channels: " << mat_LSF.channels() << ", size(width x height): " << mat_LSF.size() << std::endl;
+ //cv::Mat fft_real, fft_imag, fft_complex;//复数实部，复数虚部，完整复数
+ //fft_real = cv::Mat::zeros(1, vec_LSF[1].size(), CV_32F);
+ //fft_imag = cv::Mat::zeros(1, vec_LSF[1].size(), CV_32F);
+ //for (int i = 0; i < vec_LSF[1].size(); i++)
+ //{
+ // fft_real.at(0, i) = float(vec_LSF[1][i]);
+ // fft_imag.at(0, i) = 0;
+ //}
+ //std::cout << "1.fft_real channels: " << fft_real.channels() << ", size(width x height): " << fft_real.size() << std::endl;
+ //std::cout << "1.fft_imag channels: " << fft_imag.channels() << ", size(width x height): " << fft_imag.size() << std::endl;
+ //vectorvec_complex;
+ //vec_complex.push_back(fft_real);
+ //vec_complex.push_back(fft_imag);
+ //cv::merge(vec_complex, fft_complex);
+ //cv::dft(fft_complex, fft_complex);
+ //cv::Mat twoChannels[2];
+ ////重新拆分出实部虚部
+ //cv::split(fft_complex, twoChannels);
+ ////std::cout << "----2.fft_real----" << std::endl << twoChannels[0] << std::endl
+ //// << "channels: " << twoChannels[0].channels() << ", size(width x height): " << twoChannels[0].size() << std::endl;
+ ////std::cout << "----2.fft_imag----" << std::endl << twoChannels[1] << std::endl
+ //// << "channels: " << twoChannels[1].channels() << ", size(width x height): " << twoChannels[1].size() << std::endl;
+ //cv::Mat matMTF = cv::Mat::zeros(1, vec_LSF[1].size(), CV_32F);
+ //for (int i = 0; i < vec_LSF[1].size(); i++)
+ //{
+ // matMTF.at(0, i)
+ // = sqrt(float(twoChannels[0].at(0, i)) * float(twoChannels[0].at(0, i))
+ // + float(twoChannels[1].at(0, i)) * float(twoChannels[1].at(0, i)));
+ //}
+ //std::cout << "matMTF" << std::endl << matMTF << std::endl
+ // << "channels: " << matMTF.channels() << ", size(width x height): " << matMTF.size() << std::endl;
+ ////暂时取前1/2作为mtf，xy轴归一化，做积分算ER
+ //cout << "--------------------------------------------" << std::endl;
+ //return 0.0;
+}
+
+QPair NIIRS::calculateERy(cv::Mat img)
+{
+ QPairER_H;
+
+ cv::Mat imgGray, imgGrayX4;
+ if (img.channels() == 3)
+ cv::cvtColor(img, imgGray, CV_BGR2GRAY);
+ else
+ imgGray = img;
+ //重采样为原来4倍，用于1/4像元大小的超采样
+ cv::pyrUp(imgGray, imgGrayX4);
+ cv::pyrUp(imgGrayX4, imgGrayX4);
+ //cv::imshow("yyyimgGrayX4", imgGrayX4);
+
+ ////高斯滤波
+ //cv::Mat imgGaussBlur;
+ //cv::GaussianBlur(imgGray, imgGaussBlur, cv::Size(3, 3), 3, 3);
+ cv::Mat imgThres;
+ cv::threshold(imgGray, imgThres, 127, 255, cv::THRESH_BINARY);
+ //cv::imwrite("D:\\A_testdata\\niirs_test\\imgGrayy.png", imgThres);
+ //cv::imshow("imgThresy", imgThres);
+ //cv::waitKey(1);
+
+ //canny边缘提取
+ cv::Mat imgCanny;
+ cv::Canny(imgThres, imgCanny, 100, 200);
+ //cv::imshow("imgCannyy", imgCanny);
+ //cv::waitKey(1);
+
+ vectorR_lines;
+ cv::HoughLinesP(imgCanny, R_lines, 1, CV_PI / 180, 10, 25, 100);
+ //获取所有直线的端点坐标，寻找两两之间距离最大的作为直线端点pt1、pt2
+ vectortop_points;
+ vectorbot_points;
+ if (R_lines.size() == 0)
+ {
+ qDebug() << "find no line in edge y";
+ ui.lineEdit_INFO->setText(QString::fromLocal8Bit("无法从输入数据中提取刃边"));
+ return ER_H;
+ }
+ for (size_t z = 0; z < R_lines.size(); z++)
+ {
+ cv::Vec4i p = R_lines[z];
+ //cout << "line x0: " << cv::Point(p[0], p[1]) << ", x1: " << cv::Point(p[2], p[3]) << std::endl;
+ top_points.push_back(cv::Point(p[0], p[1]));
+ bot_points.push_back(cv::Point(p[2], p[3]));
+ }
+ //用于存放点之间距离
+ cv::Mat img_dis(top_points.size(), bot_points.size(), CV_32FC1, cv::Scalar(0));
+ for (int row = 0; row < top_points.size(); row++)
+ {
+ cv::Vec* data_Ptr = img_dis.ptr>(row);
+ for (int col = 0; col < bot_points.size(); col++)
+ {
+ if (row == col)
+ data_Ptr[col] = 0.0;
+ else
+ data_Ptr[col] = pointDistance(top_points[row], bot_points[col]);
+ }
+ }
+ double minv, maxv;
+ cv::Point pt_min, pt_max;
+ cv::minMaxLoc(img_dis, &minv, &maxv, &pt_min, &pt_max);
+ //std::cout << "img_dis = " << std::endl << img_dis << std::endl;
+ //std::cout << "maxv = " << maxv << std::endl;
+ //std::cout << "idx_max = " << pt_max << std::endl;
+
+ //刃边端点
+ cv::Point pt1 = cv::Point(top_points[pt_max.y].x * 4, top_points[pt_max.y].y * 4);
+ cv::Point pt2 = cv::Point(bot_points[pt_max.x].x * 4, bot_points[pt_max.x].y * 4);
+ //std::cout << "Edge End Point: " << pt1 << ", " << pt2 << std::endl;
+
+ cv::Mat masksrc(imgCanny.rows * 4, imgCanny.cols * 4, imgCanny.type(), cv::Scalar(0));
+ int rows = masksrc.rows;
+ int cols = masksrc.cols;
+ //cout << "row:" << rows << ", col:" << cols << std::endl;
+
+ cv::line(masksrc, pt1, pt2, cv::Scalar(255), 1, 8);
+ //cv::imshow("yyyedge", masksrc);
+ //cv::imwrite("D:\\A_testdata\\niirs_test\\imgLiney.png", masksrc);
+
+ //pair中:vector 存放所有相同距离时像元值，之后求平均
+ //pair中:int 表示离刃边距离
+ vector>>vecESF;
+
+ for (int row = 0; row < rows; row++)
+ {
+ for (int col = 0; col < cols; col++)
+ {
+ //当前采样点
+ cv::Point samplingPoint(col, row);
+ //用于计算当前点到刃边距离
+ arraysampling_to_edge = { samplingPoint, pt1, pt2 };
+ //采样点距离刃边的距离
+ int dis = 0;
+ //cout << "mask at row:" << row << ",col:" << col;
+ //求当前点到刃边的垂足，获取垂足的col，用于比较确认dis的正负
+ cv::Point samplingPointFoot = calculate_foot_point(sampling_to_edge);
+ //cout << ",edge foot at row:" << samplingPointFoot.y << ",col:" << samplingPointFoot.x;
+ if (samplingPointFoot.x < col)
+ {
+ //dis = calculate_distance(sampling_to_edge);
+ dis = pointDistance(samplingPointFoot, samplingPoint);
+ }
+ else
+ {
+ //float a1 = float(sampling_to_edge.at(1).y) - float(sampling_to_edge.at(2).y);
+ //float b1 = float(sampling_to_edge.at(2).x) - float(sampling_to_edge.at(1).x);
+ //float c1 = float(sampling_to_edge.at(1).x) * float(sampling_to_edge.at(2).y)
+ // - float(sampling_to_edge.at(2).x) * float(sampling_to_edge.at(1).y);
+ //dis = abs(a1 * sampling_to_edge.at(0).x + b1 * sampling_to_edge.at(0).y + c1)
+ // / sqrt(a1 * a1 + b1 * b1);
+ dis = 0 - pointDistance(samplingPointFoot, samplingPoint);
+ }
+ int value = imgGrayX4.at(row, col);
+ //cout << ",distance to edge:" << dis << ",value:"
+ // << value << std::endl;
+ pair>pix_dis_value;
+ //存入vecESF中
+ if (vecESF.size() == 0)//vecESF大小为零，直接放入
+ {
+ pix_dis_value.first = dis;
+ pix_dis_value.second.push_back(imgGrayX4.at(row, col));
+ vecESF.push_back(pix_dis_value);
+ }
+ else//vecESF大小不为零，查找有无距离相同的点存入
+ {
+ //遍历vecESF，寻找相同距离的值
+ for (int vecsize = 0; vecsize < vecESF.size(); vecsize++)
+ {
+ if (vecESF[vecsize].first == dis)//有，存入
+ {
+ vecESF[vecsize].second.push_back(imgGrayX4.at(row, col));
+ break;
+ }
+ if (vecsize == vecESF.size() - 1)//遍历至最后找不到，新建存入
+ {
+ pix_dis_value.first = dis;
+ pix_dis_value.second.push_back(imgGrayX4.at(row, col));
+ vecESF.push_back(pix_dis_value);
+ break;//上一行vecESF.size+1,break防止无限循环
+ }
+ }
+ }
+ }
+ }
+ //cout << "vecESF.size():" << vecESF.size() << std::endl;
+ //vecESF中存的值求平均得ESF
+ vector>ESF;
+ //ESF求斜率得LSF
+ vector>LSF;
+ //求ESF
+ for (int vecsize = 0; vecsize < vecESF.size(); vecsize++)
+ {
+ if (vecESF[vecsize].second.size() == 1)
+ {
+ pairsamedisESF(vecESF[vecsize].first, vecESF[vecsize].second[0]);
+ ESF.push_back(samedisESF);
+ }
+ else
+ {
+ vectorsamdisValue = vecESF[vecsize].second;
+ int sum = 0;
+ int mean = 0;
+ for (int valuesize = 0; valuesize < samdisValue.size(); valuesize++)
+ {
+ sum = sum + samdisValue[valuesize];
+ }
+ int sameSize = samdisValue.size();
+ mean = sum / samdisValue.size();
+ pairsamedisESF(vecESF[vecsize].first, mean);
+ ESF.push_back(samedisESF);
+ }
+ }
+ //遍历ESF，距离大于0，求平均，距离小于0，求平均
+ pairDN_max_min;
+ vectorpositive_side, negative_side;
+ // 距离刃边+0.5 距离刃边-0.5
+ double edgeValue_p05, edgeValue_n05, edgeValue_p1dot25;
+ for (int i = 0; i < ESF.size(); i++)
+ {
+ //cout << ESF[vecSize].first << ",," << ESF[vecSize].second << std::endl;
+ if (ESF[i].first > 0)
+ positive_side.push_back(ESF[i].second);
+ else if (ESF[i].first < 0)
+ negative_side.push_back(ESF[i].second);
+ if (ESF[i].first == 2)
+ edgeValue_p05 = ESF[i].second;
+ if (ESF[i].first == -2)
+ edgeValue_n05 = ESF[i].second;
+ if (ESF[i].first == 5 && ESF[i].second > 0)
+ edgeValue_p1dot25 = abs(ESF[i].second);
+ }
+ double sumpos = accumulate(begin(positive_side), end(positive_side), 0.0);
+ double sumneg = accumulate(begin(negative_side), end(negative_side), 0.0);
+ double mean_positive = sumpos / positive_side.size();
+ double mean_negative = sumneg / negative_side.size();
+
+ DN_max_min.first = mean_positive > mean_negative ? mean_positive : mean_negative;
+ DN_max_min.second = mean_positive < mean_negative ? mean_positive : mean_negative;
+ //cout << "lager side:" << DN_max_min.first << ",small side:" << DN_max_min.second << std::endl;
+ //edgeValue_n05
+ //edgeValue_p05
+ double ER_p05 = (edgeValue_p05 - DN_max_min.second) / (DN_max_min.first - DN_max_min.second);
+ double ER_n05 = (edgeValue_n05 - DN_max_min.second) / (DN_max_min.first - DN_max_min.second);
+ double Hx = double(edgeValue_p1dot25 - DN_max_min.second) / double(DN_max_min.first - DN_max_min.second);
+ double ER = ER_p05 - ER_n05;
+ ER_H.first = ER;
+ ER_H.second = Hx;
+ return ER_H;
+}
+
+float NIIRS::pointDistance(cv::Point pt1, cv::Point pt2)
+{
+ //float dis = sqrtf(((float)pt1.x - (float)pt2.x) * ((float)pt1.x - (float)pt2.x)
+ // + ((float)pt1.y - (float)pt2.y) * ((float)pt1.y - (float)pt2.y));
+ float dis = sqrt((float(pt1.x) - float(pt2.x)) * (float(pt1.x) - float(pt2.x))
+ + (float(pt1.y) - float(pt2.y)) * (float(pt1.y) - float(pt2.y)));
+ //std::cout << "pt1:" << pt1 << "pt2" << pt2 << "pointDistance:" << dis << std::endl;
+ return dis;
+}
+
+//int NIIRS::calculate_distance(array& ps)
+//{
+// int d = 0;// 距离
+// int a1 = -(ps.at(2).y - ps.at(1).y);
+// int b1 = ps.at(2).x - ps.at(1).x;
+// int c1 = (ps.at(2).y - ps.at(1).y) * ps.at(1).x - (ps.at(2).x - ps.at(1).x) * ps.at(1).y;
+// d = abs(a1 * ps.at(0).x + b1 * ps.at(0).y + c1) / sqrt(a1 * a1 + b1 * b1);
+// //allnum++;
+// //if (allnum % 100 == 0)
+// // cout << "calculate_distance:" << d << ",at:" << allnum << std::endl;
+// return d;
+//}
+
+cv::Point NIIRS::calculate_foot_point(array& ps)
+{
+ cv::Point p(0, 0); // 垂足
+ if (ps.at(1).x == ps.at(2).x) // 线与x轴垂直
+ {
+ p.x = ps.at(1).x;
+ p.y = ps.at(0).y;
+ }
+ else if (ps.at(1).y == ps.at(2).y) // 线与y轴垂直
+ {
+ p.x = ps.at(0).x;
+ p.y = ps.at(1).y;
+ }
+ else // 线与x轴，y轴都不垂直
+ {
+ float a1 = float(ps.at(1).y) - float(ps.at(2).y);
+ float b1 = float(ps.at(2).x) - float(ps.at(1).x);
+ float c1 = float(ps.at(1).x) * float(ps.at(2).y) - float(ps.at(2).x) * float(ps.at(1).y);
+ //pFoot.x = (B * B * pA.x - A * B * pA.y - A * C) / (A * A + B * B);
+ //pFoot.y = (A * A * pA.y - A * B * pA.x - B * C) / (A * A + B * B);
+ float footx = (b1 * b1 * float(ps.at(0).x) - a1 * b1 * float(ps.at(0).y) - a1 * c1) / (a1 * a1 + b1 * b1);
+ float footy = (a1 * a1 * float(ps.at(0).y) - a1 * b1 * float(ps.at(0).x) - b1 * c1) / (a1 * a1 + b1 * b1);
+ //p.x = (b1 * b1 * ps.at(0).x - a1 * b1 * ps.at(0).y - a1 * c1) / (a1 * a1 + b1 * b1);
+ //p.y = (a1 * a1 * ps.at(0).y - a1 * b1 * ps.at(0).x - b1 * c1) / (a1 * a1 + b1 * b1);
+ p.x = footx;
+ p.y = footy;
+ }
+ return p;
+}
+
+//cv::Point NIIRS::line_intersection(array& line1, array& line2)
+//{
+// cv::Point inter;
+//
+// double ka, kb;
+// // line1 y2 line1 y1 line1 x2 line1 x1
+// //ka = (double)(LineA[3] - LineA[1]) / (double)(LineA[2] - LineA[0]); //求出LineA斜率
+// // line2 y2 line2 y1 line2 x2 line2 x1
+// //kb = (double)(LineB[3] - LineB[1]) / (double)(LineB[2] - LineB[0]); //求出LineB斜率
+//
+// // line1 x1 line1 y1 line2 x1 line2 y1
+// //inter.x = (ka * LineA[0] - LineA[1] - kb * LineB[0] + LineB[1]) / (ka - kb);
+// // line1 x1 line2 x1 line2 y1 line1 y1
+// //inter.y = (ka * kb * (LineA[0] - LineB[0]) + ka * LineB[1] - kb * LineA[1]) / (ka - kb);
+// ka = (double(line1[1].y) - double(line1[0].y)) / (double(line1[1].x) - double(line1[0].x)); //求出LineA斜率
+// kb = (double(line2[1].y) - double(line2[0].y)) / (double(line2[1].x) - double(line2[0].x)); //求出LineB斜率
+//
+// inter.x = (ka * line1[0].x - line1[0].y - kb * line2[0].x + line2[0].y) / (ka - kb);
+// inter.y = (ka * kb * (double(line1[0].x) - double(line2[0].x)) + ka * line2[0].y - kb * line1[0].y) / (ka - kb);
+//
+// return inter;
+//}
+
+//cv::Point NIIRS::point_k_edge_intersection(cv::Point pt, double k, array& line)
+//{
+// cv::Point inter;
+// // 根据点斜式求直线
+// // 根据直线是水平或垂直，求两线交点
+// //cout << "point_k_edge_intersection point: " << pt << std::endl;
+// if (line[0].x == line[1].x)//x相等，相当于x已知，求y = k * (x - x1) + y1
+// {
+// //cout << "line[0].x == line[1].x: " << line[0] << "," << line[1] << std::endl;
+// float interpoint = k * (double(line[0].x) - double(pt.x)) + double(pt.y);
+// inter.x = line[0].x;
+// inter.y = round(interpoint);
+// }
+// else if (line[0].y == line[1].y)//y相等，相当于y已知，求x = (y - y1) / k + x1
+// {
+// //cout << "line[0].y == line[1].y: " << line[0] << "," << line[1] << std::endl;
+// float interpoint = (double(line[0].y) - double(pt.y)) / k + double(pt.x);
+// inter.x = round(interpoint);
+// inter.y = line[0].y;
+// }
+// //cout << "intersection: " << inter << std::endl;
+//
+// return inter;
+//}
+
+//基于sobel算子求梯度、陡度
+double NIIRS::gradient_Sobel(cv::Mat img)
+{
+ //assert(img.empty());
+ cv::Mat gray_img, sobel_x, abs_x, sobel_y, abs_y, sobel_mean;
+ if (img.channels() == 3)
+ cv::cvtColor(img, gray_img, CV_BGR2GRAY);
+ else
+ gray_img = img;
+ //分别计算x、y方向梯度
+ cv::Sobel(gray_img, sobel_x, CV_32FC1, 1, 0);
+ cv::Sobel(gray_img, sobel_y, CV_32FC1, 0, 1);
+
+ //cv::multiply(sobel_x, sobel_x, sobel_x);
+ //cv::multiply(sobel_y, sobel_y, sobel_y);
+ //cv::Mat sqrt_mat = sobel_x + sobel_y;
+ //cv::sqrt(sqrt_mat, sobel_mean);
+ cv::convertScaleAbs(sobel_x, abs_x);
+ cv::convertScaleAbs(sobel_y, abs_y);
+ cv::addWeighted(abs_x, 0.5, abs_y, 0.5, 0, sobel_mean);
+
+ double mean = cv::mean(sobel_mean)[0];
+ //cv::resize(sobel_mean, sobel_mean, cv::Size(600, 400));
+
+ return mean;
+}
+//基于Laplacian算子求二阶导，求清晰度
+double NIIRS::definition_Laplacian(cv::Mat img)
+{
+ cv::Mat gray_img, lap_image;
+ if (img.channels() == 3)
+ cv::cvtColor(img, gray_img, CV_BGR2GRAY);
+ else
+ gray_img = img;
+
+ cv::Laplacian(gray_img, lap_image, CV_32FC1);
+ lap_image = cv::abs(lap_image);
+
+ return cv::mean(lap_image)[0];
+}
+
+void NIIRS::get_GLCM_0deg(cv::Mat& input, cv::Mat& dst)
+{
+ cv::Mat src;
+ input.copyTo(src);
+ src.convertTo(src, CV_32S);
+ int height = src.rows;
+ int width = src.cols;
+ int max_gray_level = 0;
+ for (int j = 0; j < height; j++)//寻找像素灰度最大值
+ {
+ int* srcdata = src.ptr(j);
+ for (int i = 0; i < width; i++)
+ {
+ if (srcdata[i] > max_gray_level)
+ {
+ max_gray_level = srcdata[i];
+ }
+ }
+ }
+ max_gray_level++;//像素灰度最大值加1即为该矩阵所拥有的灰度级数
+ dst.create(max_gray_level, max_gray_level, CV_32SC1);
+ dst = cv::Scalar::all(0);
+ for (int i = 0; i < height; i++)
+ {
+ int* srcdata = src.ptr(i);
+ for (int j = 0; j < width - 1; j++)
+ {
+ int rows = srcdata[j];
+ int cols = srcdata[j + 1];
+ dst.ptr(rows)[cols]++;
+ }
+ }
+}
+
+void NIIRS::get_GLCM_90deg(cv::Mat& input, cv::Mat& dst)
+{
+ cv::Mat src;
+ input.copyTo(src);
+ src.convertTo(src, CV_32S);
+ int height = src.rows;
+ int width = src.cols;
+ int max_gray_level = 0;
+ for (int j = 0; j < height; j++)
+ {
+ int* srcdata = src.ptr(j);
+ for (int i = 0; i < width; i++)
+ {
+ if (srcdata[i] > max_gray_level)
+ {
+ max_gray_level = srcdata[i];
+ }
+ }
+ }
+ max_gray_level++;
+ dst.create(max_gray_level, max_gray_level, CV_32SC1);
+ dst = cv::Scalar::all(0);
+ for (int i = 0; i < height - 1; i++)
+ {
+ int* srcdata = src.ptr(i);
+ int* srcdata1 = src.ptr(i + 1);
+ for (int j = 0; j < width; j++)
+ {
+ int rows = srcdata[j];
+ int cols = srcdata1[j];
+ dst.ptr(rows)[cols]++;
+ }
+ }
+}
+
+void NIIRS::get_GLCM_45deg(cv::Mat& input, cv::Mat& dst)
+{
+ cv::Mat src;
+ input.copyTo(src);
+ src.convertTo(src, CV_32S);
+ int height = src.rows;
+ int width = src.cols;
+ int max_gray_level = 0;
+ for (int j = 0; j < height; j++)
+ {
+ int* srcdata = src.ptr(j);
+ for (int i = 0; i < width; i++)
+ {
+ if (srcdata[i] > max_gray_level)
+ {
+ max_gray_level = srcdata[i];
+ }
+ }
+ }
+ max_gray_level++;
+ dst.create(max_gray_level, max_gray_level, CV_32SC1);
+ dst = cv::Scalar::all(0);
+ for (int i = 0; i < height - 1; i++)
+ {
+ int* srcdata = src.ptr(i);
+ int* srcdata1 = src.ptr(i + 1);
+ for (int j = 0; j < width - 1; j++)
+ {
+ int rows = srcdata[j];
+ int cols = srcdata1[j + 1];
+ dst.ptr(rows)[cols]++;
+ }
+ }
+}
+
+void NIIRS::get_GLCM_135deg(cv::Mat& input, cv::Mat& dst)
+{
+ cv::Mat src;
+ input.copyTo(src);
+ src.convertTo(src, CV_32S);
+ int height = src.rows;
+ int width = src.cols;
+ int max_gray_level = 0;
+ for (int j = 0; j < height; j++)
+ {
+ int* srcdata = src.ptr(j);
+ for (int i = 0; i < width; i++)
+ {
+ if (srcdata[i] > max_gray_level)
+ {
+ max_gray_level = srcdata[i];
+ }
+ }
+ }
+ max_gray_level++;
+ dst.create(max_gray_level, max_gray_level, CV_32SC1);
+ dst = cv::Scalar::all(0);
+ for (int i = 0; i < height - 1; i++)
+ {
+ int* srcdata = src.ptr(i);
+ int* srcdata1 = src.ptr(i + 1);
+ for (int j = 1; j < width; j++)
+ {
+ int rows = srcdata[j];
+ int cols = srcdata1[j - 1];
+ dst.ptr(rows)[cols]++;
+ }
+ }
+}
+
+double NIIRS::calculateASM(cv::Mat& src)
+{
+ //cv::Mat src;
+ //cv::cvtColor(input, src, CV_BGR2GRAY);//转为灰度图
+
+ int height = src.rows;
+ int width = src.cols;
+ int total = 0;
+ for (int i = 0; i < height; i++)
+ {
+ int* srcdata = src.ptr(i);
+ for (int j = 0; j < width; j++)
+ {
+ total += srcdata[j];//求图像所有像素的灰度值的和
+ }
+ }
+
+ cv::Mat copy;
+ copy.create(height, width, CV_64FC1);
+ for (int i = 0; i < height; i++)
+ {
+ int* srcdata = src.ptr(i);
+ double* copydata = copy.ptr(i);
+ for (int j = 0; j < width; j++)
+ {
+ copydata[j] = (double)srcdata[j] / (double)total;//图像每一个像素的的值除以像素总和
+ }
+ }
+ //纹理特征的因子
+ double Asm = 0.0, Eng = 0.0, Con = 0.0;
+ for (int i = 0; i < height; i++)
+ {
+ double* srcdata = copy.ptr(i);
+ for (int j = 0; j < width; j++)
+ {
+ Asm += srcdata[j] * srcdata[j];//能量
+ if (srcdata[j] > 0)
+ Eng -= srcdata[j] * log(srcdata[j]);//熵//+0.0000001
+ Con += (double)(i - j) * (double)(i - j) * srcdata[j];//对比度
+ //Idm += srcdata[j] / (1 + (double)(i - j) * (double)(i - j));//逆差矩
+ }
+ }
+ //qDebug() << "--1. Angular Sec Moment" << Asm;
+ //qDebug() << "--2. Entropy-----------" << Eng;
+ //qDebug() << "--3. contranst radio---" << Con;
+
+ return Asm;
+}
+
+double NIIRS::calculatePIQEScore(cv::Mat& img)
+{
+ int blockSize = 16;
+ double activityThreshold = 0.1;
+ double blockImpairedThreshold = 0.1;
+ int windowSize = 6;
+ int nSegments = blockSize - windowSize + 1;
+ double distBlockScores = 0.0;
+ double NHSA = 0.0;
+ double c = 1.0;//常数
+
+ cv::Mat gray_img;
+ if (img.channels() == 3)
+ cv::cvtColor(img, gray_img, CV_BGR2GRAY);
+ else
+ gray_img = img;
+ //cv::Mat copy = gray_img;
+ //cv::resize(copy, copy, cv::Size(960, 540));
+ //cv::imshow("gray_img", copy);
+
+ gray_img.convertTo(gray_img, CV_64FC1);
+ int rows = gray_img.rows;
+ int columns = gray_img.cols;
+ int rowsPad = rows % blockSize;
+ int columnsPad = columns % blockSize;
+ //qDebug() << "rowsPad:" << rowsPad << "columnsPad:" << columnsPad;
+ bool isPadded = false;
+ if (rowsPad > 0 || columnsPad > 0)
+ {
+ if (rowsPad > 0)
+ rowsPad = blockSize - rowsPad;
+ if (columnsPad > 0)
+ columnsPad = blockSize - columnsPad;
+ isPadded = true;
+ //扩充图像，向右向下扩充，使图像行列能够被block整除，[rowsPad向下 columnsPad向右]
+ cv::copyMakeBorder(gray_img, gray_img, 0, rowsPad, 0, columnsPad, cv::BORDER_REPLICATE);
+ //rows = gray_img.rows;
+ //columns = gray_img.cols;
+ //rowsPad = rows % blockSize;
+ //columnsPad = columns % blockSize;
+ //qDebug() << "rowsPad:" << rowsPad << "columnsPad:" << columnsPad;
+ }
+ cv::Mat muGaussfilt, gray_Dot_gray, sigmaGaussfilt;
+ cv::GaussianBlur(gray_img, muGaussfilt, cv::Size(7, 7), 7 / 6, 7 / 6);
+ //sigma = sqrt(abs(imgaussfilt(ipImage.*ipImage, 7 / 6, 'FilterSize', 7, 'Padding', 'replicate') - mu.*mu));
+ //1.原图与自己点乘gray_Dot_gray=gray_img.*gray_img
+ //2.点乘结果进行高斯滤波GaussianBlur(gray_Dot_gray, gray_Dot_gray, cv::Size(7, 7), 7 / 6, 7 / 6)
+ //3.sqrt(abs(gray_Dot_gray- muGaussfilt.*muGaussfilt))
+ //1.原图与自己点乘
+ gray_Dot_gray = gray_img.mul(gray_img);
+ //2.点乘结果进行高斯滤波
+ cv::GaussianBlur(gray_Dot_gray, gray_Dot_gray, cv::Size(7, 7), 7 / 6, 7 / 6);
+ cv::Mat absDiff, imnorm;
+ //3.1
+ cv::absdiff(gray_Dot_gray, muGaussfilt.mul(muGaussfilt), absDiff);
+ //3.2
+ cv::sqrt(absDiff, sigmaGaussfilt);
+ //imnorm = (ipImage - mu). / (sigma + 1);
+ imnorm = (gray_img - muGaussfilt) / (sigmaGaussfilt + 1);
+
+ //NoticeableArtifactsMask = false(size(imnorm, 1), size(imnorm, 2));
+ //NoiseMask = false(size(imnorm, 1), size(imnorm, 2));
+ //ActivityMask = false(size(imnorm, 1), size(imnorm, 2));
+ //qDebug() << imnorm.rows << "," << imnorm.cols << ","<< imnorm.rows - blockSize * 2<<","
+ // << "rows%blockSize" << imnorm.rows % blockSize << "cols%blockSize" << imnorm.cols % blockSize;
+ for (int i = 0; i < imnorm.rows; i = i + blockSize)
+ {
+ for (int j = 0; j < imnorm.cols; j = j + blockSize)
+ {
+ int WNDC = 0, WNC = 0;
+ cv::Mat Block = imnorm(cv::Rect(j, i, blockSize, blockSize));
+ cv::Scalar meanScalar, stdDevScalar;
+ cv::meanStdDev(Block, meanScalar, stdDevScalar);
+ double blockStdDev = stdDevScalar.val[0];
+ double blockVar = stdDevScalar.val[0] * stdDevScalar.val[0];
+ if (blockVar > activityThreshold)
+ {
+ double WHSA = 1.0;
+ NHSA = NHSA + 1;
+ //qDebug() << i << j << "var > activityThreshold"<< blockVar;
+
+ bool blockImpaired = noticeDistCriterion(Block, nSegments, blockSize - 1, windowSize,
+ blockImpairedThreshold, blockSize);
+ if (blockImpaired)
+ WNDC = 1;
+ double blockBeta = noiseCriterion(Block, blockSize - 1, blockVar, blockStdDev);
+ if (blockStdDev > 2 * blockBeta)
+ WNC = 1;
+ distBlockScores = distBlockScores + WHSA * WNDC * (1 - blockVar) + WHSA * WNC * (blockVar);
+ }
+ }
+ }
+ //cv::resize(imnorm, imnorm, cv::Size(960, 540));
+ //cv::imshow("imnorm", imnorm);
+ double Score = ((distBlockScores + c) / (c + NHSA));//* 100
+
+ return Score;
+}
+
+bool NIIRS::noticeDistCriterion(cv::Mat Block, int nSegments, int blockSize_1,
+ int windowSize, float blockImpairedThreshold, int blockSize)
+{
+ cv::Mat topEdge, segTopEdge;
+ topEdge = Block.row(0);
+ segTopEdge = segmentEdge(topEdge, nSegments, blockSize_1, windowSize);
+
+ cv::Mat bottomEdge, segBottomEdge;
+ bottomEdge = Block.row(15);
+ segBottomEdge = segmentEdge(bottomEdge, nSegments, blockSize_1, windowSize);
+
+ cv::Mat leftEdge, segLeftEdge;
+ leftEdge = Block.col(0).t();//并转置
+ segLeftEdge = segmentEdge(leftEdge, nSegments, blockSize_1, windowSize);
+
+ cv::Mat rightEdge, segRightEdge;
+ rightEdge = Block.col(15).t();//并转置
+ segRightEdge = segmentEdge(rightEdge, nSegments, blockSize_1, windowSize);
+
+ bool blockImpaired = false;
+ //按行求标准差 [6col×11row]
+ //std::cout << segTopEdge.rows<< std::endl;
+ cv::Scalar meanScalar, stdDevScalar;
+ cv::Mat segRowTop, segRowBottom, segRowLeft, segRowRight;
+ vectorstdTop, stdBottom, stdLeft, stdRight;
+ for (int row = 0; row < segTopEdge.rows; row++)
+ {
+ //获取行
+ segRowTop = segTopEdge.row(row);
+ segRowBottom = segBottomEdge.row(row);
+ segRowLeft = segLeftEdge.row(row);
+ segRowRight = segRightEdge.row(row);
+ //求StdDev，并存在vector中
+ cv::meanStdDev(segRowTop, meanScalar, stdDevScalar);
+ stdTop.push_back(stdDevScalar.val[0]);
+ cv::meanStdDev(segRowBottom, meanScalar, stdDevScalar);
+ stdBottom.push_back(stdDevScalar.val[0]);
+ cv::meanStdDev(segRowLeft, meanScalar, stdDevScalar);
+ stdLeft.push_back(stdDevScalar.val[0]);
+ cv::meanStdDev(segRowRight, meanScalar, stdDevScalar);
+ stdRight.push_back(stdDevScalar.val[0]);
+ }
+ for (int vecsize = 0; vecsize < stdTop.size(); vecsize++)
+ {
+ if (stdTop[vecsize] < blockImpairedThreshold ||
+ stdBottom[vecsize] < blockImpairedThreshold ||
+ stdLeft[vecsize] < blockImpairedThreshold ||
+ stdRight[vecsize] < blockImpairedThreshold)
+ {
+ blockImpaired = true;
+ break;
+ }
+ }
+ return blockImpaired;
+}
+
+cv::Mat NIIRS::segmentEdge(cv::Mat blockEdge, int nSegments, int blockSize_1, int windowSize)
+{
+ cv::Mat segments(nSegments, windowSize, CV_64FC1);//空
+ //matlab:
+ //for i = 1:nSegments
+ // segments(i, :) = blockEdge(i : windowSize);
+ // if (windowSize <= (blockSize + 1))
+ // windowSize = windowSize + 1;
+ // end
+ //end
+ for (int row = 0; row < segments.rows; row++)
+ {
+ for (int col = 0; col < segments.cols; col++)
+ {
+ segments.at(row, col) = blockEdge.at(0, col + row);
+ }
+ }
+ //std::cout << "blockEdge" << std::endl << blockEdge << std::endl;
+ //std::cout << "segments" << std::endl << segments << std::endl;
+ return segments;
+}
+
+float NIIRS::noiseCriterion(cv::Mat Block, int blockSize_1, float blockVar, float blockStdDev)
+{
+ float blockBeta = 0.0, cenSurDev = 0.0;
+ cenSurDev = centerSurDev(Block, blockSize_1);
+ blockBeta = (abs(blockStdDev - cenSurDev)) / (max(blockStdDev, cenSurDev));
+
+ return blockBeta;
+}
+
+float NIIRS::centerSurDev(cv::Mat Block, int blockSize_1)
+{
+ int center1 = (blockSize_1 + 1) / 2;
+ int center2 = center1 + 1;
+
+ cv::Mat center(16, 2, CV_64FC1), surround(16, 14, CV_64FC1);
+
+ int surround_col = 0;
+ for (int col = 0; col < Block.cols; col++)
+ {
+ if (col == center1 - 1)//center col1
+ {
+ for (int row = 0; row < Block.rows; row++)
+ center.at(row, 0) = Block.at(row, col);
+ }
+ else if (col == center2 - 1)//center col2
+ {
+ for (int row = 0; row < Block.rows; row++)
+ center.at(row, 1) = Block.at(row, col);
+ }
+ else//surround
+ {
+ for (int row = 0; row < Block.rows; row++)
+ surround.at(row, surround_col) = Block.at(row, col);
+ surround_col = surround_col + 1;
+ }
+ }
+ //std::cout << "Block" << std::endl << Block << std::endl
+ // << "center" << std::endl << center << std::endl
+ // << "surround" << std::endl << surround << std::endl;
+ cv::Scalar meanScalar, stdDevScalar;
+ float stdDevCenter, stdDevSurround;
+ cv::meanStdDev(center, meanScalar, stdDevScalar);
+ stdDevCenter = stdDevScalar.val[0];
+ cv::meanStdDev(surround, meanScalar, stdDevScalar);
+ stdDevSurround = stdDevScalar.val[0];
+ float cenSurDev = stdDevCenter / stdDevSurround;
+
+ return cenSurDev;
+}
+
+double NIIRS::calculateMTFC_G(cv::Mat kernel)
+{
+ double sqr_sum = 0.0;
+ cv::Mat floatkernel = kernel;
+ floatkernel.convertTo(floatkernel, CV_64FC1);
+
+ for (int row = 0; row < floatkernel.rows; row++)
+ for (int col = 0; col < floatkernel.cols; col++)
+ sqr_sum = sqr_sum + floatkernel.at(row, col) * floatkernel.at(row, col);
+
+ double MTFC_G = sqrt(sqr_sum);
+
+ return MTFC_G;
+}
diff --git a/NiirsLinux/niirs.h b/NiirsLinux/niirs.h
new file mode 100644
index 0000000..8b28355
--- /dev/null
+++ b/NiirsLinux/niirs.h
@@ -0,0 +1,199 @@
+#pragma once
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+#include "ui_niirs.h"
+
+#include
+#include
+#include
+
+#include
+#include
+//#include
+//#include
+
+#include "gdal_priv.h"
+
+using namespace std;
+
+struct ImgInfo
+{
+ GDALDataType dataType; //标示数据类型
+ int bandCount; //数据波段数目
+ int xSize; //数据在x方向上的尺寸
+ int ySize; //数据在y方向上的尺寸
+ int xOffset; //x方向上的偏移；
+ int yOffset; //y方向上的偏移；
+ char projWkt[1024]; //数据的投影方式
+ double geoTransform[6]; //数据的仿射坐标
+ double noData; //数据的背景值
+};
+
+class NIIRS : public QMainWindow
+{
+ Q_OBJECT
+
+public:
+ NIIRS(QWidget* parent = Q_NULLPTR);
+
+ void setXmlInterface(QString strXmlInterface);
+
+ /* 计算NIIRS
+ * @return NIIRS */
+ double calculateNiirsScore();
+ //计算10个客观参数
+ void calculateObjectivePara();
+ //void calculateEvaluationImage();
+ double calculateRER();
+
+ /*
+ * @brief 获取图像基础信息并保存至结构体
+ * @param srcImgPath 输入文件路径
+ * @param imgInfo 结构体，保存获取的tiff文件信息
+ * @return 返回错误信息 */
+ int getImageInfo(const char* srcImgPath, struct ImgInfo* imgInfo);
+
+ /*
+ * @brief 按波段读取tiff
+ * @param srcImgPath 输入文件路径
+ * @param imgInfo 结构体，存储tiff文件信息
+ * @param imgArray 读取的tiff数组
+ * @param bandIndex 要读取的波段，bandIndex=0,则读取所有波段
+ * @return 返回错误信息 */
+ int readTiffFile(const char* srcImgPath, ImgInfo imgInfo, void* imgArray, int bandIndex);
+
+ void initTableWidget(int rowNum, int colNum);
+ void setTableWidgetValue(double val, int row, int col);
+
+ /*
+ * @brief 遍历文件夹下文件
+ * @param path 输入文件路径
+ * @param fileType 文件格式
+ * @return 文件路径absoluteFilePath,返回空值时为存在错误 */
+ QStringList getAllFiles(QString path, QString fileType);
+
+public slots:
+ void onRadioBtn_Visi();
+ void onRadioBtn_Infra();
+
+ void on_pbtOpenImg_clicked();
+ //void on_pbtOpenEdgeX_clicked();
+ //void on_pbtOpenEdgeY_clicked();
+
+private:
+ Ui::NIIRSClass ui;
+
+ /* 计算噪声增益
+ * @param kernel MTFC所使用的滤波器
+ * @return MTFC_G噪声增益 */
+ double calculateMTFC_G(cv::Mat kernel);
+ ////计算信噪比，mean/StdDev
+ //double calculateSNR(cv::Mat img);
+
+ /* 计算数字图像亮度
+ * @return 灰度图片平均亮度 */
+ double calculateMeanStdDev(cv::Mat img);
+ ////计算图像标准差方差
+ //float calculateStdDev(cv::Mat img);
+ //计算图像对比度
+ double calculateContrastRatio(cv::Mat img);
+ ////计算信息熵
+ //double calculateEntropy(cv::Mat img);
+ /*
+ * @brief 计算单波段的图像熵
+ * @param imgArray tiff数组
+ * @param maxVal tiff波段中的最大值
+ * @return 该波段的图像熵 */
+ double calculateTiffEntropy(unsigned short* bandBuf, int bufSize, int maxVal);
+ //计算相对边缘相应
+// double calculateERx(cv::Mat img);
+// double calculateERy(cv::Mat img);
+ //return pair.first:ER, pair.second:H
+ QPair calculateERx(cv::Mat img);
+ QPair calculateERy(cv::Mat img);
+
+ QPair calculate_rer_h(cv::Mat img, cv::Point pointA, cv::Point pointB);
+
+ cv::Point2d calculateFootPoint(array& ps);
+ //float SFRCalculation(cv::Mat& ROI, double gamma);
+ //void de_Gamma(cv::Mat& Src, double gamma);
+ //std::vector CentroidFind(cv::Mat& Src, std::vector& y_shifts, double* CCoffset);
+ //void SLR(std::vector& Cen_Shifts, std::vector& y_shifts, double* a, double* b);
+ //void ReduceRows(double slope, int* ImgHeight);
+ //std::vector OverSampling(cv::Mat& Src, double slope, double CCoffset, int height, int width, int* SamplingLen);
+ //std::vector HammingWindows(std::vector& deSampling, int SamplingLen);
+ //void DFT(std::vector& data, int size);
+
+
+ //陡度
+ //基于sobel算子，计算平均梯度
+ double gradient_Sobel(cv::Mat img);
+ //基于Laplacian算子，拉普拉斯导数，评价清晰度
+ double definition_Laplacian(cv::Mat img);
+ //四个方向的角二阶矩
+ void get_GLCM_0deg(cv::Mat& input, cv::Mat& dst);
+ void get_GLCM_90deg(cv::Mat& input, cv::Mat& dst);
+ void get_GLCM_45deg(cv::Mat& input, cv::Mat& dst);
+ void get_GLCM_135deg(cv::Mat& input, cv::Mat& dst);
+ //计算二阶矩的特征值
+ double calculateASM(cv::Mat& src);
+ //PIQE
+ double calculatePIQEScore(cv::Mat& img);
+ bool noticeDistCriterion(cv::Mat Block, int nSegments, int blockSize_1,
+ int windowSize, float blockImpairedThreshold, int blockSize);
+ cv::Mat segmentEdge(cv::Mat blockEdge, int nSegments, int blockSize_1, int windowSize);
+ float noiseCriterion(cv::Mat Block, int blockSize_1, float blockVar, float blockStdDev);
+ float centerSurDev(cv::Mat Block, int blockSize_1);
+
+ // 计算垂足的坐标，array：点、直线两端点
+ cv::Point calculate_foot_point(array& ps);
+ // 计算两点距离
+ float pointDistance(cv::Point pt1, cv::Point pt2);
+ //// 计算两直线交点
+ //cv::Point line_intersection(array& line1, array& line2);
+ //// 根据一条线的点pt、斜率k、另一条水平或垂直的直线的端点，求交点
+ //cv::Point point_k_edge_intersection(cv::Point pt, double k, array& line);
+ //// 计算点到直线距离
+ //int calculate_distance(array& ps);
+
+ bool jpegTest = false;
+ //int mBandCount = 0;
+
+ QString mXmlInterface = "";
+
+ //QString mXmlImgPath = "/home/testdata/19data/VALID-ALIGACE/JB19-2_MSS_000811113_MMB04_001_01_001_003_L1/JB19-2_MSS_000811113_MMB04_001_01_001_003_L1.tiff";
+ QString mXmlImgPath = "";//D:\\A_testdata\\niirs_test\\data\\GF2_mask1024.tif
+ QString mXmlTopLeftX = "0";//311,336
+ QString mXmlTopLeftY = "0";//324,293
+ QString mXmlBottomRightX = "0";//335,348
+ QString mXmlBottomRightY = "0";//337,311
+
+ QString mXmlProductDir = "";
+ QString mResultXmlDir = "";
+
+ cv::Point mRerEdgePointA;
+ cv::Point mRerEdgePointB;
+
+
+ double imgLight = 0.0; //亮度
+ double imgStdDev = 0.0; //标准差
+ double imgContrastRadio = 0.0; //对比度
+ double imgEntropy = 0.0; //熵
+ double imgSNR = 0.0; //信噪比
+ double imgGradient = 0.0; //陡度
+ double imgDefinition = 0.0; //清晰度
+ double imgASM = 0.0; //角二阶矩
+ double imgPIQE = 0.0; //PIQE
+
+ double imgRER = 0.986; //相对边缘响应
+
+};
diff --git a/NiirsLinux/niirs.qrc b/NiirsLinux/niirs.qrc
new file mode 100644
index 0000000..a4d4aa6
--- /dev/null
+++ b/NiirsLinux/niirs.qrc
@@ -0,0 +1,4 @@
+
+
+
+
diff --git a/NiirsLinux/niirs.ui b/NiirsLinux/niirs.ui
new file mode 100644
index 0000000..bbf83e7
--- /dev/null
+++ b/NiirsLinux/niirs.ui
@@ -0,0 +1,680 @@
+
+
+ NIIRSClass
+
+
+
+ 0
+ 0
+ 780
+ 520
+
+
+
+
+ 780
+ 520
+
+
+
+
+ 16777215
+ 16777215
+
+
+
+ NIIRS
+
+
+
+ -
+
+
-
+
+
+
+ 0
+ 0
+
+
+
+ 客观参数计算
+
+
+
-
+
+
+ Qt::Horizontal
+
+
+ QSizePolicy::Fixed
+
+
+
+ 40
+ 20
+
+
+
+
+ -
+
+
-
+
+
+
+ 0
+ 0
+
+
+
+
+ 16777215
+ 23
+
+
+
+ 左上角X坐标
+
+
+
+ -
+
+
+
+ 60
+ 16777215
+
+
+
+ 浏览
+
+
+
+ -
+
+
-
+
+
+ -
+
+
+
+ 0
+ 0
+
+
+
+
+ 16777215
+ 23
+
+
+
+ 左上角Y坐标
+
+
+
+ -
+
+
+
+
+ -
+
+
+
+ 0
+ 0
+
+
+
+
+ 16777215
+ 23
+
+
+
+ 待评价影像
+
+
+
+ -
+
+
-
+
+
+ -
+
+
+
+ 0
+ 0
+
+
+
+
+ 16777215
+ 23
+
+
+
+ 右下角Y坐标
+
+
+
+ -
+
+
+
+
+ -
+
+
+
+ 0
+ 0
+
+
+
+
+ 16777215
+ 23
+
+
+
+ 右下角X坐标
+
+
+
+ -
+
+
+
+
+ -
+
+
+ Qt::Horizontal
+
+
+
+ 40
+ 20
+
+
+
+
+ -
+
+
+
+ 0
+ 0
+
+
+
+
+ -
+
+
+ Qt::Horizontal
+
+
+
+ 40
+ 20
+
+
+
+
+ -
+
+
+
+ 90
+ 0
+
+
+
+
+ 90
+ 16777215
+
+
+
+ 计算客观参数
+
+
+
+
+
+
+ -
+
+
-
+
+
+
+ 0
+ 20
+
+
+
+
+ -
+
+
+ 异常信息:
+
+
+
+
+
+ -
+
+
+
+ 0
+ 0
+
+
+
+
+ 340
+ 16777215
+
+
+
+ NIIRS计算
+
+
+
-
+
+
-
+
+
+ NIIRS等级计算结果:
+
+
+
+ -
+
+
+
+ 0
+ 20
+
+
+
+
+ -
+
+
+ 计算
+
+
+
+
+
+ -
+
+
+ 调制传递函数补偿(MTFC)影像参数
+
+
+
-
+
+
-
+
+
+ 噪 声 增 益G
+
+
+ Qt::AlignLeading|Qt::AlignLeft|Qt::AlignVCenter
+
+
+
+ -
+
+
+
+ 0
+ 20
+
+
+
+
+
+
+
+
+
+ -
+
+
+ 信噪比(SNR)
+
+
+
-
+
+
-
+
+
+
+ 0
+ 20
+
+
+
+
+ -
+
+
+ SNR
+
+
+
+
+
+
+
+
+ -
+
+
+ 分辨率(GSD)
+
+
+
-
+
+
-
+
+
+ GSD_X
+
+
+
+ -
+
+
+
+ 0
+ 20
+
+
+
+
+ -
+
+
+ 米
+
+
+
+ -
+
+
+ GSD_Y
+
+
+
+ -
+
+
+
+ 0
+ 20
+
+
+
+
+ -
+
+
+ 米
+
+
+
+
+
+
+
+
+ -
+
+
+ 边缘响应及过冲
+
+
+
-
+
+
-
+
+
+
+ 0
+ 20
+
+
+
+
+ 60
+ 16777215
+
+
+
+
+ -
+
+
+
+ 0
+ 20
+
+
+
+
+ 60
+ 16777215
+
+
+
+
+ -
+
+
+
+ 40
+ 16777215
+
+
+
+ 过冲H
+
+
+
+ -
+
+
+
+ 40
+ 16777215
+
+
+
+ RER
+
+
+
+
+
+
+
+
+ -
+
+
+ 参数(PIQE)
+
+
+
-
+
+
-
+
+
+ PIQE
+
+
+
+ -
+
+
+
+ 0
+ 20
+
+
+
+
+
+
+
+
+
+ -
+
+
-
+
+
+ 可见光数据
+
+
+ true
+
+
+ buttonGroup
+
+
+
+ -
+
+
+ 中红外数据
+
+
+ buttonGroup
+
+
+
+ -
+
+
+ 长波红外数据
+
+
+ buttonGroup
+
+
+
+
+
+ -
+
+
-
+
+
+
+ 16777215
+ 23
+
+
+
+ 计算公式:
+
+
+
+ -
+
+
+
+ 0
+ 60
+
+
+
+
+ 16777215
+ 150
+
+
+
+ <!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0//EN" "http://www.w3.org/TR/REC-html40/strict.dtd">
+<html><head><meta name="qrichtext" content="1" /><style type="text/css">
+p, li { white-space: pre-wrap; }
+</style></head><body style=" font-family:'SimSun'; font-size:9pt; font-weight:400; font-style:normal;">
+<p style="-qt-paragraph-type:empty; margin-top:0px; margin-bottom:0px; margin-left:0px; margin-right:0px; -qt-block-indent:0; text-indent:0px;"><br /></p></body></html>
+
+
+
+ -
+
+
-
+
+
+
+ 60
+ 16777215
+
+
+
+ 波段选择:
+
+
+
+ -
+
+
+
+ 0
+ 0
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ lineEdit_img
+ pbtOpenImg
+ radioVisible
+ radioMidInfrared
+ textFormula
+ lineEdit_GSDX
+ lineEdit_GSDY
+ lineEdit_RER
+ lineEdit_H
+ lineEdit_MTFC_G
+ lineEdit_SNR
+ lineEdit_PIQE
+ lineEdit_NIIRS
+ pbtCalculate
+ lineEdit_INFO
+
+
+
+
+
+
+