Previously, I wrote about performing convolution in SkiaSharp. I used the **CreateMatrixConvolution** method from the **SKImageFilter** class to convolve kernels with the source image. This method allows you to specify kernels of arbitrary size, allows you to specify how the edge pixels in the image are handled, and is lightning fast. I was particularly impressed with the execution speed of the method when considering the size of the source image (4032x3024), and the amount of processing that’s performed when convolving an image with a kernel, particularly for larger kernels.

However, convolution is still considered relatively basic in image processing terms, despite the amount of processing performed. Therefore, I decided to move things up a gear and see what the execution speed would be like when performing frequency filtering, which substantially increases the amount of processing performed on an image.

In this blog post I’ll discuss implementing frequency filtering in SkiaSharp. The sample this code comes from can be found on GitHub.

### Implementing frequency filtering

As any computer science undergrad can tell you, a Fourier transform take a signal from the time domain and transforms it into the frequency domain. What this means is it takes a signal, and decomposes it into its constituent frequencies. The world of signal processing is full of applications for the Fourier transform, with one of the stanard applications being musical instrument tuners, which compare the frequency of sound input to the known frequencies of specific notes. Fourier transforms can be performed using analog electronics, digital electronics, and even at the speed of light using optics. For more info about Fourier transforms, see Fourier transform.

Another standard application of the Fourier transform is to perform frequency filtering. A signal can be transformed to the frequency domain, where specific frequencies, or frequency ranges, can be removed. This enables a class of filters to be implemented, known as pass through filters: specifically high-pass, low-pass, and band-pass filters. In terms of image processing, a low-pass filter attentuates high frequencies and retains low frequencies unchanged, with the result being similar to that of a smoothing filter. A high-pass filter attentuates low frequencies, and retains high frequencies unchanged, with the results being similar to that of edge detection. A band-pass filter attentuates very low and very high frequencies, but retains a middle band of frequencies, and is often used to enhance edges while reducing noise.

In terms of image processing, the Fourier transform actually transforms an image from the spatial domain to the frequency domain, where frequency filtering can then be performed. The process of performing frequency filtering on an image is:

- Fourier transform the image.
- Filter the image frequencies in the frequency domain.
- Inverse Fourier transform the frequency data back into the spatial domain.

In reality, the process of Fourier transforming an image is more complex than the bullet points suggest. The frequency data that results from a Fourier transform is expressed as a complex number (a number that has a real, and imaginary component), whose magnitude represents the amount of that frequency present in the original signal, with the phase offset being the basic sinusoid in that frequency. For more information about complex numbers, see complex numbers.

In addition, as well as handling complex numbers, there are other issues to be dealt with when Fourier transforming an image. Firstly, computing a Fourier transform is too slow to be practical, because it has a complexity of O(N^{2}). However, this issue can be overcome by implementing a Fast Fourier Transform (FFT), which has a complexity of O(N log N). For more information about this algorithm, see Fast Fourier transform. The second issue is that Fourier transforms only operate on image dimensions that are a power of 2 (512x512, 1024x1024 etc.). While there are techniques to allow arbitrarily dimensoined images to be Fourier transformed, they are beyond the scope of this blog post. Therefore, any images processed by the algorithm must first be manually resized. The final issue is which image data do you Fourier transform? The answer is to Fourier transform a greyscale representation of the image, as it’s the intensity data that caries more information than colour channels.

For clarity, restating all this leads to the following assumptions for my implementation:

- The input image for the Fourier transform will be a 32-bit image, with RGBA channels.
- The input image dimensions must be a power of 2. This is due to the FFT algorithm (Radix-2 FFT) being used.
- The input image will first be converted to a greyscale representation.
- The greyscale representation will then be converted to a complex number representation.
- The complex number data will undergo a 2D FFT.
- Filtering will then be performed in the frequency domain.
- The filtered complex data will be inverse FFT’d.
- The filtered complex data will be converted back from complex numbers to pixel data for display.
- The output image will be in greyscale, but still stored as a RGBA image.

### Implementation

It takes more code to implement a 2D FFT than can be covered in a blog post, particularly as there are complexities I haven’t outlined here. Rather than show every last bit of code, I’ll just state that there’s a **Complex** type that represents a complex number, and a **ComplexImage** type that represents an image as a 2D array of complex numbers. The full source code can be found on GitHub.

The following code example shows a high level overview of how the frequency filtering process is implemented:

```
public static unsafe SKPixmap FrequencyFilter(this SKImage image, int min, int max)
{
ComplexImage complexImage = image.ToComplexImage();
complexImage.FastFourierTransform();
FrequencyFilter filter = new FrequencyFilter(new FrequencyRange(min, max));
filter.Apply(complexImage);
complexImage.ReverseFastFourierTransform();
SKPixmap pixmap = complexImage.ToSKPixmap(image);
return pixmap;
}
```

The **FrequencyFilter** method converts the image to a **ComplexImage** object, which then undergoes a 2D FFT. A **FrequencyFilter** object is then created, based on the values of the **Slider** objects displayed on the UI.The filter is applied to the **ComplexImage** object, and then inverse FFT’d, before being converted back to pixel data.

A Fourier transform produces a complex number valued output image that can be displayed as two images, for the real and imaginary coefficients, or as their magnitude and phase. In image processing, it’s usually the magnitude of the Fourier transform that’s displayed, as it contains most of the information about the geometric structure of the source image. However, to inverse transform the data after processing in the frequency domain, both the magnitude and phase of the Fourier data is required, and so must be preserved.

It’s possible to view the Fourier transformed image (known as a frequency spectra) by commenting out a couple of lines in the **FrequencyFilter** method shown above. I mentioned earlier that are additional complexities when implementing a FFT, and one of them is that the dynamic range of Fourier coefficients is too large to be displayed in an image, and the resulting image would appear all black. However, if a logarithmic transformation is applied to the coefficients (which the source code does), the Fourier transformed image can be displayed:

Frequency filtering is performed by the **Apply** method in the **FrequencyFilter** class:

```
public void Apply(ComplexImage complexImage)
{
if (!complexImage.IsFourierTransformed)
{
throw new ArgumentException("The source complex image should be Fourier transformed.");
}
int width = complexImage.Width;
int height = complexImage.Height;
int halfWidth = width >> 1;
int halfHeight = height >> 1;
int min = frequencyRange.Min;
int max = frequencyRange.Max;
Complex[,] data = complexImage.Data;
for (int i = 0; i < height; i++)
{
int y = i - halfHeight;
for (int j = 0; j < width; j++)
{
int x = j - halfWidth;
int d = (int)Math.Sqrt(x * x + y * y);
if ((d > max) || (d < min))
{
data[i, j].Re = 0;
data[i, j].Im = 0;
}
}
}
}
```

This method iterates over the complex image data, and zeros the real and imaginery values that lie outside the frequency range specified by the **min** and **max** values. Conversely, frequency data within the **min** and **max** values is passed through. This method therefore implements a band-pass filter, which can be configured to operate at any frequency range.

It should therefore follow from this, that if a frequency filter with a min value of 0 and a max value of 1024 is applied, the resulting inverse transformed frequency filtered image should be a perfect greyscale representation of the original source image. The following screenshot shows this:

Furthermore, because the earlier frequency spectra shows that the image is largely comprised of low frequency data, a frequency filter with a min value of 0 and a max value of 128 still results (after inverse FFT) in a perfect greyscale representation of the original source image. The following screenshot shows this:However, a frequency filter with a min value of 10 and a max value of 128 yields the following image:In this image, because some of the low frequency data has been removed, only sharp changes in intensity values are being preserved, with the resulting image beginning to look as if it’s been edge detected. Similarly, a frequency filter with a min value of 20 and a max value of 128 furthers this effect:Again, the output is now looking even more like the output of an edge detector.While this example is of no immediate practical use, it hopefully shows the potential of what can be achieved with frequency filtering. One of the main uses of frequency filtering in image processing is to remove noise. If the frequency range of the noise can be identified, which it often can be, that frequency range can be removed, resulting in a denoised image. Another real world application of the Fourier transform in imaging, is producing images of steel rebar quality inside concrete (think steel rebar inside concrete walls, bridges etc.). In this case, the transducer data can be deconvolved (in the frequency domain) with the point spread function of the transducer, to yield images of steel rebar quality inside concrete, from which deterioration can be identified.

### Wrapping up

What I set out to address here is, is the combination of Xamarin.Forms and SkiaSharp a viable platform for writing cross-platform imaging apps, when performing substantial image processing? My main criteria for answering this question are:

- Can most of the app be written in shared code?
- Can imaging algorithms be implemented so that they have a fast execution speed, particularly on Android?

The answer to both questions, is yes. I was happy with the speed of execution of the Fourier transform on both platform, especically when considering the size of the source image (2048x2048), and the sheer amount of processing that’s performed when frequency filtering an image. In addition, there are plenty of opportunities to further optimise my implementation, as my implementation focus was clarity rather than optimisation. In particular, a 2D FFT naturally lends itself to parallelisation.

While consumer imaging apps don’t typically contain operations that Fourier transform image data, it’s a mainstay of scientific imaging. Therefore, it’s safe to say that Xamarin.Forms and SkiaSharp also make a good combination for scientific imaging apps.

The sample this code comes from can be found on GitHub.

## No comments:

## Post a Comment