FFT Convolution

I have written the following routines to convolve two images in the frequency domain which are represented as 2d Complex arrays.

How can I optimize my routines for better performance?

public static class Convolution
{
public static Complex[,] Convolve(Complex[,] image, Complex[,] mask)
{
Complex[,] convolve = null;

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

{
FourierTransform ftForImage = new FourierTransform(image); ftForImage.ForwardFFT();

Complex[,] fftImage = ftForImage.FourierImageComplex;

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 = ftForConv.GrayscaleImageComplex;

Rescale(convolve);

convolve = FourierShifter.FFTShift(convolve);
}
else
{
}

return convolve;
}

//Rescale values between 0 and 255.
private static void Rescale(Complex[,] convolve)
{
int imageWidth = convolve.GetLength(0);
int imageHeight = convolve.GetLength(1);

double maxAmp = 0.0;
for (int j = 0; j < imageHeight; j++)
{
for (int i = 0; i < imageWidth; i++)
{
maxAmp = Math.Max(maxAmp, convolve[i, j].Magnitude);
}
}
double scale = 255.0 / maxAmp;

for (int j = 0; j < imageHeight; j++)
{
for (int i = 0; i < imageWidth; i++)
{
convolve[i, j] = new Complex(convolve[i, j].Real * scale, convolve[i, j].Imaginary * scale);
maxAmp = Math.Max(maxAmp, convolve[i, j].Magnitude);
}
}
}
}


A few detail remarks:

Complex[,] convolve = null;


This variable declaration should be moved further down, to the line where it is actually needed. Initializing it with null is misleading.

maxAmp = Math.Max(maxAmp, convolve[i, j].Magnitude);


Since Magnitude involves a square root, it is expensive to calculate. Prefer calculating the square magnitude only and calculate the square root only once at the end. This saves width * height - 1 square root calculations.

In the second loop, maxAmp is calculated again but never used. Removing it saves another width * height calculations.

Have you tried to use anything from the TPL yet?

You could try to parallelize one of the loops and see if the performance increases or drops: see In a nested loop, should Parallel.For be used on the outer or inner loop?

Question Review

• If you work with types that are not included in the standard .NET Framework, include the library and namespace for these types in the question: FourierTransform, FourierShifter, Complex.

Code Review

• Check arguments against null to avoid the nasty NullReferenceException.
• Give meaningful names to variables. Do not use ft and fft prefixes. This kills readability.
• Exit early on invalid state of the arguments. This simplifies the number of nested code blocks. Invert if (imageWidth == maskWidth && imageHeight == maskHeight) and throw error.
• Don't put multiple statements on a single line; FourierTransform ftForMask = new FourierTransform(mask); ftForMask.ForwardFFT();.
• Don't introduce unnecessary white lines between simple statements.
• Prefer the use of var for declaring variables.
• I presume maskeHeight has a typo in it -> maskHeight
• I don't see much room for optimizing the performance, since each step requires intermediate results from the previous one; perhaps as suggested in another answer, making one of the loops use Parallel, would benefit you.

If you use a certain pattern time and again, you might want to make a utility method for it. As suggested by Roland in the comments, we could implement a possible performance gain here. I would opt to call AsParallel() on the items. You could also try calling Parallel.For instead.

for (int j = 0; j < imageHeight; j++)
{
for (int i = 0; i < imageWidth; i++)
{
// do something ..
}
}

static void Walk(int height, int width, Action<int, int> visit)
{
foreach (var point in (
from j in Enumerable.Range(0, height)
from i in Enumerable.Range(0, width)
select (i, j)).AsParallel())
{
visit(point.i, point.j);
}
}


Method 'Convolve' refactored

public static Complex[,] Convolve(Complex[,] image, Complex[,] mask)
{
image = image ?? throw new ArgumentNullException(nameof(image));

var imageWidth = image.GetLength(0);
var imageHeight = image.GetLength(1);

{
}

var imageTransform = new FourierTransform(image);

imageTransform.ForwardFFT();

var imageComplex = imageTransform.FourierImageComplex;
var convolvedComplex = new Complex[imageWidth, imageHeight];

Walk(imageHeight, imageWidth,
(i, j) => convolvedComplex[i, j] = imageComplex[i, j] * maskComplex[i, j]);

var convolvedTransform = new FourierTransform();
convolvedTransform.InverseFFT(convolvedComplex);
var convolve = convolvedTransform.GrayscaleImageComplex;
Rescale(convolve);
convolve = FourierShifter.FFTShift(convolve);

return convolve;
}

• Why add explicit code for NullArgumentException if it would be thrown anyway a few lines further down, unconditionally, easy to relate to the parameter? – Roland Illig Jul 23 '19 at 5:45
• Does your Walk function have any chance of being faster than a simple for loop? The OP asked about making the code a bit faster. – Roland Illig Jul 23 '19 at 5:48
• @RolandIllig (1) a few lines further a NullReferenceException would be thrown, which is harder to introspect than a NullArgumentException. It has a much more generic error message and less accurate stack trace. – dfhwze Jul 23 '19 at 5:52
• @RolandIllig (2) valid point, the Parallel suggested by t3chb0t could be implemented here. – dfhwze Jul 23 '19 at 5:52