Destripe satellite image using FFT in python

Question

I am currently in the midst of processing a satellite image, and I've come to a point where I need to remove those faint stripes caused by the sensor, as you can see below:

The image is RGB. I converted it to grayscale and performed FFT to visualise it in the frequency domain. These are the results:

What I've come to realize, mainly from here, is that this cross created represents the stripes I want to remove. I don't know how to remove them though, any help would be appreciated.

Claude.AI provided me with:

import numpy as np
import matplotlib.pyplot as plt
from PIL import Image
from scipy.fftpack import fft2, ifft2, fftshift, ifftshift

# Load the image (replace with your image path)
image_path = 'input.tiff'
image = np.array(Image.open(image_path).convert('L'))

# Apply FFT to transform to frequency domain
f_transform = fft2(image)
f_shift = fftshift(f_transform)

# Create the notch filter
rows, cols = image.shape
crow, ccol = rows // 2, cols // 2  # center of the frequency domain

# Create a mask for the filter - start with all ones (no filtering)
mask = np.ones((rows, cols), dtype=np.float32)

# Define the width of the vertical line to filter out
line_width = 3  # Adjust this based on your image

# Create a vertical band filter that removes the vertical line in the FFT
# (which corresponds to horizontal stripes in the image)
mask[:, ccol-line_width:ccol+line_width+1] = 0
mask[crow-line_width:crow+line_width+1, :] = 0

# Apply the filter
f_shift_filtered = f_shift * mask

# Inverse FFT to get back to the spatial domain
f_ishift = ifftshift(f_shift_filtered)
img_back = np.abs(ifft2(f_ishift))

# Save the destriped image
Image.fromarray(img_back.astype(np.uint8)).save('destriped_satellite_image.tiff')

But it provides me with a result I don't get:

Answer 1

I feel that the other answers address only part of the issue and leave out the elephant in the room. Like TimWescott points out in one of the comments, the faint diagonal stripes in the frequency domain are the ones responsible for the diagonal stripes in the original image, not the cross in the middle. In general, any line pattern in an image is translated to perpendicular line patterns in the frequency domain, and you'll see that those diagonal patterns are perpendicular (or at least close to perpendicular)

Also, I'm not sure where you got the idea that zeroing the cross in the middle would get rid of the noise, I don't think the link you provided advocates for that nuclear choice, but by a more localized intervention in a specific region of the FT.

It's been a long time since I've worked with 2D FFTs and I wish I could give you more guidance, but I feel that it was worth it to give more visibility to some of the comments in the form of an answer.

[EDIT 2025-04-26] After seeing some of the comments, I've noticed that an essential property of the 2D continuous FT (which is great to get an intuition over the results in the discrete FFT case) is not so well-known (and I'm surprised that I remember it at all!): lines in the spatial domain result in perpendicular lines in the frequency domain. I wrote that in my original answer, but I think this didn't get enough attention.

Let's see first that a single line results in one single perpendicular line in the frequency domain (notice my non-compromising claim). A line can be represented with a delta such as $ \delta(Ax + By + C) $. For the sake of simplicity, I'll assume that the line goes through the origin, i.e. $ C = 0 $, in the following derivations. It's not important, $ C \neq 0 $ results only in a phase shift. Let's calculate the 2D FT:

$$ \mathcal{F}\left\{ \delta(Ax + By) \right\} = \int_{-\infty}^\infty \int_{-\infty}^\infty \delta(Ax + By) \exp\left\{ -2 \pi i (ux + vy) \right\} \mathrm{d}x \mathrm{d}y = $$ $$ \int_{-\infty}^\infty \exp\left\{ -2 \pi i (u - v \frac{A}{B}) x \right\} \mathrm{d}x $$

We know from the relation between the 1D delta and its transform (a complex exponential), that the FT of a complex exponential must be a delta (why? Because the FT and the IFT are virtually identical operations). Specifically, $ \delta(u - \frac{A}{B} v) $. And the line represented by this delta is perpendicular to the line in the spatial domain. One way to see this is to take the normal vectors (A,B) and (1,-A/B) and calculate the scalar product: $ 1 \cdot A - \frac{A}{B} \cdot B = 0 $.

Now, some empirical evidence for this fact. Here's the FFT of a straight line:

Naturally, the FFT is not a perfect line for two reasons: the effect of discretization and the effect of the rectangular window.

Answer 2

I think you're basically seeing what you'd expect: since the DFT is circular on the image borders, removing the cross in the center is like multiplying the spatial frequency domain with a very box window that's the image dimensions - 1px on each side (just shifted so that the center of the box is at the corners of the image).

Now something you've probably heard before:

Multiplication in one domain is equivalent to convolution of the Fourier-transformed factors in the time domain

(or in your case, spatial domain)

That's why Why is it a bad idea to filter by zeroing out FFT bins? is such a popula question :)

You're convolving your image with a sinc shaped kernel, and having removed the center bin, you do that on an image where the average has been removed, which probably perceptually amplifies noise in the formerly dark regions of the image.

Instead, you probably want to design a more sensible filter that you convolve with :)

Answer 3

First, you should for spectral analysis purposes generally apply a window before taking the FFT. This avoids artifacts from the image edges from obscuring the actually interesting properties (in this case, the square cross which the AI tried to remove).

The feature we should actually get rid of are the tilted stripes in the frequency domain. But you should not brutally erase them from the spectrum, in particular not delete a whole line through the center – that removes also low-frequency components that are important part of the clean signal.
Instead, you should only reduce the amplitudes enough to make the spectrum isotropic.

One way to accomplish that is to transform the spectrum into polar coordinates^† and compute the median of each row (i.e. of each circle of direction-independent frequency). Broadcast this to the original shape and transform it back to cartesian to get your target spectrum.

Then you can compute the quotient of the target spectrum and the actually observed one, to figure out how the image needs to be tweaked. This you can then apply as a multiplier to the (complex) spectrum of the image without window function^‡. iFFT back and the stripes should be mostly removed, without too much change to the signal.

---

^†_{A problem with the polar transform is that it's only defined in a disc-shaped region. It is probably ok to just zero the spectrum outside of the disc, but you can also pad the spectrum before going to polar coordinates (taking the median makes the processing somewhat stable to this).}

^‡_{Alternatively, do the entire processing with windowed chunks. This is the preferred strategy if the image is too big for a single FFT. Hann- or Tukey windows have the property that if you add together two overlapping neighbouring windowed chunks, the amplitudes add to unity. This avoids getting visible mismatches along the edges of the processed tiles.}

Answer 4

As @Asher mentioned, do not remove the cross in the middle, since it contains the image information you want to preserve. Instead, remove the lines perpendicular to your artifacts.

I drew over the diagonal stripes in black with GIMP and got this mask:

I then filtered the image by multiplying its FFT with the mask.

import numpy as np
from PIL import Image

image = np.array(Image.open("ONkM3J18.png").convert("RGB")) / 255.0
mask = np.array(Image.open("mask.png").convert("L")) / 255.0

z = [np.fft.fftshift(np.fft.fft2(c)) for c in image.transpose(2, 0, 1)]

if 0:
    # Save FFT for painting over mask
    z = np.log(np.abs(np.dstack(z)) + 1)
    z = (z - z.min()) / (z.max() - z.min())
    Image.fromarray((z * 255).astype(np.uint8)).save("fft.png")
else:
    # Save FFT-filtered image
    image = np.dstack([np.fft.ifft2(np.fft.ifftshift(c * mask)) for c in z])
    image = np.clip(np.abs(image) * 255, 0, 255).astype(np.uint8)
    Image.fromarray(image).save("filtered.png")

The filtered image has much fewer stripe artifacts. For comparison, open both the original and the filtered image in a new browser tab and then flip back and forth between the two tabs to see the difference.

Answer 5

You erased the DC component, which is the center of the cross.

That causes the picture to become "black", or rather, its mean will be 0, and about half the image content will be below 0, i.e. darker than black.

Do not erase it. Then it'll look right. Or better anyway.