9

マイクロメートルスケールでポジショニングするためのオートフォーカスルーチンを開発しているので、画像間のフォーカス/ブラーの非常に小さな違いを見つける必要があります。幸いなことに、画像パターンは常に同じです(これらは、元の2 MP画像の256x256の中央の切り抜きです)。

0 µmオフ 50 µmオフ

         Perfect focus         |           50 µm off

上記の2つのより焦点の合った画像を見つけることは問題ではありません、私はほとんどのアルゴリズムがそうするだろうと思います。しかし、以下のように、焦点の差がはるかに少ない画像を比較する必要があります。

5 µmオフ 10 µmオフ

           5 µm off            |           10 µm off

最適な焦点にどんどん近づいていく代わりに、焦点面の反対側で同じ量のぼけがある2つの画像を見つけることもできます。たとえば、-50 µmの画像を保存してから、ブラーが等しい+50 µm付近の画像を見つけようとすることができます。画像が+58µmで見つかったとすると、焦点面は+4 µmに配置する必要があります。

適切なアルゴリズムのアイデアはありますか?

4

3 に答える 3

14

驚くべきことに、非常に単純なオートフォーカス アルゴリズムの多くが、実際にこの問題に対して非常にうまく機能しました。Liu、Wang、Sun による論文「 Dynamic Evaluation of autofocusing for automatic microscopic analysis of Blood smear and pap smear」で概説されている 16 のアルゴリズムのうち 11 を実装しました。しきい値を設定するための推奨事項を見つけるのに苦労したため、しきい値のないバリアントもいくつか追加しました。また、こちらのSOにあるシンプルだが巧妙な提案を追加しました。JPEG圧縮画像のファイルサイズを比較してください(サイズが大きい=詳細=焦点が合っている)。

私のオートフォーカスルーチンは次のことを行います:

  • 焦点距離 2 μm、全範囲 ±20 μm の間隔で 21 枚の画像を保存します。
  • 各画像のフォーカス値を計算します。
  • 結果を 2 次多項式にあてはめます。
  • 多項式の最大値を与える位置を見つけます。

ヒストグラム範囲を除くすべてのアルゴリズムで良好な結果が得られました。一部のアルゴリズムはわずかに変更されています。たとえば、X 方向と Y 方向の両方の明るさの差を使用します。また、StdevBasedCorrelation、Entropy、ThresholdedContent、ImagePower、および ThresholdedImagePower アルゴリズムの符号を変更して、フォーカス位置で最小値ではなく最大値を取得する必要がありました。アルゴリズムは、R = G = B である 24 ビットのグレースケール画像を想定しています。カラー画像で使用すると、青のチャネルのみが計算されます (もちろん、簡単に修正できます)。

最適なしきい値は、しきい値 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 μm の範囲のすべてのアルゴリズムのプロット:

±20μm

0μm 20μm

              0 µm             |             20 µm

±50 µm の範囲では非常によく似ています。

±50μm

0μm 50μm

              0 µm             |             50 µm

±500 µm の範囲を使用すると、さらに興味深い結果が得られます。4 つのアルゴリズムは 4 次多項式の形状を示し、他のアルゴリズムはよりガウス関数のように見え始めます。また、ヒストグラム範囲アルゴリズムは、範囲が狭い場合よりもパフォーマンスが向上し始めます。

±500μm

0μm 500μm

              0 µm             |             500 µm

全体として、これらの単純なアルゴリズムのパフォーマンスと一貫性には非常に感銘を受けました。肉眼では、50 µm の画像でさえ焦点が合っていないと判断するのは非常に困難ですが、アルゴリズムはわずか数ミクロン離れた画像を比較しても問題ありません。

于 2013-03-13T22:20:32.290 に答える
3

元の回答に対するNindzAlのコメントに対する追加の回答:

Extreme Optimizationライブラリを使用して、シャープネス値を 2 次多項式に適合させます。次に、多項式の一次導関数を使用して、最大のシャープネスの距離が抽出されます。

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();
于 2013-07-20T05:24:56.880 に答える
0

無料のライブラリについては、Math.Net がその目的で機能します。

于 2014-08-22T13:19:56.243 に答える