带通滤波器组 [英] Bandpass filter bank
问题描述
我已经实现了一个定向的带通滤波器库
源代码：
其中是图像的中心。这可以通过以下方式来实现：
pre $
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））;;
返回1 / div;
$ b现在你必须意识到这些方程是给出频率中的滤波器表示而您的 Convolution.Convolve
期望在空间域中提供滤波器内核（尽管计算的核心是在频域中完成的）。
应用这些过滤器的最简单的方法（并且仍然在空间域中获得正确的填充）是：
- 设置频域内核大小到填充大小以计算频域中的内核
- 将频域内核转换到空间域
- 在那里填充已经被添加
- 将内核转换回频域
使用以下修改版本的 KassWitkinKernel.Pad ：
private Complex [，] cPaddedKernel;
$ b public void Pad（int unpaddedWidth，int unpaddedHeight，int newWidth，int newHeight）
{
Complex [，] unpaddedKernelFrequencyDomain = ImageDataConverter.ToComplex（（double [，]）Kernel 。克隆（））;
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; $（b）b int startPointY =（int）Math.Ceiling（（double）（newHeight - unpaddedHeight）/（double）2） - 1;
for（int j = 0; j
for（int i = 0; i
cKernel [i，j] = 0;
for（int i = startPointX + unpaddedWidth; i
cKernel [i，j] = 0; （int i = startPointX; i
}
for（int j = 0; j
cKernel [i，j] = 0;
for（int j = startPointY + unpaddedHeight; j
cKernel [i，j] = 0;
}
}
FourierTransform ftForward = new FourierTransform（cKernel）; ftForward.ForwardFFT（）;
cPaddedKernel = ftForward.FourierImageComplex;
}
public Complex [，] ToComplexPadded（）
{
return cPaddedKernel;
$ b 当计算卷积时，你会忽略卷积内核的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]; （int j = 0; j
（int i = 0; i
fftConvolved [i，j] = fftImage [i，j] * fftKernel [i，j];
}
}
FourierTransform ftForConv = new FourierTransform（）;
ftForConv.InverseFFT（fftConvolved）;
convolve = FourierShifter.FFTShift（ftForConv.GrayscaleImageComplex）;
重新标定（卷积）;
}
else
{
throw new Exception（padding needed）;
}
return convolve;
$ b $ KBSWitkinFilterBank.Apply 如下：
位图图像=（位图）bitmap.Clone（）;
Complex [，] cImagePadded = ImageDataConverter.ToComplex（image）;
FourierTransform ftForImage = new FourierTransform（cImagePadded）; ftForImage.ForwardFFT（）;
Complex [，] fftImage = ftForImage.FourierImageComplex;
foreach（Kernel中的KassWitkinKernel k）
{
Complex [，] cKernelPadded = k.ToComplexPadded（）;
Complex [，] convolved = Convolution.ConvolveInFrequencyDomain（fftImage，cKernelPadded）;
位图temp = ImageDataConverter.ToBitmap（convolved）;
list.Add（temp）;
$ b $ p
$ b 所以这应该让你通过在你的问题中指出的凹凸。
当然，如果你想重现论文的结果，你还有其他障碍要克服。
第一个实际使用锐化的图像作为输入到滤波器组。
当你这样做的时候，你可能需要首先平滑图像的边缘，以避免在图像周围产生一个强大的边缘
，这会歪曲线条检测算法的结果。
I have implemented a bank of oriented bandpass filters described in this article.
See the last paragraph of the section named "2.1 Pre-processing".
We selected 12 not overlapping filters, to analyze 12 different directions, rotated with respect to 15° each other.
I am having the following issue,
The filter-bank was supposed to generate 12 filtered images. But, in reality, I am having only 03 outputs as seen in the following snapshot,
Source Code:
Here is the complete VS2013 solution as a zipped file.
Here is the most relevant part of the source code,
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;
}
}
解决方案 As I pointed out earlier in comments, most of the filter outputs are blank because they contain NaNs. These are caused by
the implementation of equations (1) and (2) from
your reference article.
Getting in touch with the original authors would probably have the best chance of reproducing the original results, but at the very least you could ensure that no NaNs are produced with:
double arg = uStarDu + vStarDv;
double div = 1 + 0.414 * Math.Pow(Math.Abs(arg), N);
On the other hand, given the general form of the equation being reminescent of a Butterworth filter
(together with the mention of bandpass filtering),
and the seemingly unecessary square root followed by exponentiation (which suggest either a missed obvious simplification, or more likely in my opinion an error
in rendering of the equation), I would suggest to instead use the following equation:
where is the center of the image. This could be implemented with:
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;
}
Now you must realize that those equations are given for the filter representation in the frequency domain, whereas your Convolution.Convolve
expect the filter kernel to be provided in the spatial domain (despite the core of the computation being done in the frequency domain).
The easiest way to apply these filters (and still get the correct padding in the spatial domain) is to:
- set the frequency domain kernel size to the padded size to compute the kernel in the frequency domain
- transform the frequency domain kernel to spatial domain
- zero out the kernel where the padding would have been added
- transform the kernel back to frequency domain
This can be achieved with the following modified version of 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;
}
Latter when computing the convolution you would skip the FFT for the kernel in the convolution.
Note that you could similarly avoid recomputing the image's FFT for every filter in the filter bank.
If you precompute the FFT of the image, the remaining computations required to get the convolution
is reduced to the frequency domain multiplication and the final inverse transform:
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;
}
}
Which would be used in KassWitkinFilterBank.Apply as follow:
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);
}
So that should get you past the bump indicated in your question.
Of course if the intention is to reproduce the results of the paper you still have other hurdles to get over.
The first one being to actually use a sharpened image as input to the filter bank.
While you do this, you may want to first smooth the edges of the image to avoid generating a strong edge all
around the image, which would skew the results of the line detection algorithm.
这篇关于带通滤波器组的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
