Sample Code

Halftone Descreening with 2D Fast Fourier Transform

Reduce or remove periodic artifacts from images.

Download

Overview

Accelerate’s vDSP module provides functions to perform 2D fast Fourier transforms (FFTs) on matrices of data, such as images. You can exploit the amplitude peaks in the frequency domain of periodic patterns, such as halftone screens, to reduce or remove such artifacts from images. The example below shows an image with halftone artifacts (left) and the same image with the halftone artifacts reduced (right):

Photographs showing before and after images.

This sample walks you through the following steps to implement halftone descreening, using a halftone screen sample that you’ve created programmatically or taken from source material:

  1. Convert the image data to a split complex vector.

  2. Prepare the FFT setup.

  3. Perform forward 2D FFT on image data.

  4. Zero the peaks in the halftone sample magnitude.

  5. Descreen the source image.

  6. Perform inverse 2D FFT on the source image frequency domain data.

  7. Generate an image from a split complex vector.

The sample works on a single-channel, monochrome image.

Convert Image Data to a Split Complex Vector

Create a split complex vector—suitable for use with vDSP’s 2D FFT—by copying odd pixels to the real parts and even pixels to the imaginary parts of an array of complex numbers. Use the following code to convert the 8-bit unsigned integer image data to floating-point data with which the 2D FFT routine works:

let pixelCount = Int(image.size.width * image.size.height)

let pixelData = cgImage.dataProvider?.data

let pixelsArray = Array(UnsafeBufferPointer(start: CFDataGetBytePtr(pixelData),
                                            count: pixelCount))

let floatPixels = vDSP.integerToFloatingPoint(pixelsArray,
                                              floatingPointType: Float.self)

With the floating-point values populated, create an array of complex numbers, interleavedPixels, that you pass to convert(interleavedComplexVector:toSplitComplexVector:), which copies the values into a split complex vector:

let interleavedPixels = stride(from: 1, to: floatPixels.count, by: 2).map {
    return DSPComplex(real: floatPixels[$0.advanced(by: -1)],
                      imag: floatPixels[$0])
}

vDSP.convert(interleavedComplexVector: interleavedPixels,
             toSplitComplexVector: &splitComplexOut)

Prepare the 2D FFT Setup

Create a setup object that contains all the information required to perform the forward and inverse 2D FFT operations. Creating this setup object can be expensive, so do it only once—for example, when your app is starting—and reuse it.

The following code creates a setup object suitable for performing forward and inverse 2D FFTs on a 1024 x 1024 pixel image:

static let fftSetUp = vDSP.FFT2D(width: 1024,
                                 height: 1024,
                                 ofType: DSPSplitComplex.self)

Prepare Arrays for Transformed Image Data

Create arrays to receive the forward transformed source image and halftone sample:

var sourceImage_floatPixelsReal_frequency = [Float](repeating: 0,
                                                  count: n)
var sourceImage_floatPixelsImag_frequency = [Float](repeating: 0,
                                                    count: n)

var halftoneSample_floatPixelsReal_frequency = [Float](repeating: 0,
                                                       count: n)
var halftoneSample_floatPixelsImag_frequency = [Float](repeating: 0,
                                                       count: n)
var halftoneSampleAmplitude = [Float](repeating: 0,
                                        count: n)

Perform Forward 2D FFTs on Image Data

Use the transform(input:output:direction:) function to perform a forward 2D FFT on the image data, creating the frequency domain representation of the image. Pass transform(input:output:direction:) a split complex structure as the destination, with the same length as the source structure.

The following example shows the code required to perform the FFT on the source image and halftone sample data. After the forward FFT of the halftone sample is is complete, use squareMagnitudes(_:result:) function to compute the magnitudes of the complex values representing the halftone sample:

sourceImage_floatPixelsReal_frequency.withUnsafeMutableBufferPointer { sourceRealPtr in
    sourceImage_floatPixelsImag_frequency.withUnsafeMutableBufferPointer { sourceImagPtr in
        halftoneSample_floatPixelsReal_frequency.withUnsafeMutableBufferPointer { halftoneRealPtr in
            halftoneSample_floatPixelsImag_frequency.withUnsafeMutableBufferPointer { halftoneImagPtr in
                
                var sourceImage_floatPixels_frequency = DSPSplitComplex(
                    realp: sourceRealPtr.baseAddress!,
                    imagp: sourceImagPtr.baseAddress!)
                
                fftSetUp?.transform(input: sourceImageSplitComplex,
                                    output: &sourceImage_floatPixels_frequency,
                                    direction: .forward)
                                    
                var halftoneSample_floatPixels_frequency = DSPSplitComplex(
                    realp: halftoneRealPtr.baseAddress!,
                    imagp: halftoneImagPtr.baseAddress!)
                
                fftSetUp?.transform(input: halftoneSampleSplitComplex,
                                    output: &halftoneSample_floatPixels_frequency,
                                    direction: .forward)
                
                vDSP.squareMagnitudes(halftoneSample_floatPixels_frequency,
                                      result: &halftoneSampleAmplitude)
            }
        }
    }
}

Zero the Peaks in the Halftone Sample Magnitude

You can reduce the halftone screen artifacts by manipulating the magnitude of the frequency domain data for the halftone sample. Zero all the samples above a specified threshold in the magnitudes and clamp the data to 0…1. Then, multiply the frequency domain data of the source image by the manipulated magnitudes.

Use the threshold(_:to:with:) function to set all magnitude values that are over the threshold to -1, and all magnitude values that are less than or equal to the threshold to 1:

let outputConstant: Float = -1

vDSP.threshold(halftoneSampleAmplitude,
               to: threshold,
               with: .signedConstant(outputConstant),
               result: &halftoneSampleAmplitude)

You can now clip the magnitude data between 0 and 1. After clip(_:to:result:) returns, all the high-magnitude values in halftoneSampleAmplitude are set to 0, and all the low-magnitude values are set to 1:

vDSP.clip(halftoneSampleAmplitude,
          to: 0 ... 1,
          result: &halftoneSampleAmplitude)

Descreen the Source Image

Multiply the frequency domain data of the source image by the values in halftoneSampleAmplitude, to remove or reduce the halftone screen:

vDSP.multiply(sourceImage_floatPixels_frequency,
              by: halftoneSampleAmplitude,
              result: &sourceImage_floatPixels_frequency)

Perform Inverse 2D FFT on Frequency Domain Data

You can now perform an inverse 2D FFT to generate a spatial domain version of the image. Use the same fftSetup pointer as you used for the forward 2D FFT, but specify the inverse direction:

fftSetUp?.transform(input: sourceImage_floatPixels_frequency,
                    output: &floatPixels_spatial,
                    direction: .inverse)

Generate an Image from a Split Complex Vector

Your last step is to create a displayable image from the spatial domain representation of the treated source image. The final image is generated from 8-bit, unsigned integers. Because single-precision values can exceed the range of an 8-bit, unsigned integer, clamp the values before converting them:

var low: Float = 0
var high: Float = 255

vDSP_vclip(pixelSource.realp,
           stride,
           &low,
           &high,
           pixelSource.realp,
           stride,
           n)

vDSP_vclip(pixelSource.imagp,
           stride,
           &low,
           &high,
           pixelSource.imagp,
           stride, n)

Use convertElements(of:to:) to convert the floating-point values back to 8-bit unsigned integers that are suitable for generating an image:

var uIntPixels = [UInt8](repeating: 0,
                         count: pixelCount)

let floatPixels = [Float](fromSplitComplex: pixelSource,
                          scale: 1,
                          count: pixelCount)

vDSP.convertElements(of: floatPixels,
                     to: &uIntPixels,
                     rounding: .towardZero)

Now that uIntPixels is populated with the pixel values, generate a CGImage instance from those values:


let result: UIImage? = uIntPixels.withUnsafeMutableBufferPointer { uIntPixelsPtr in
    
    let buffer = vImage_Buffer(data: uIntPixelsPtr.baseAddress!,
                               height: vImagePixelCount(height),
                               width: vImagePixelCount(width),
                               rowBytes: width)
    
    if
        let format = vImage_CGImageFormat(bitsPerComponent: 8,
                                          bitsPerPixel: 8,
                                          colorSpace: CGColorSpaceCreateDeviceGray(),
                                          bitmapInfo: bitmapInfo),
        let cgImage = try? buffer.createCGImage(format: format) {
        
        return UIImage(cgImage: cgImage)
    } else {
        print("Unable to create `CGImage`.")
        return nil
    }
}

See Also

Fourier and Cosine Transforms

Finding the Component Frequencies in a Composite Sine Wave

Use 1D fast Fourier transform to compute the frequency components of a signal.

Using Windowing with Discrete Fourier Transforms

Multiply signal data by window sequence values to reduce spectral leakage.

Signal Extraction from Noise

Use Accelerate’s discrete cosine transform to remove noise from a signal.

Fast Fourier Transforms

Transform vectors and matrices of temporal and spatial domain complex values to the frequency domain and vice versa.

Discrete Fourier Transforms

Transform vectors of temporal and spatial domain complex values to the frequency domain and vice versa.

Discrete Cosine Transforms

Transform vectors of temporal and spatial domain real values to the frequency domain and vice versa.