c# - 带通滤波器组

我已经实现了一组定向带通滤波器 described in this article .

请参阅“2.1 预处理”部分的最后一段。

We selected 12 not overlapping filters, to analyze 12 different directions, rotated with respect to 15° each other.

我遇到以下问题，

滤波器组应该生成 12 个滤波后的图像。但是，实际上，我只有 03 个输出，如下面的快照所示，

源代码:

Here is the complete VS2013 solution as a zipped file.

这是源代码中最相关的部分，

public class KassWitkinFunction
{
    /*
     *  tx = centerX * cos 
     *  ty = centerY * sin
     *  
     *  u* =   cos . (u + tx) + sin . (v + ty)
     *  v* = - sin . (u + tx) + cos . (v + ty)
     *  
     */
    //#region MyRegion
    public static double tx(int centerX, double theta)
    {
        double costheta = Math.Cos(theta);
        double txx = centerX * costheta;
        return txx;
    }

    public static double ty(int centerY, double theta)
    {
        double sintheta = Math.Sin(theta);
        double tyy = centerY * sintheta;
        return tyy;
    }

    public static double uStar(double u, double v, int centerX, int centerY, double theta)
    {
        double txx = tx(centerX, theta);
        double tyy = ty(centerY, theta);
        double sintheta = Math.Sin(theta);
        double costheta = Math.Cos(theta);

        double cosThetaUTx = costheta * (u + txx);
        double sinThetaVTy = sintheta * (v + tyy);

        double returns = cosThetaUTx + sinThetaVTy;

        return returns;
    }
    //#endregion

    public static double vStar(double u, double v, int centerX, int centerY, double theta)
    {
        double txx = tx(centerX, theta);
        double tyy = ty(centerY, theta);
        double sintheta = Math.Sin(theta);
        double costheta = Math.Cos(theta);

        double sinThetaUTx = (-1) * sintheta * (u + txx);
        double cosThetaVTy = costheta * (v + tyy);

        double returns = sinThetaUTx + cosThetaVTy;

        return returns;
    }

    public static double ApplyFilterOnOneCoord(double u, double v, double Du, double Dv, int CenterX, int CenterY, double Theta, int N)
    {
        double uStar = KassWitkinFunction.uStar(u, v, CenterX, CenterY, Theta);
        double vStar = KassWitkinFunction.vStar(u, v, CenterX, CenterY, Theta);

        double uStarDu = uStar / Du;
        double vStarDv = vStar / Dv;

        double sqrt = Math.Sqrt(uStarDu + vStarDv);
        double _2n = 2 * N;
        double pow = Math.Pow(sqrt, _2n);
        double div = 1 + 0.414 * pow;

        double returns = 1/div;

        return returns;
    }
}

public class KassWitkinKernel
{
    public readonly int N = 4;
    public int Width { get; set; }
    public int Height { get; set; }
    public double[,] Kernel { get; private set; }
    public double[,] PaddedKernel { get; private set; }
    public double Du { get; set; }
    public double Dv { get; set; }
    public int CenterX { get; set; }
    public int CenterY { get; set; }
    public double ThetaInRadian { get; set; }

    public void SetKernel(double[,] value)
    {
        Kernel = value;
        Width = Kernel.GetLength(0);
        Height = Kernel.GetLength(1);
    }

    public void Pad(int newWidth, int newHeight)
    {
        double[,] temp = (double[,])Kernel.Clone();
        PaddedKernel = ImagePadder.Pad(temp, newWidth, newHeight);
    }

    public Bitmap ToBitmap()
    {
        return ImageDataConverter.ToBitmap(Kernel);
    }

    public Bitmap ToBitmapPadded()
    {
        return ImageDataConverter.ToBitmap(PaddedKernel);
    }

    public Complex[,] ToComplex()
    {
        return ImageDataConverter.ToComplex(Kernel);
    }

    public Complex[,] ToComplexPadded()
    {
        return ImageDataConverter.ToComplex(PaddedKernel);
    }

    public void Compute()
    {
        Kernel = new double[Width, Height];

        for (int i = 0; i < Width; i++)
        {
            for (int j = 0; j < Height; j++)
            {
                Kernel[i, j] = (double)KassWitkinFunction.ApplyFilterOnOneCoord(i, j,
                                                                            Du,
                                                                            Dv,
                                                                            CenterX,
                                                                            CenterY,
                                                                            ThetaInRadian,
                                                                            N);

                //Data[i, j] = r.NextDouble();
            }
        }

        string str = string.Empty;
    }
}

public class KassWitkinBandpassFilter
{
    public Bitmap Apply(Bitmap image, KassWitkinKernel kernel)
    {
        Complex[,] cImagePadded = ImageDataConverter.ToComplex(image);
        Complex[,] cKernelPadded = kernel.ToComplexPadded();
        Complex[,] convolved = Convolution.Convolve(cImagePadded, cKernelPadded);

        return ImageDataConverter.ToBitmap(convolved);
    }
}

public class KassWitkinFilterBank
{
    private List<KassWitkinKernel> Kernels;
    public int NoOfFilters { get; set; }
    public double FilterAngle { get; set; }
    public int WidthWithPadding { get; set; }
    public int HeightWithPadding { get; set; }
    public int KernelDimension { get; set; }

    public KassWitkinFilterBank()
    {}

    public List<Bitmap> Apply(Bitmap bitmap)
    {
        Kernels = new List<KassWitkinKernel>();

        double degrees = FilterAngle;

        KassWitkinKernel kernel;
        for (int i = 0; i < NoOfFilters; i++)
        {
            kernel = new KassWitkinKernel();
            kernel.Width = KernelDimension;
            kernel.Height = KernelDimension;
            kernel.CenterX = (kernel.Width) / 2;
            kernel.CenterY = (kernel.Height) / 2;
            kernel.Du = 2;
            kernel.Dv = 2;
            kernel.ThetaInRadian = Tools.DegreeToRadian(degrees);
            kernel.Compute();
            kernel.Pad(WidthWithPadding, HeightWithPadding);

            Kernels.Add(kernel);

            degrees += degrees;
        }

        List<Bitmap> list = new List<Bitmap>();

        foreach (KassWitkinKernel k in Kernels)
        {
            Bitmap image = (Bitmap)bitmap.Clone();

            Complex[,] cImagePadded = ImageDataConverter.ToComplex(image);
            Complex[,] cKernelPadded = k.ToComplexPadded();
            Complex[,] convolved = Convolution.Convolve(cImagePadded, cKernelPadded);

            Bitmap temp = ImageDataConverter.ToBitmap(convolved);

            list.Add(temp);
        }

        return list;
    }
}

最佳答案

正如我之前在评论中指出的，大多数过滤器输出都是空白，因为它们包含 NaN。这些都是由方程(1)和(2)的实现 your reference article 。与原始作者取得联系可能有最好的机会重现原始结果，但至少您可以确保不会产生 NaN:

double arg = uStarDu + vStarDv;
double div = 1 + 0.414 * Math.Pow(Math.Abs(arg), N);

另一方面，考虑到方程的一般形式让人想起 Butterworth filter (连同提到的带通滤波)，以及看似不必要的平方根后跟幂(这表明要么错过了明显的简化，要么在我看来更有可能是一个错误在方程的渲染中)，我建议改用以下方程:

$\begin{align} u' &= (u - C_x) \cos\theta + (v - C_y) \sin\theta &, 0 \le u \lt \mbox{width} \ v' &= - (u - C_x) \sin\theta + (v - C_y) \cos\theta &, 0 \le v \lt \mbox{height} \ H(u,v) &= \frac{1}{\sqrt{1 + \left(\frac{u'}{D_u} + \frac{v'}{D_v}\right)^{2N}}}\end{align}$

哪里是图像的中心。这可以通过以下方式实现:

public static double uStar(double u, double v, int centerX, int centerY, double theta)
{
    double sintheta = Math.Sin(theta);
    double costheta = Math.Cos(theta);

    return costheta * (u - centerX) + sintheta * (v - centerY);
}

public static double vStar(double u, double v, int centerX, int centerY, double theta)
{
    double sintheta = Math.Sin(theta);
    double costheta = Math.Cos(theta);

    return (-1) * sintheta * (u - centerX) + costheta * (v - centerY);
}

public static double ApplyFilterOnOneCoord(double u, double v, double Du, double Dv, int CenterX, int CenterY, double Theta, int N)
{
    double uStarDu = KassWitkinFunction.uStar(u, v, CenterX, CenterY, Theta) / Du;
    double vStarDv = KassWitkinFunction.vStar(u, v, CenterX, CenterY, Theta) / Dv;

    double arg = uStarDu + vStarDv;
    double div = Math.Sqrt(1 + Math.Pow(arg, 2*N));;

    return 1/div;
}

现在您必须意识到这些方程是针对频域中的滤波器表示给出的，而您的Convolution.Convolve 期望在空间域中提供滤波器内核(尽管计算的核心是在频域中完成的)。应用这些过滤器的最简单方法(并且仍然在空间域中获得正确的填充)是:

将频域内核大小设置为填充大小以计算频域内核
将频域核变换到空间域
将要添加填充的内核清零
将内核变换回频域

这可以通过 KassWitkinKernel.Pad 的以下修改版本来实现:

private Complex[,] cPaddedKernel;

public void Pad(int unpaddedWidth, int unpaddedHeight, int newWidth, int newHeight)
{
  Complex[,] unpaddedKernelFrequencyDomain = ImageDataConverter.ToComplex((double[,])Kernel.Clone());
  FourierTransform ftInverse = new FourierTransform();
  ftInverse.InverseFFT(FourierShifter.RemoveFFTShift(unpaddedKernelFrequencyDomain));

  Complex[,] cKernel = FourierShifter.FFTShift(ftInverse.GrayscaleImageComplex);

  int startPointX = (int)Math.Ceiling((double)(newWidth - unpaddedWidth) / (double)2) - 1;
  int startPointY = (int)Math.Ceiling((double)(newHeight - unpaddedHeight) / (double)2) - 1;
  for (int j = 0; j < newHeight; j++)
  {
    for (int i=0; i<startPointX; i++)
    {
      cKernel[i, j] = 0;
    }
    for (int i = startPointX + unpaddedWidth; i < newWidth; i++)
    {
      cKernel[i, j] = 0;
    }
  }
  for (int i = startPointX; i < startPointX + unpaddedWidth; i++)
  {
    for (int j = 0; j < startPointY; j++)
    {
      cKernel[i, j] = 0;
    }
    for (int j = startPointY + unpaddedHeight; j < newHeight; j++)
    {
      cKernel[i, j] = 0;
    }
  }

  FourierTransform ftForward = new FourierTransform(cKernel); ftForward.ForwardFFT();
  cPaddedKernel = ftForward.FourierImageComplex;
}

public Complex[,] ToComplexPadded()
{
  return cPaddedKernel;
}

稍后在计算卷积时，您将跳过卷积中内核的 FFT。请注意，您可以类似地避免为滤波器组中的每个滤波器重新计算图像的 FFT。如果预先计算图像的 FFT，则获得卷积所需的剩余计算简化为频域乘法和最终的逆变换:

public partial class Convolution
{
  public static Complex[,] ConvolveInFrequencyDomain(Complex[,] fftImage, Complex[,] fftKernel)
  {
    Complex[,] convolve = null;

    int imageWidth = fftImage.GetLength(0);
    int imageHeight = fftImage.GetLength(1);

    int maskWidth = fftKernel.GetLength(0);
    int maskHeight = fftKernel.GetLength(1);

    if (imageWidth == maskWidth && imageHeight == maskHeight)
    {
      Complex[,] fftConvolved = new Complex[imageWidth, imageHeight];

      for (int j = 0; j < imageHeight; j++)
      {
        for (int i = 0; i < imageWidth; i++)
        {
          fftConvolved[i, j] = fftImage[i, j] * fftKernel[i, j];
        }
      }

      FourierTransform ftForConv = new FourierTransform();
      ftForConv.InverseFFT(fftConvolved);
      convolve = FourierShifter.FFTShift(ftForConv.GrayscaleImageComplex);

      Rescale(convolve);
    }
    else
    {
      throw new Exception("padding needed");
    }

    return convolve;
  }
}

它将在 KassWitkinFilterBank.Apply 中使用，如下所示:

Bitmap image = (Bitmap)bitmap.Clone();

Complex[,] cImagePadded = ImageDataConverter.ToComplex(image);
FourierTransform ftForImage = new FourierTransform(cImagePadded); ftForImage.ForwardFFT();
Complex[,] fftImage = ftForImage.FourierImageComplex;

foreach (KassWitkinKernel k in Kernels)
{
  Complex[,] cKernelPadded = k.ToComplexPadded();
  Complex[,] convolved = Convolution.ConvolveInFrequencyDomain(fftImage, cKernelPadded);

  Bitmap temp = ImageDataConverter.ToBitmap(convolved);
  list.Add(temp);
}

因此，这应该可以帮助您克服问题中指出的障碍。当然，如果目的是重现论文的结果，您仍然需要克服其他障碍。第一个是实际使用锐化图像作为滤波器组的输入。执行此操作时，您可能需要首先平滑图像的边缘，以避免生成强边缘图像周围，这会扭曲线条检测算法的结果。

关于c# - 带通滤波器组，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/39209951/

c# - 带通滤波器组

上一篇：c# - 是什么允许 .NET 任务在等待 I/O 完成后继续执行？

下一篇：c# - 使用 lambda 表达式比较两个列表