CLAHE算法的OpenCV重构及详细解读
opencv
OpenCV: 开源计算机视觉库
项目地址：https://gitcode.com/gh_mirrors/opencv31/opencv
免费下载资源
·
OpenCV封装的CLAHE算法阅读起来比较麻烦，作者使用OpenCV在完全遵循原代码逻辑的基础上对该算法进行了重构，完整实现了OpenCV4.6.0中的CLAHE算法（包括16位图像和8位图像），代码如下：
void CLAHE16UC1(const cv::Mat& src,cv::Mat& dst,double clip = 40.0, cv::Size tiles = cv::Size(8,8));
void CLAHE8UC1(const cv::Mat& src,cv::Mat& dst,double clip = 40.0, cv::Size tiles = cv::Size(8,8));
void CLAHE16UC1(const cv::Mat& src,cv::Mat& dst,double clip, cv::Size tiles)
{
if(src.type() != CV_16UC1)
return;
int histSize = 65536;
//图像被划分为tilesX×tilesY个分块
int tilesX = tiles.width;
int tilesY = tiles.height;
cv::Size tileSize;
cv::Mat srcForLut;
if (src.size().width % tilesX == 0 && src.size().height % tilesY == 0)
{
tileSize = cv::Size(src.size().width / tilesX, src.size().height / tilesY);
srcForLut = src.clone();
}
else
{
cv::Mat srcExt;
//图像边缘填充，保证用于计算LUT的图像其宽度和高度可以分别被tilesX和tilesY整除
cv::copyMakeBorder(src, srcExt,0, tilesY - (src.size().height % tilesY),0, tilesX - (src.size().width % tilesX), cv::BORDER_REFLECT_101);
tileSize = cv::Size(srcExt.size().width / tilesX, srcExt.size().height / tilesY);
srcForLut = srcExt.clone();
}
//计算直方图“削峰填谷”的阈值，计算clipLimit的过程可以理解为：局部分块图像直方图平均频数的clip倍
int clipLimit = static_cast<int>(clip * tileSize.area() / histSize);
clipLimit = std::max(clipLimit, 1);
cv::Mat lut(tilesX*tilesY, histSize,src.type()); //存储每个分块的LUT，则共有tilesX*tilesY个LUT，每个LUT有histSize个元素
lut.setTo(0);
int* tileHist = new int[histSize];
//对每个块分别计算LUT，存在上面的lut里
for(int ty=0;ty<tilesY;ty++)
{
for(int tx=0;tx<tilesX;tx++)
{
memset(tileHist,0,histSize*sizeof(int));
ushort* tileLut = lut.ptr<ushort>(ty*tilesX+tx);
cv::Rect tileRect;
tileRect.x = tx * tileSize.width;
tileRect.y = ty * tileSize.height;
tileRect.width = tileSize.width;
tileRect.height = tileSize.height;
const cv::Mat tile = srcForLut(tileRect).clone();
const ushort* pTile = tile.ptr<ushort>();
//统计局部分块的直方图
for(int i=0;i<tile.rows;i++)
{
for(int j=0;j<tile.cols;j++)
{
tileHist[pTile[i*tile.cols+j]]++;
}
}
if (clipLimit > 0)
{
//直方图“削峰”，并统计削去的像素数
int clipped = 0;
for (int i = 0; i < histSize; ++i)
{
if (tileHist[i] > clipLimit)
{
clipped += (tileHist[i] - clipLimit);
tileHist[i] = clipLimit;
}
}
//将削去的像素个数平均分配给直方图的每个灰阶
int redistBatch = clipped / histSize;
int residual = clipped - redistBatch * histSize;
//将削去的像素个数平均分配给直方图的每个灰阶，1）首先平均分配
for (int i = 0; i < histSize; ++i)
tileHist[i] += redistBatch;
//将削去的像素个数平均分配给直方图的每个灰阶，2）仍然有剩余的，直方图上等间距分配
if (residual != 0)
{
int residualStep = MAX(histSize / residual, 1);
for (int i = 0; i < histSize && residual > 0; i += residualStep, residual--)
tileHist[i]++;
}
}
const float lutScale = cv::saturate_cast<float>(histSize - 1) / (tile.cols*tile.rows);
//计算LUT
int sum = 0;
for (int i = 0; i < histSize; ++i)
{
sum += tileHist[i];
tileLut[i] = cv::saturate_cast<ushort>(sum * lutScale);
}
}
}
delete [] tileHist;
tileHist = NULL;
float invTileW = 1.0f/tileSize.width;
float invTileH = 1.0f/tileSize.height;
dst.create(src.size(),CV_16UC1);
ushort* pDst = dst.ptr<ushort>();
const ushort* pSrc = src.ptr<ushort>();
//对每个像素点，找到其4个近邻分块的LUT所映射的像素值，并进行双线性插值
for(int i=0;i<src.rows;i++)
{
for(int j=0;j<src.cols;j++)
{
ushort pix = pSrc[i*src.cols+j];
float txf = j*invTileW - 0.5f;
int tx1 = cvFloor(txf);
int tx2 = tx1 + 1; //注意这里，说明4个近邻分块是紧挨着的
float px1 = txf - tx1;
float px2 = 1.0f - px1;
tx1 = std::max(tx1, 0);
tx2 = std::min(tx2, tilesX - 1);
float tyf = i*invTileH - 0.5f;
int ty1 = cvFloor(tyf);
int ty2 = ty1 + 1; //注意这里，说明4个近邻分块是紧挨着的
float py1 = tyf - ty1;
float py2 = 1.0f - py1;
ty1 = std::max(ty1, 0);
ty2 = std::min(ty2, tilesY - 1);
ushort* tileLuty1x1 = lut.ptr<ushort>(ty1*tilesX+tx1);
ushort* tileLuty1x2 = lut.ptr<ushort>(ty1*tilesX+tx2);
ushort* tileLuty2x1 = lut.ptr<ushort>(ty2*tilesX+tx1);
ushort* tileLuty2x2 = lut.ptr<ushort>(ty2*tilesX+tx2);
//4个邻块的映射灰度值线性插值，先x方向，再y方向
//注意：前面提到px1+px2=1.0和py1+py2=1.0，这里的px1，px2，py1，py2都是距离。在作为双线性插值的权重时，距离近的权重应大；反之亦然。
//以x方向插值举例：距离为px1时，其权重应取1.0-px1，即px2；距离为px2时，其权重应取1.0-px2，即px1。
pDst[i*src.cols+j] = cv::saturate_cast<ushort>(
(tileLuty1x1[pix]*px2+tileLuty1x2[pix]*px1)*py2+
(tileLuty2x1[pix]*px2+tileLuty2x2[pix]*px1)*py1);
}
}
}
void CLAHE8UC1(const cv::Mat& src,cv::Mat& dst,double clip, cv::Size tiles)
{
if(src.type() != CV_8UC1)
return;
int histSize = 256;
//图像被划分为tilesX×tilesY个分块
int tilesX = tiles.width;
int tilesY = tiles.height;
cv::Size tileSize;
cv::Mat srcForLut;
if (src.size().width % tilesX == 0 && src.size().height % tilesY == 0)
{
tileSize = cv::Size(src.size().width / tilesX, src.size().height / tilesY);
srcForLut = src.clone();
}
else
{
cv::Mat srcExt;
//图像边缘填充，保证用于计算LUT的图像其宽度和高度可以分别被tilesX和tilesY整除
cv::copyMakeBorder(src, srcExt,0, tilesY - (src.size().height % tilesY),0, tilesX - (src.size().width % tilesX), cv::BORDER_REFLECT_101);
tileSize = cv::Size(srcExt.size().width / tilesX, srcExt.size().height / tilesY);
srcForLut = srcExt.clone();
}
//计算直方图“削峰填谷”的阈值，计算clipLimit的过程可以理解为：局部分块图像直方图平均频数的clip倍
int clipLimit = static_cast<int>(clip * tileSize.area() / histSize);
clipLimit = std::max(clipLimit, 1);
cv::Mat lut(tilesX*tilesY, histSize,src.type()); //存储每个分块的LUT，则共有tilesX*tilesY个LUT，每个LUT有histSize个元素
lut.setTo(0);
int* tileHist = new int[histSize];
//对每个块分别计算LUT，存在上面的lut里
for(int ty=0;ty<tilesY;ty++)
{
for(int tx=0;tx<tilesX;tx++)
{
memset(tileHist,0,histSize*sizeof(int));
uchar* tileLut = lut.ptr<uchar>(ty*tilesX+tx);
cv::Rect tileRect;
tileRect.x = tx * tileSize.width;
tileRect.y = ty * tileSize.height;
tileRect.width = tileSize.width;
tileRect.height = tileSize.height;
const cv::Mat tile = srcForLut(tileRect).clone();
const uchar* pTile = tile.ptr<uchar>();
//统计局部分块的直方图
for(int i=0;i<tile.rows;i++)
{
for(int j=0;j<tile.cols;j++)
{
tileHist[pTile[i*tile.cols+j]]++;
}
}
if (clipLimit > 0)
{
//直方图“削峰”，并统计削去的像素数
int clipped = 0;
for (int i = 0; i < histSize; ++i)
{
if (tileHist[i] > clipLimit)
{
clipped += (tileHist[i] - clipLimit);
tileHist[i] = clipLimit;
}
}
//将削去的像素个数平均分配给直方图的每个灰阶
int redistBatch = clipped / histSize;
int residual = clipped - redistBatch * histSize;
//将削去的像素个数平均分配给直方图的每个灰阶，1）首先平均分配
for (int i = 0; i < histSize; ++i)
tileHist[i] += redistBatch;
//将削去的像素个数平均分配给直方图的每个灰阶，2）仍然有剩余的，直方图上等间距分配
if (residual != 0)
{
int residualStep = MAX(histSize / residual, 1);
for (int i = 0; i < histSize && residual > 0; i += residualStep, residual--)
tileHist[i]++;
}
}
const float lutScale = cv::saturate_cast<float>(histSize - 1) / (tile.cols*tile.rows);
//计算LUT
int sum = 0;
for (int i = 0; i < histSize; ++i)
{
sum += tileHist[i];
tileLut[i] = cv::saturate_cast<uchar>(sum * lutScale);
}
}
}
delete [] tileHist;
tileHist = NULL;
float invTileW = 1.0f/tileSize.width;
float invTileH = 1.0f/tileSize.height;
dst.create(src.size(),CV_8UC1);
uchar* pDst = dst.ptr<uchar>();
const uchar* pSrc = src.ptr<uchar>();
//对每个像素点，找到其4个近邻分块的LUT所映射的像素值，并进行双线性插值
for(int i=0;i<src.rows;i++)
{
for(int j=0;j<src.cols;j++)
{
uchar pix = pSrc[i*src.cols+j];
float txf = j*invTileW - 0.5f;
int tx1 = cvFloor(txf);
int tx2 = tx1 + 1; //注意这里，说明4个近邻分块是紧挨着的
float px1 = txf - tx1;
float px2 = 1.0f - px1;
tx1 = std::max(tx1, 0);
tx2 = std::min(tx2, tilesX - 1);
float tyf = i*invTileH - 0.5f;
int ty1 = cvFloor(tyf);
int ty2 = ty1 + 1; //注意这里，说明4个近邻分块是紧挨着的
float py1 = tyf - ty1;
float py2 = 1.0f - py1;
ty1 = std::max(ty1, 0);
ty2 = std::min(ty2, tilesY - 1);
uchar* tileLuty1x1 = lut.ptr<uchar>(ty1*tilesX+tx1);
uchar* tileLuty1x2 = lut.ptr<uchar>(ty1*tilesX+tx2);
uchar* tileLuty2x1 = lut.ptr<uchar>(ty2*tilesX+tx1);
uchar* tileLuty2x2 = lut.ptr<uchar>(ty2*tilesX+tx2);
//4个邻块的映射灰度值线性插值，先x方向，再y方向
//注意：前面提到px1+px2=1.0和py1+py2=1.0，这里的px1，px2，py1，py2都是距离。在作为双线性插值的权重时，距离近的权重应大；反之亦然。
//以x方向插值举例：距离为px1时，其权重应取1.0-px1，即px2；距离为px2时，其权重应取1.0-px2，即px1。
pDst[i*src.cols+j] = cv::saturate_cast<uchar>(
(tileLuty1x1[pix]*px2+tileLuty1x2[pix]*px1)*py2+
(tileLuty2x1[pix]*px2+tileLuty2x2[pix]*px1)*py1);
}
}
}
代码中比较难理解的是双线性插值时4个近邻块的选择方式，这里根据作者的理解给出示意图（以图像被均分为6*6块为例，其中红色小块表示当前处理像素点所在位置，绿色大块表示所选择的用于双线性插值的4个近邻块的位置），如下所示：
使用上述代码里的默认参数，8位图像处理效果如下：
GitHub 加速计划 / opencv31 / opencv
77.39 K
55.71 K
下载
OpenCV: 开源计算机视觉库
最近提交(Master分支：2 个月前 )
48668119
dnn: use dispatching for Winograd optimizations 6 天前
3dace76c
flann: remove unused hdf5 header 6 天前
旨在为数千万中国开发者提供一个无缝且高效的云端环境，以支持学习、使用和贡献开源项目。
更多推荐
3
0- 0
已为社区贡献1条内容
所有评论(0)