目錄
參考資料：
- MinFilter - Wolfram 語言與系統參考資料中心
- ImageFilter - Wolfram 語言與系統參考資料中心
- Streaming Maximum-Minimum Filter Using No More than Three Comparisons per Element
- [SSE圖像算法優化系列七：基於SSE實現的極速的矩形核腐蝕和膨脹（最大值和最小值）算法。]
1、算法簡述
1.1、MinFilter(MaxFilter) 算法簡述
MinFilter（MaxFilter）算法是用於對一維或多維數據進行濾波的算法，濾波的結果為原數據中對應位置領域r內的最小（最大）值。在數據的邊界處，使用較小（較大）的鄰域.。
1.2、MinFilter(MaxFilter) 快速算法簡述
對於MinFilter(MaxFilter)的快速算法，思想來自於這篇論文Streaming Maximum-Minimum Filter Using No More than Three Comparisons per Element。在網上找到了這張圖，但這個圖也沒有什么文字說明，並不是很清楚。
下面按照我實現的時候的思路，來說一下我的理解。
首先，對於一個多維的數據，都可以逐個維度進行處理。比如說一個圖片，也就是二維數據，可以先對每一行進行處理，然后再對每一列進行處理，這樣得到的結果與行列同時處理是一樣的。
假設r=1
原始數據 --> 逐行處理 --> 逐列處理
5 2 1 3 4 2 1 1 1 3 2 1 1 1 3
6 9 8 4 7 6 6 4 4 4 2 1 1 0 0
7 3 8 2 0 3 3 2 0 0 0 0 0 0 0
9 0 1 5 6 0 0 0 1 5 0 0 0 0 0
原始數據 --> 逐行列處理
5 2 1 3 4 2 1 1 1 3
6 9 8 4 7 2 1 1 0 0
7 3 8 2 0 0 0 0 0 0
9 0 1 5 6 0 0 0 0 0
因此算法的關鍵在於提高一行數據處理的效率。
這個算法的過程大概是這樣的：
1、首先遍歷一行數據中最左邊的r*2+1個數據，獲取最小值和最小值的位置。然后對左邊邊界部分的處理，直接賦最小值。
2、從r+1位置開始向后遍歷，一直到右邊界部分。
3、遍歷的時候，判斷上一次獲取的最小值索引minIndex，是否在當前位置的領域r以內。
如果不在，則遍歷當前位置的領域r范圍，找出最小值的位置。也可以先與當前位置領域r內最右邊的比較，如果最右邊的小於minIndex位置的值，則minIndex就是這最右邊的這個，否則就需要遍歷當前位置領域r范圍內。
如果在，則說明當前位置領域r內，除了最右邊的元素，肯定都小於minIndex處的值。因為minIndex是當前位置上一個的領域r內的最小值，而上一個位置的領域r范圍與當前位置的領域r范圍只偏移了一個位置。
2、實現代碼
我這里實現了對一行數據的過濾，然后在一行數據過濾的基礎上實現對二維矩陣進行過濾。
對於MaxFilter的相關實現，只需要將下面對應的>=改為<=即可。
2.1、MinFilterOneRow 單行濾波代碼
/**********************************************************************//**
* @brief 對一行數據進行濾波，每個值用鄰域 r 內的最小值替換.
* @author solym@sohu.com/ymwh@foxmail.com
* @date 2019/4/23
* @param srcData 待濾波數據地址.
* @param srcChanelCount 待濾波數據每個像素的通道數.
* @param srcChanelIndex 待濾波數據要進行濾波的通道[0,srcChanelCount).
* @param dstData 濾波后輸出數據地址.
* @param dstChanelCount 濾波后輸出數據每個像素的通道數.
* @param dstChanelIndex 濾波后數據要輸出的通道索引[0,srcChanelCount).
* @param colnumCount 該行數據要濾波的像素數.
* @param radius 濾波的半徑大小.
*************************************************************************/
template<typename PixelDataType>
void MinFilterOneRow(
PixelDataType* srcData, const size_t srcChanelCount, const size_t srcChanelIndex,
PixelDataType* dstData, const size_t dstChanelCount, const size_t dstChanelIndex,
const size_t colnumCount, const size_t radius)
{
PixelDataType* pSrc = srcData;
PixelDataType* pDst = dstData;
size_t minIndex = 0; // 記錄最小值的下標
size_t blockSize = radius * 2 + 1; // 塊大小，以當前點為中心，左右各radius的寬度
PixelDataType minValue = pSrc[srcChanelIndex]; // 比較中獲取最小值進行記錄
// 對第一個塊進行處理（i從1開始，比較(i-1,i)位置像素值）
// 找出最小值(第一個塊內的最小值，就是r位置(塊中心)處的輸出值)
for(size_t iPixel=1; iPixel < blockSize; ++iPixel){
PixelDataType value = pSrc[iPixel*srcChanelCount + srcChanelIndex];
// 使用 >= 比 > 更快推進minIndex向前走
if(minValue >= value){
minValue = value;
minIndex = iPixel;
}
}
// 輸出到第一個塊中心(r)位置處的值。
// 它已經是第一個塊內的最小值，也就是該塊左邊都只能是這個值
for(size_t i=0;i<=radius;++i){
pDst[i*dstChanelCount + dstChanelIndex] = minValue;
}
// 開始處理r+1位置之后的值
for (size_t iPixel = radius + 1; iPixel < colnumCount - radius; ++iPixel) {
/* i-r i i+r
* |____________|___________|_
* └min
* 當前最小的索引在當前位置為中心的塊的內(一定位於當前塊內或前一個)
* iPixel是當前塊的中心，下面說的當前位置都指iPixel
*/
if(minIndex >= (iPixel - radius)) {
// 當前最小索引位置值與當前位置為中心的塊的最后一個值比較
// 根據下面的代碼可知，如果mIndex在塊的內部，它所在位置的值一定是最小的
// 進入本次循環時，minIndex是上次比較的值，而上一個塊與當前塊等長，位置差一位
// 所以可以直接和當前塊最后一個像素值進行比較了，當前塊也就完全比較完了
size_t nextBlockFirstIndex = iPixel + radius;
if(pSrc[minIndex*srcChanelCount + srcChanelIndex] >
pSrc[nextBlockFirstIndex*srcChanelCount + srcChanelIndex]){
// 賦值當前最小值索引和值
minIndex = nextBlockFirstIndex;
minValue = pSrc[nextBlockFirstIndex*srcChanelCount + srcChanelIndex];
}
}else{
// 如果不在當前位置為中心的塊內，則對當前塊進行查找最小值
// 則將minIndex設置該塊的最左邊位置
minIndex = iPixel - radius;
// 獲取當前位置為中心的塊的最小值和索引
minValue = pSrc[minIndex*srcChanelCount + srcChanelIndex];
size_t blockEnd = minIndex + blockSize;
for (size_t iBPixel = minIndex; iBPixel < blockEnd; ++iBPixel) {
PixelDataType value = pSrc[iBPixel*srcChanelCount + srcChanelIndex];
if(minValue >= value){
minIndex = iBPixel;
minValue = value;
}
}
} // end if minIndex > ...
pDst[iPixel*dstChanelCount + dstChanelIndex] = minValue;
} // end for iPixel
// 最后一個塊中心位置的右邊，一定都是和它中心位置的值是一樣的
for (size_t i = colnumCount-radius; i < colnumCount; ++i) {
pDst[i*dstChanelCount + dstChanelIndex] = minValue;
}
}
2.2、MinFilterOneMatrix 單個二維矩陣濾波代碼
對於二維矩陣進行濾波，實際上是先進行行濾波，然后結果進行行列轉置，對轉置的結果再次進行行濾波，然后再行列轉置輸出。
/**********************************************************************//**
* @brief 對一個矩陣數據進行濾波，每個值用鄰域 r 內的最小值替換.
* @author solym@sohu.com/ymwh@foxmail.com
* @date 2019/4/23
* @param srcData 待濾波數據地址.
* @param srcBytePerRow 待濾波數據每行的字節數.
* @param srcChanelCount 待濾波數據每個像素的通道數.
* @param srcChanelIndex 待濾波數據要進行濾波的通道[0,srcChanelCount).
* @param dstData 濾波后輸出數據地址.
* @param dstBytePerRow 濾波后輸出數據每行的字節數.
* @param dstChanelCount 濾波后輸出數據每個像素的通道數.
* @param dstChanelIndex 濾波后數據要輸出的通道索引[0,srcChanelCount).
* @param rowCount 矩陣的行數.
* @param colCount 矩陣的列數.
* @param radius 濾波的半徑大小.
*************************************************************************/
template<typename PixelDataType>
void MinFilterOneMatrix(
PixelDataType* srcData, const size_t srcBytePerRow,
const size_t srcChanelCount, const size_t srcChanelIndex,
PixelDataType* dstData, const size_t dstBytePerRow,
const size_t dstChanelCount, const size_t dstChanelIndex,
const size_t rowCount, const size_t colCount,
const size_t radius)
{
unsigned char* pSrc = reinterpret_cast<unsigned char*>(srcData);
unsigned char* pDst = reinterpret_cast<unsigned char*>(dstData);
// 保存中間結果
std::vector<PixelDataType> tmpData(rowCount * colCount);
// 逐行進行濾波
for (size_t row = 0; row < rowCount; ++row) {
// 獲取輸入和輸出每行的行首位置
PixelDataType* pSrcRowFirst = (PixelDataType*)(pSrc + row * srcBytePerRow);
PixelDataType* pDstRowFirst = tmpData.data() + row * colCount;
// 對當前行進行濾波
MinFilterOneRow<PixelDataType>(pSrcRowFirst,srcChanelCount,srcChanelIndex,
pDstRowFirst,1,0,
colCount,radius);
}
// 將行濾波后的結果進行 行列轉置（進行列濾波）
std::vector<PixelDataType> tmpDataT(rowCount * colCount);
for (size_t row = 0; row < rowCount; ++row) {
for(size_t col = 0; col < colCount; ++col){
tmpDataT[col*rowCount + row] = tmpData[row*colCount + col];
}
}
// 對轉置后的矩陣進行 逐行濾波（就是原行濾波后結果進行列濾波）
for (size_t col = 0; col < colCount; ++col) {
PixelDataType* pSrcColFirst = tmpDataT.data() + col * rowCount;
PixelDataType* pDstColFirst = tmpData.data() + col * rowCount;
// 對當前行進行濾波
MinFilterOneRow<PixelDataType>(pSrcColFirst,1,0,
pDstColFirst,1,0,
rowCount,radius);
}
// 將行列濾波后的結果輸出
for (size_t row = 0; row < rowCount; ++row) {
PixelDataType* pDstRowFirst = (PixelDataType*)(pDst + row * dstBytePerRow);
for(size_t col = 0; col < colCount; ++col){
pDstRowFirst[col*dstChanelCount+dstChanelIndex] = tmpData[col*rowCount + row];
}
}
}
3、測試
3.1 測試截圖
使用Qt寫了一個簡單的測試程序進行測試，測試結果如下：
3.2 測試代碼
#include "filter.hpp"
#include <QApplication>
#include <QWidget>
#include <QLineEdit>
#include <QPushButton>
#include <QComboBox>
#include <QVBoxLayout>
#include <QHBoxLayout>
#include <QFileDialog>
#include <QWebEngineView>
#include <QXmlStreamWriter>
#include <QBuffer>
int main(int argc, char *argv[])
{
QApplication a(argc, argv);
QImage srcImage,dstImage;
// 創建窗口
QWidget widget;
// 添加控件
QWebEngineView *wevView = new QWebEngineView(&widget);
QLineEdit* leSrcImagePath = new QLineEdit(&widget);
QPushButton* pbSelectSrcFile = new QPushButton(QStringLiteral("選擇圖片"),&widget);
QComboBox* cbSelectFilterAlg = new QComboBox(&widget);
cbSelectFilterAlg->addItems(
{ QStringLiteral("MinFilter"),QStringLiteral("MaxFilter")});
QPushButton* pbRunFilter = new QPushButton(QStringLiteral("執行濾波"),&widget);
//pbRunDetect->setEnabled(false);
QHBoxLayout* hbLayout = new QHBoxLayout;
// 設置布局
hbLayout->addWidget(leSrcImagePath);
hbLayout->addWidget(pbSelectSrcFile);
hbLayout->addWidget(cbSelectFilterAlg);
hbLayout->addWidget(pbRunFilter);
QVBoxLayout* vbLayout = new QVBoxLayout(&widget);
vbLayout->addLayout(hbLayout);
vbLayout->addWidget(wevView);
// 添加處理操作
std::function<void(QString,const QImage&,const QImage&)>
updateHtmlView =
[wevView,leSrcImagePath](QString filterName,const QImage& srcImage,const QImage& resultimage)
{
QString tmpPath;
QByteArray html;
{
QXmlStreamWriter writer(&html);
writer.setAutoFormatting(true);
writer.writeStartDocument();
writer.writeStartElement("html");
writer.writeStartElement("body");
writer.writeAttribute("bgcolor","gray");
if(!srcImage.isNull()){
writer.writeTextElement("h2",QStringLiteral("原圖"));
writer.writeStartElement("img");
QBuffer buffer;
srcImage.save(&buffer,"PNG");
writer.writeAttribute("src","data:image/png;base64," + buffer.data().toBase64());
// writer.writeAttribute("src",QUrl(leSrcImagePath->text()).toString());
writer.writeEndElement();
}
if(!resultimage.isNull()){
writer.writeTextElement("h2",filterName + QStringLiteral("濾波結果圖"));
writer.writeStartElement("img");
QBuffer buffer;
resultimage.save(&buffer,"PNG");
writer.writeAttribute("src","data:image/png;base64," + buffer.data().toBase64());
// tmpPath = QDir::tempPath() + QString::fromUtf8(".png");
// resultimage.save(tmpPath,"PNG");
// writer.writeAttribute("src",QUrl(tmpPath).toString());
writer.writeEndElement();
}
writer.writeEndElement();
writer.writeEndElement();
}
wevView->setHtml(QString::fromUtf8(html));
// if(!resultimage.isNull()){
// qDebug()<<tmpPath;
// QFile::remove(tmpPath);
// }
};
// 文件選擇按鈕單擊信號處理
QObject::connect(pbSelectSrcFile,&QPushButton::clicked,
[leSrcImagePath,&srcImage,&widget,&updateHtmlView]()
{
static QString path(".");
path = QFileDialog::getOpenFileName(&widget,
QStringLiteral("選擇待濾波圖片"),
path,
QString("Images (*.png *.jpg *.jpeg *.jfif)"));
if(path.isEmpty()){return;}
QImage image;
if(!image.load(path)){return;}
srcImage = image.convertToFormat(QImage::Format_RGBA8888);
leSrcImagePath->setText(path);
updateHtmlView(QString(),srcImage,QImage());
});
// 執行濾波按鈕單擊信號處理
QObject::connect(pbRunFilter,&QPushButton::clicked,
[&cbSelectFilterAlg,&srcImage,&dstImage,&updateHtmlView]
{
// 獲取濾波算法名稱
QString filterName = cbSelectFilterAlg->currentText();
// 最小值濾波 https://reference.wolfram.com/language/ref/MinFilter.html
if(filterName == QStringLiteral("MinFilter")){
dstImage = srcImage;
unsigned int raduis = 2;
uchar* pSrc = srcImage.bits();
uchar* pDst = dstImage.bits();
unsigned int colCount = srcImage.width();
unsigned int rowCount = srcImage.height();
unsigned int chanel = 4;
// 分別對RGB通道進行濾波
MinFilterOneMatrix<uchar>(pSrc,srcImage.bytesPerLine(),chanel,0,
pDst,dstImage.bytesPerLine(),chanel,0,
rowCount,colCount,raduis);
MinFilterOneMatrix<uchar>(pSrc,srcImage.bytesPerLine(),chanel,1,
pDst,dstImage.bytesPerLine(),chanel,1,
rowCount,colCount,raduis);
MinFilterOneMatrix<uchar>(pSrc,srcImage.bytesPerLine(),chanel,2,
pDst,dstImage.bytesPerLine(),chanel,2,
rowCount,colCount,raduis);
}
else if(filterName == QStringLiteral("MaxFilter")){
dstImage = srcImage;
unsigned int raduis = 2;
uchar* pSrc = srcImage.bits();
uchar* pDst = dstImage.bits();
unsigned int colCount = srcImage.width();
unsigned int rowCount = srcImage.height();
unsigned int chanel = 4;
// 分別對RGB通道進行濾波
MaxFilterOneMatrix<uchar>(pSrc,srcImage.bytesPerLine(),chanel,0,
pDst,dstImage.bytesPerLine(),chanel,0,
rowCount,colCount,raduis);
MaxFilterOneMatrix<uchar>(pSrc,srcImage.bytesPerLine(),chanel,1,
pDst,dstImage.bytesPerLine(),chanel,1,
rowCount,colCount,raduis);
MaxFilterOneMatrix<uchar>(pSrc,srcImage.bytesPerLine(),chanel,2,
pDst,dstImage.bytesPerLine(),chanel,2,
rowCount,colCount,raduis);
}
updateHtmlView(filterName,srcImage,dstImage);
});
widget.resize(1024,768);
widget.show();
return a.exec();
}