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.