Удивительно, но многие достаточно простые алгоритмы автофокусировки на самом деле неплохо справились с этой задачей. Я реализовал 11 из 16 алгоритмов, описанных в статье Динамическая оценка автофокусировки для автоматизированного микроскопического анализа крови. мазок и мазок Папаниколау Liu, Wang & Sun. Поскольку у меня возникли проблемы с поиском рекомендаций по установке пороговых значений, я также добавил несколько вариантов без порогов. Я также добавил простое, но умное предложение, найденное здесь, на SO: сравните размер файла сжатых изображений JPEG (больший размер = больше деталей = лучший фокус).
Моя процедура автофокуса делает следующее:
- Сохраните 21 изображение с интервалом фокусного расстояния 2 мкм, общий диапазон ±20 мкм.
- Рассчитайте значение фокуса каждого изображения.
- Соответствуйте результату полиному второй степени.
- Найдите позицию, которая дает максимальное значение многочлена.
Все алгоритмы, кроме Histogram Range, дали хорошие результаты. Некоторые алгоритмы немного изменены, например, они используют разницу яркости в обоих направлениях X и Y. Мне также пришлось изменить знак алгоритмов StdevBasedCorrelation, Entropy, ThresholdedContent, ImagePower и ThresholdedImagePower, чтобы получить максимум вместо минимума в положении фокуса. Алгоритмы ожидают 24-битное изображение в градациях серого, где R = G = B. При использовании на цветном изображении будет рассчитываться только синий канал (конечно, легко корректируется).
Оптимальные пороговые значения были найдены путем запуска алгоритмов с пороговыми значениями 0, 8, 16, 24 и т. д. до 255 и выбора наилучшего значения для:
- Стабильное положение фокуса
- Большой отрицательный коэффициент x² приводит к узкому пику в положении фокусировки
- Низкая остаточная сумма квадратов полиномиальной подгонки
Интересно отметить, что алгоритмы ThresholdedSquaredGradient и ThresholdedBrennerGradient имеют почти плоскую линию положения фокуса, коэффициент x² и остаточную сумму квадратов. Они очень нечувствительны к изменению порогового значения.
Реализация алгоритмов:
public unsafe List<Result> CalculateFocusValues(string filename)
{
using(Bitmap bmp = new Bitmap(filename))
{
int width = bmp.Width;
int height = bmp.Height;
int bpp = Bitmap.GetPixelFormatSize(bmp.PixelFormat) / 8;
BitmapData data = bmp.LockBits(new Rectangle(0, 0, width, height), ImageLockMode.ReadOnly, bmp.PixelFormat);
long sum = 0, squaredSum = 0;
int[] histogram = new int[256];
const int absoluteGradientThreshold = 148;
long absoluteGradientSum = 0;
long thresholdedAbsoluteGradientSum = 0;
const int squaredGradientThreshold = 64;
long squaredGradientSum = 0;
long thresholdedSquaredGradientSum = 0;
const int brennerGradientThreshold = 184;
long brennerGradientSum = 0;
long thresholdedBrennerGradientSum = 0;
long autocorrelationSum1 = 0;
long autocorrelationSum2 = 0;
const int contentThreshold = 35;
long thresholdedContentSum = 0;
const int pixelCountThreshold = 76;
long thresholdedPixelCountSum = 0;
const int imagePowerThreshold = 40;
long imagePowerSum = 0;
long thresholdedImagePowerSum = 0;
for(int row = 0; row < height - 1; row++)
{
for(int col = 0; col < width - 1; col++)
{
int current = *((byte *) (data.Scan0 + (row + 0) * data.Stride + (col + 0) * bpp));
int col1 = *((byte *) (data.Scan0 + (row + 0) * data.Stride + (col + 1) * bpp));
int row1 = *((byte *) (data.Scan0 + (row + 1) * data.Stride + (col + 0) * bpp));
int squared = current * current;
sum += current;
squaredSum += squared;
histogram[current]++;
int colDiff1 = col1 - current;
int rowDiff1 = row1 - current;
int absoluteGradient = Math.Abs(colDiff1) + Math.Abs(rowDiff1);
absoluteGradientSum += absoluteGradient;
if(absoluteGradient >= absoluteGradientThreshold)
thresholdedAbsoluteGradientSum += absoluteGradient;
int squaredGradient = colDiff1 * colDiff1 + rowDiff1 * rowDiff1;
squaredGradientSum += squaredGradient;
if(squaredGradient >= squaredGradientThreshold)
thresholdedSquaredGradientSum += squaredGradient;
if(row < bmp.Height - 2 && col < bmp.Width - 2)
{
int col2 = *((byte *) (data.Scan0 + (row + 0) * data.Stride + (col + 2) * bpp));
int row2 = *((byte *) (data.Scan0 + (row + 2) * data.Stride + (col + 0) * bpp));
int colDiff2 = col2 - current;
int rowDiff2 = row2 - current;
int brennerGradient = colDiff2 * colDiff2 + rowDiff2 * rowDiff2;
brennerGradientSum += brennerGradient;
if(brennerGradient >= brennerGradientThreshold)
thresholdedBrennerGradientSum += brennerGradient;
autocorrelationSum1 += current * col1 + current * row1;
autocorrelationSum2 += current * col2 + current * row2;
}
if(current >= contentThreshold)
thresholdedContentSum += current;
if(current <= pixelCountThreshold)
thresholdedPixelCountSum++;
imagePowerSum += squared;
if(current >= imagePowerThreshold)
thresholdedImagePowerSum += current * current;
}
}
bmp.UnlockBits(data);
int pixels = width * height;
double mean = (double) sum / pixels;
double meanDeviationSquared = (double) squaredSum / pixels;
int rangeMin = 0;
while(histogram[rangeMin] == 0)
rangeMin++;
int rangeMax = histogram.Length - 1;
while(histogram[rangeMax] == 0)
rangeMax--;
double entropy = 0.0;
double log2 = Math.Log(2);
for(int i = rangeMin; i <= rangeMax; i++)
{
if(histogram[i] > 0)
{
double p = (double) histogram[i] / pixels;
entropy -= p * Math.Log(p) / log2;
}
}
return new List<Result>()
{
new Result("AbsoluteGradient", absoluteGradientSum),
new Result("ThresholdedAbsoluteGradient", thresholdedAbsoluteGradientSum),
new Result("SquaredGradient", squaredGradientSum),
new Result("ThresholdedSquaredGradient", thresholdedSquaredGradientSum),
new Result("BrennerGradient", brennerGradientSum),
new Result("ThresholdedBrennerGradient", thresholdedBrennerGradientSum),
new Result("Variance", meanDeviationSquared - mean * mean),
new Result("Autocorrelation", autocorrelationSum1 - autocorrelationSum2),
new Result("StdevBasedCorrelation", -(autocorrelationSum1 - pixels * mean * mean)),
new Result("Range", rangeMax - rangeMin),
new Result("Entropy", -entropy),
new Result("ThresholdedContent", -thresholdedContentSum),
new Result("ThresholdedPixelCount", thresholdedPixelCountSum),
new Result("ImagePower", -imagePowerSum),
new Result("ThresholdedImagePower", -thresholdedImagePowerSum),
new Result("JpegSize", new FileInfo(filename).Length),
};
}
}
public class Result
{
public string Algorithm { get; private set; }
public double Value { get; private set; }
public Result(string algorithm, double value)
{
Algorithm = algorithm;
Value = value;
}
}
Чтобы иметь возможность отображать и сравнивать значения фокуса различных алгоритмов, они были масштабированы до значения от 0 до 1 (scaled = (value - min)/(max - min)
).
График всех алгоритмов для диапазона ±20 мкм:
![±20 мкм](https://i.stack.imgur.com/eJw9Q.png)
![0 мкм](https://i.stack.imgur.com/EWbYT.jpg)
![20 мкм](https://i.stack.imgur.com/7hiAh.jpg)
0 µm | 20 µm
Все выглядит очень похоже для диапазона ± 50 мкм:
![±50 мкм](https://i.stack.imgur.com/vkyE2.png)
![0 мкм](https://i.stack.imgur.com/SVrqQ.jpg)
![50 мкм](https://i.stack.imgur.com/jjb46.jpg)
0 µm | 50 µm
При использовании диапазона ±500 мкм все становится интереснее. Четыре алгоритма демонстрируют форму полинома четвертой степени, а остальные начинают больше походить на функции Гаусса. Кроме того, алгоритм диапазона гистограммы начинает работать лучше, чем для меньших диапазонов.
![±500 мкм](https://i.stack.imgur.com/cZytA.png)
![0 мкм](https://i.stack.imgur.com/AOdLI.jpg)
![500 мкм](https://i.stack.imgur.com/PSDQX.jpg)
0 µm | 500 µm
В целом я очень впечатлен производительностью и согласованностью этих простых алгоритмов. Невооруженным глазом довольно сложно сказать, что даже изображение размером 50 мкм не в фокусе, но у алгоритмов нет проблем со сравнением изображений, отстоящих друг от друга всего на несколько микрон.
person
Anlo
schedule
13.03.2013