Bandpass filter bank

Question

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;
    }
}

Answer 1

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.

Answer 2

The problem is here:

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;
}

What I don't understand is why should we take the square root, before computing the Math.Pow, especially when we know that the power is an even number. The only thing it does (besides making the code more complex and slow) is to generate NaN for negative values.

I'm not sure if the entire computation is right now, but now all the 12 filtered images appear!

This is used in pre-processing, and is claimed to come from the paper by Kass and Within. I tried to read the original paper, but the quality is very low and hard to read. Do you happen to have a link to a better quality scan of their [15] reference?