Процедура автофокусировки, обнаруживающая очень небольшие различия в размытии

Я разрабатываю процедуру автофокусировки для позиционирования по микрометрической шкале, поэтому мне нужно найти очень небольшие различия в фокусе/размытии между изображениями. К счастью, шаблон изображения всегда будет одним и тем же (это 256x256 центральных кадров исходных 2-мегапиксельных изображений):

0 µm off50 мкм выкл.

         Perfect focus         |           50 µm off

Найти более сфокусированное изображение из двух вышеперечисленных не проблема, я думаю, подойдет большинство алгоритмов. Но мне действительно нужно сравнить изображения с гораздо меньшей разницей в фокусе, как показано ниже:

5 мкм offОтклонение 10 мкм

           5 µm off            |           10 µm off

Альтернативой шагу все ближе и ближе к оптимальной фокусировке является поиск двух изображений с одинаковым уровнем размытия на противоположных сторонах плоскости фокусировки. Например, можно сохранить изображение от -50 мкм, а затем попытаться найти изображение около +50 мкм, где размытие равно. Предположим, что изображение было найдено на +58 мкм, тогда плоскость фокусировки должна располагаться на +4 мкм.

Любые идеи для подходящего алгоритма?


person Anlo    schedule 12.03.2013    source источник
comment
Что происходит, когда вы пробуете базовые алгоритмы? Слишком мало сигнал-шум?   -  person Kendall Frey    schedule 12.03.2013
comment
Ну, я просто предположил, что базовый алгоритм не обнаружит никакой достоверной разницы между изображениями 0, 5 и 10 мкм. Но я только что попробовал несколько и получил довольно многообещающие результаты. Я получу больше изображений с интервалом всего в 1 мкм и посмотрю, удовлетворительны ли результаты.   -  person Anlo    schedule 12.03.2013


Ответы (3)


Удивительно, но многие достаточно простые алгоритмы автофокусировки на самом деле неплохо справились с этой задачей. Я реализовал 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 мкм

0 мкм20 мкм

              0 µm             |             20 µm

Все выглядит очень похоже для диапазона ± 50 мкм:

±50 мкм

0 мкм50 мкм

              0 µm             |             50 µm

При использовании диапазона ±500 мкм все становится интереснее. Четыре алгоритма демонстрируют форму полинома четвертой степени, а остальные начинают больше походить на функции Гаусса. Кроме того, алгоритм диапазона гистограммы начинает работать лучше, чем для меньших диапазонов.

±500 мкм

0 мкм500 мкм

              0 µm             |             500 µm

В целом я очень впечатлен производительностью и согласованностью этих простых алгоритмов. Невооруженным глазом довольно сложно сказать, что даже изображение размером 50 мкм не в фокусе, но у алгоритмов нет проблем со сравнением изображений, отстоящих друг от друга всего на несколько микрон.

person Anlo    schedule 13.03.2013
comment
Привет, я знаю, что это старый вопрос, но я работаю над ОЧЕНЬ похожей проблемой, поэтому я надеялся, что вы сможете уточнить этот вопрос. Подгоните результат к полиному второй степени. - person NindzAI; 18.07.2013
comment
Привет, НиндзАл, см. мой дополнительный ответ (stackoverflow.com/a/17758994/306074) - person Anlo; 20.07.2013
comment
@Anlo Я знаю, что это очень старый вопрос, но, может быть, вы также можете попробовать очень простой алгоритм, основанный на размере JPEG, используемый в NASA Curiosity Mars Rover, см. мой ответ stackoverflow.com/a/32951113/15485 - person Alessandro Jacopson; 05.10.2015
comment
Да, размер Jpeg был одним из алгоритмов, которые я тестировал. Он работал достаточно хорошо. - person Anlo; 06.10.2015

Дополнительный ответ на комментарий Ниндзала к исходному ответу:

Я использую библиотеку Extreme Optimization, чтобы согласовать значения резкости с полиномом второй степени. Затем расстояние максимальной резкости извлекается с использованием первой производной полинома.

Библиотека Extreme Optimization стоит 999 долларов США за одну лицензию разработчика, но я уверен, что существуют математические библиотеки с открытым исходным кодом, которые могут выполнить настройку так же хорошо.

// Distances (in µm) where the images were saved
double[] distance = new double[]
{
  -50,
  -40,
  -30,
  -20,
  -10,
    0,
  +10,
  +20,
  +30,
  +40,
  +50,
};

// Sharpness value of each image, as returned by CalculateFocusValues()
double[] sharpness = new double[]
{
  3960.9,
  4065.5,
  4173.0,
  4256.1,
  4317.6,
  4345.4,
  4339.9,
  4301.4,
  4230.0,
  4131.1,
  4035.0,
};

// Fit data to y = ax² + bx + c (second degree polynomial) using the Extreme Optimization library
const int SecondDegreePolynomial = 2;
Extreme.Mathematics.Curves.LinearCurveFitter fitter = new Extreme.Mathematics.Curves.LinearCurveFitter();
fitter.Curve = new Extreme.Mathematics.Curves.Polynomial(SecondDegreePolynomial);
fitter.XValues = new Extreme.Mathematics.LinearAlgebra.GeneralVector(distance,  true);
fitter.YValues = new Extreme.Mathematics.LinearAlgebra.GeneralVector(sharpness, true);
fitter.Fit();

// Find distance of maximum sharpness using the first derivative of the polynomial
// Using the sample data above, the focus point is located at distance +2.979 µm
double focusPoint = fitter.Curve.GetDerivative().FindRoots().First();
person Anlo    schedule 20.07.2013

Что касается бесплатной библиотеки, для этой цели подойдет Math.Net.

person KorsaR    schedule 22.08.2014