带通滤波器组

我已经实现了本文中描述的一组定向带通滤波器。

请参阅“ 2.1预处理 ”一节中的最后一段。

我们选择了12个不重叠的滤镜，分析了12个不同的方向，相对于15°相互旋转。

我有以下问题，

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

---

源代码：

以下是完整的VS2013解决方案作为压缩文件。

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

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 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 Apply(Bitmap bitmap) { Kernels = new List(); 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 list = new List(); 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; } } 

Answer 1

正如我在前面的评论中指出的那样，大多数滤波器输出都是空白的，因为它们包含NaN 。 这些是由参考文章中的等式（1）和（2）的实现引起的。 与原作者取得联系可能最有可能再现原始结果，但至少可以确保不会产生任何NaN ：

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

另一方面，考虑到方程的一般forms是巴特沃斯滤波器 （同时提到带通滤波），以及看似不必要的平方根然后取幂（这表明错过了明显的简化，或更可能在我认为在渲染方程时出现错误），我建议改为使用以下等式：

哪里 是图像的中心。 这可以用以下方式实现：

 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期望滤波器内核在空间域中提供（尽管计算的核心是在频域中完成的）。 应用这些filter的最简单方法（并且仍在空间域中获得正确的填充）是：

• 将频域内核大小设置为填充大小以计算频域中的内核
• 将频域内核转换为空间域
• 将已添加填充的内核清零
• 将内核转换回频域

这可以通过以下KassWitkinKernel.Pad修改版本来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 

在计算卷积时，您将在卷积中跳过内核的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); } 

所以这应该让你超越你的问题中指出的凹凸。 当然，如果打算重现论文的结果，你还有其他障碍可以克服。 第一个是实际使用锐化图像作为滤波器组的输入。 执行此操作时，您可能需要首先平滑图像的边缘，以避免在图像周围生成强边缘，这会使线检测算法的结果发生偏差。

Answer 2

问题出在这里：

 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); if (!double.IsNaN(sqrt) && Math.Abs(pow - Math.Pow(uStarDu + vStarDv, N)) > 1e-7) { //execution will never reach here!! } pow = Math.Pow(uStarDu + vStarDv, N); double div = 1 + 0.414 * pow; double returns = 1 / div; return returns; } 

我不明白的是为什么在计算Math.Pow之前我们应该采用平方根，特别是当我们知道幂是偶数时。 它唯一能做的事情（除了使代码更复杂和更慢）是为负值生成NaN。

我不确定整个计算是否正确，但现在所有12个过滤后的图像都出现了！

这用于预处理，据称来自Kass和Within的论文。 我试着阅读原始论文，但质量非常低，难以阅读。 您是否碰巧有更好质量的[15]参考扫描链接？