fix(processors): 修复对比度与阈值算子16位分支的降位精度损失

- ContrastProcessor: CLAHE分支改用原生CvInvoke.CLAHE(支持CV_16U),不再16→8→16往返
- ThresholdProcessor: Otsu分支在16位完整直方图(65536 bin)上原生计算最优阈值
This commit is contained in:
wei.lw.li
2026-07-21 10:05:13 +08:00
parent 87cd56b87e
commit 85045d45e2
2 changed files with 70 additions and 122 deletions
@@ -71,25 +71,29 @@ public class ThresholdProcessor<TDepth> : ImageProcessorBase<TDepth>
if (useOtsu)
{
// Otsu 仅支持 8 位,16 位需先降位
if (typeof(TDepth) == typeof(ushort))
{
using var img8 = PixelDepthHelper.ToByteImage(inputImage);
using var res8 = new Image<Gray, byte>(inputImage.Size);
CvInvoke.Threshold(img8, res8, minThreshold >> 8, 255, ThresholdType.Otsu);
// 将 8 位二值结果升位到 16 位
var res16 = PixelDepthHelper.FromByteImage<TDepth>(res8);
_logger.Debug("Process: UseOtsu=true (16→8→16)");
return res16;
}
else
// 8 位直接用 OpenCV Otsu;16 位在完整位深直方图上原生计算 Otsu 阈值,避免降位
if (typeof(TDepth) == typeof(byte))
{
using var res8 = new Image<Gray, byte>(inputImage.Size);
CvInvoke.Threshold(inputImage as Image<Gray, byte> ?? inputImage.Convert<Gray, byte>(),
res8, minThreshold, 255, ThresholdType.Otsu);
_logger.Debug("Process: UseOtsu=true");
_logger.Debug("Process: UseOtsu=true (8bit native)");
return res8 as Image<Gray, TDepth> ?? PixelDepthHelper.FromByteImage<TDepth>(res8);
}
else
{
int otsuThreshold = ComputeOtsuThreshold16((inputImage as Image<Gray, ushort>)!, maxVal);
Parallel.For(0, height, y =>
{
for (int x = 0; x < width; x++)
{
int val = PixelDepthHelper.ReadPixel(inputImage, y, x);
PixelDepthHelper.WritePixel(result, y, x, val > otsuThreshold ? maxVal : 0);
}
});
_logger.Debug("Process: UseOtsu=true (16bit native), threshold={Threshold}", otsuThreshold);
return result;
}
}
// 手工双阈值分割(支持全位深)
@@ -106,4 +110,51 @@ public class ThresholdProcessor<TDepth> : ImageProcessorBase<TDepth>
_logger.Debug("Process: MinThreshold={Min}, MaxThreshold={Max}", minThreshold, maxThreshold);
return result;
}
/// <summary>
/// 在完整 16 位直方图(65536 bin)上计算 Otsu 最优阈值(最大化类间方差)。
/// 返回值为灰度阈值:像素值 &gt; 阈值 判为前景。
/// </summary>
private static int ComputeOtsuThreshold16(Image<Gray, ushort> image, int maxVal)
{
int levels = maxVal + 1;
var histogram = new long[levels];
int h = image.Height, w = image.Width;
var data = image.Data;
for (int y = 0; y < h; y++)
for (int x = 0; x < w; x++)
histogram[data[y, x, 0]]++;
long totalPixels = (long)w * h;
double totalSum = 0;
for (int i = 0; i < levels; i++)
totalSum += (double)i * histogram[i];
double bgSum = 0;
long bgPixels = 0;
double maxVariance = -1;
int bestThreshold = 0;
for (int t = 0; t < levels; t++)
{
bgPixels += histogram[t];
if (bgPixels == 0) continue;
long fgPixels = totalPixels - bgPixels;
if (fgPixels == 0) break;
bgSum += (double)t * histogram[t];
double bgMean = bgSum / bgPixels;
double fgMean = (totalSum - bgSum) / fgPixels;
double variance = (double)bgPixels * fgPixels * (bgMean - fgMean) * (bgMean - fgMean);
if (variance > maxVariance)
{
maxVariance = variance;
bestThreshold = t;
}
}
return bestThreshold;
}
}
@@ -95,19 +95,12 @@ public class ContrastProcessor<TDepth> : ImageProcessorBase<TDepth>
if (useCLAHE)
{
// CLAHE 仅对 8 位有原生 CvInvoke 支持,16 位降位处理
if (typeof(TDepth) == typeof(ushort))
{
using var img8 = PixelDepthHelper.ToByteImage(inputImage);
using var res8 = ApplyCLAHE(img8, clipLimit);
_logger.Debug("Process (CLAHE 16→8→16)");
return PixelDepthHelper.FromByteImage<TDepth>(res8);
}
else
{
var res = ApplyCLAHE(inputImage as Image<Gray, byte> ?? inputImage.Convert<Gray, byte>(), clipLimit);
return res as Image<Gray, TDepth> ?? PixelDepthHelper.FromByteImage<TDepth>(res);
}
// OpenCV CLAHE 原生支持 CV_8UC1 和 CV_16UC1,直接按位深处理,无需降位
var result = new Image<Gray, TDepth>(inputImage.Size);
CvInvoke.CLAHE(inputImage, clipLimit, new Size(8, 8), result);
_logger.Debug("Process (CLAHE native {Depth}bit): ClipLimit={ClipLimit}",
typeof(TDepth) == typeof(ushort) ? 16 : 8, clipLimit);
return result;
}
else if (autoContrast)
{
@@ -140,100 +133,4 @@ public class ContrastProcessor<TDepth> : ImageProcessorBase<TDepth>
_logger.Debug("AutoContrastStretch: min={Min}, max={Max}", minVal, maxVal);
return PixelDepthHelper.FromFloatImage<TDepth>(floatImage);
}
private Image<Gray, byte> ApplyCLAHE(Image<Gray, byte> inputImage, double clipLimit)
{
int tileSize = 8;
int width = inputImage.Width;
int height = inputImage.Height;
byte[,,] srcData = inputImage.Data;
// 计算分块数
int tilesX = (width + tileSize - 1) / tileSize;
int tilesY = (height + tileSize - 1) / tileSize;
int actualTileW = (width + tilesX - 1) / tilesX;
int actualTileH = (height + tilesY - 1) / tilesY;
// 为每个 tile 计算带 clip limit 的均衡化映射表
var luts = new byte[tilesY, tilesX, 256];
for (int ty = 0; ty < tilesY; ty++)
{
for (int tx = 0; tx < tilesX; tx++)
{
int x0 = tx * actualTileW;
int y0 = ty * actualTileH;
int x1 = Math.Min(x0 + actualTileW, width);
int y1 = Math.Min(y0 + actualTileH, height);
int tilePixels = (x1 - x0) * (y1 - y0);
// 构建直方图
var hist = new int[256];
for (int y = y0; y < y1; y++)
for (int x = x0; x < x1; x++)
hist[srcData[y, x, 0]]++;
// Clip limit 裁剪并重新分配
int clipThreshold = (int)(clipLimit * tilePixels / 256);
if (clipThreshold > 0)
{
int excess = 0;
for (int i = 0; i < 256; i++)
{
if (hist[i] > clipThreshold)
{
excess += hist[i] - clipThreshold;
hist[i] = clipThreshold;
}
}
int avgInc = excess / 256;
int remainder = excess - avgInc * 256;
for (int i = 0; i < 256; i++)
hist[i] += avgInc + (i < remainder ? 1 : 0);
}
// 构建 CDF 映射表
int sum = 0;
for (int i = 0; i < 256; i++)
{
sum += hist[i];
luts[ty, tx, i] = (byte)Math.Clamp(sum * 255 / tilePixels, 0, 255);
}
}
}
// 双线性插值生成结果
var result = new Image<Gray, byte>(width, height);
byte[,,] dstData = result.Data;
for (int y = 0; y < height; y++)
{
for (int x = 0; x < width; x++)
{
// 计算当前像素所在 tile 的中心坐标
double fx = (double)x / actualTileW - 0.5;
double fy = (double)y / actualTileH - 0.5;
int tx0 = Math.Max(0, (int)fx);
int ty0 = Math.Max(0, (int)fy);
int tx1 = Math.Min(tx0 + 1, tilesX - 1);
int ty1 = Math.Min(ty0 + 1, tilesY - 1);
double ax = fx - tx0;
double ay = fy - ty0;
ax = Math.Clamp(ax, 0, 1);
ay = Math.Clamp(ay, 0, 1);
byte val = srcData[y, x, 0];
double v00 = luts[ty0, tx0, val];
double v10 = luts[ty0, tx1, val];
double v01 = luts[ty1, tx0, val];
double v11 = luts[ty1, tx1, val];
double interpolated = v00 * (1 - ax) * (1 - ay) + v10 * ax * (1 - ay)
+ v01 * (1 - ax) * ay + v11 * ax * ay;
dstData[y, x, 0] = (byte)Math.Clamp((int)(interpolated + 0.5), 0, 255);
}
}
_logger.Debug("ApplyCLAHE: ClipLimit={ClipLimit}, TileSize={TileSize}", clipLimit, tileSize);
return result;
}
}