11

I am trying to do speckle noise removal in satellite SAR image.I am not getting any package which does speckle noise removal in SAR image. I have tried pyradar but it works with python 2.7 and I am working on Anaconda with python 3.5 on windows. Also Rsgislib is available but it is on Linux. Joseph meiring has also given a Lee filter code on github but it fails to work. : https://github.com/reptillicus/LeeFilter

Kindly, can anyone share the python script for Speckle Filter or how to proceed for speckle filter design in python.

Shubham_geo
  • 368
  • 1
  • 4
  • 16

2 Answers2

20

This is a fun little problem. Rather than try to find a library for it, why not write it from the definition?

from scipy.ndimage.filters import uniform_filter
from scipy.ndimage.measurements import variance

def lee_filter(img, size):
    img_mean = uniform_filter(img, (size, size))
    img_sqr_mean = uniform_filter(img**2, (size, size))
    img_variance = img_sqr_mean - img_mean**2

    overall_variance = variance(img)

    img_weights = img_variance / (img_variance + overall_variance)
    img_output = img_mean + img_weights * (img - img_mean)
    return img_output

If you don't want the window to be a square of size x size, just replace uniform_filter with something else (convolution with a disk, gaussian filter, etc). Any type of (weighted) averaging filter will do, as long as it is the same for calculating both img_mean and img_square_mean.

The Lee filter seems rather old-fashioned as a filter. It won't behave well at edges because for any window that has an edge in it, the variance is going to be much higher than the overall image variance, and therefore the weights (of the unfiltered image relative to the filtered image) are going to be close to 1.

An example:

from pylab import *
import numpy as np
img = np.random.normal(0.5, 0.1, (100,100))
img[:,:50] += 0.25
imshow(img, vmin=0, vmax=1, cmap='gray')
imshow(lee_filter(img, 20), vmin=0, vmax=1, cmap='gray')

Noisy image with edge Lee filtered

As you can see the noise reduction is very good in general, but much weaker along the edge.

I'm not familiar with SAR so I don't know if Lee filter has some features that make it particularly good for speckle in SAR, but you may want to look into modern edge-aware denoisers, like guided filter or bilateral filter.

Alex I
  • 19,689
  • 9
  • 86
  • 158
  • Thanks a lot!!!! Alex, actually I am new to this domain and was not able to find a good literature.Thanks a lot!!!!!!!! – Shubham_geo Sep 30 '16 at 11:25
  • 1
    @Shubham_geo You're welcome. I added some more notes and an example. – Alex I Nov 28 '16 at 19:26
  • pyradar package seem to be an issue for the images i have, could you point us to right place where we can find definitions for all filters here just like the one you have above ? https://pyradar-tools.readthedocs.io/en/latest/examples.html#example-of-filtros – bicepjai Nov 08 '17 at 08:24
  • 1
    This is a great, simple implementation. I think there is one small error in your code though. According to [this article](https://angeljohnsy.blogspot.com/2014/08/lee-filter.html), the weight parameter is calculated as the kernel variance divided by the sum of the kernel variance plus the overall image variance. When you calculate the weights, you square the variances instead. The filter works essentially the same but smoothing is more severe in homogeneous regions. – Addison Mar 12 '18 at 06:52
  • @Addison: Yes, you are correct, variance is "sigma squared", not sigma. Well spotted! It's fixed now. – Alex I Apr 11 '18 at 21:27
  • @bicepjai: If you'd like to use multiple filters from pyradar, why not fix the problem in the package? I'm sure they would appreciate some patches/contributions. If you'd like, you can post a new question here with the specific problem. – Alex I Apr 11 '18 at 21:30
  • @Benjamin - Perhaps, but my implementation should be identical to pyradar: https://pyradar-tools.readthedocs.io/en/latest/_modules/pyradar/filters/lee.html I have not read the original Lee paper that closely: https://www.researchgate.net/profile/Patrick_Wambacq/publication/239659062_Speckle_filtering_of_synthetic_aperture_radar_images_A_Review/links/56a5d56508ae232fb20976a5.pdf please let me know if my code (and pyradar) are wrong according to that. – Alex I Feb 22 '19 at 02:09
  • @AlexI - how would you change the code to account also for `np.nan` values? For example when you clip a raster according to a polygon then you will have values outside of the polygon which will be `np.nan`, you can't remove them, so you need somehow to ignore them. – user88484 Feb 15 '21 at 13:46
6

In general it is very difficult to see the effect of a noise filter by eye in a 2D plot. Let me demonstrate this by example. Suppose we have this noisy picture:

noise example

Let me now convert this image to a 3d mesh plot. Then it will look like this. The noise becomes very clear but also the differences in depth between the left and right side of the picture.

mesh plot

The library findpeaks contains many filters which are utilized from various (old python 2) libraries and rewritten to python 3. Applying the filters is very easy as shown below. Note that this example seems not very representative for a SAR image as there is no speckle noise. A mean or median filter seems to perform very well in this example. In speckle noise images where local heights are important, such mean/median filters can remove the peaks and thus destroy the signal of interest.

Install by:

pip install findpeaks

Run by:

from findpeaks import findpeaks

# Read image
img = cv2.imread('noise.png')

filters = [None, 'lee','lee_enhanced','kuan', 'fastnl','bilateral','frost','median','mean']

for getfilter in filters:
    fp = findpeaks(method='topology', scale=False, denoise=getfilter, togray=True, imsize=False, window=15)
    fp.fit(img)
    fp.plot_mesh(wireframe=False, title=str(getfilter), view=(30,30))

Comparison methods

If you directly want to use the denoising filters, it can be done as following:

import findpeaks
import matplotlib.pyplot as plt

# Read image
img = cv2.imread('noise.png')

# filters parameters
# window size
winsize = 15
# damping factor for frost
k_value1 = 2.0
# damping factor for lee enhanced
k_value2 = 1.0
# coefficient of variation of noise
cu_value = 0.25
# coefficient of variation for lee enhanced of noise
cu_lee_enhanced = 0.523
# max coefficient of variation for lee enhanced
cmax_value = 1.73

# Some pre-processing
# Make grey image
img = findpeaks.stats.togray(img)
# Scale between [0-255]
img = findpeaks.stats.scale(img)

# Denoising
# fastnl
img_fastnl = findpeaks.stats.denoise(img, method='fastnl', window=winsize)
# bilateral
img_bilateral = findpeaks.stats.denoise(img, method='bilateral', window=winsize)
# frost filter
image_frost = findpeaks.frost_filter(img, damping_factor=k_value1, win_size=winsize)
# kuan filter
image_kuan = findpeaks.kuan_filter(img, win_size=winsize, cu=cu_value)
# lee filter
image_lee = findpeaks.lee_filter(img, win_size=winsize, cu=cu_value)
# lee enhanced filter
image_lee_enhanced = findpeaks.lee_enhanced_filter(img, win_size=winsize, k=k_value2, cu=cu_lee_enhanced, cmax=cmax_value)
# mean filter
image_mean = findpeaks.mean_filter(img, win_size=winsize)
# median filter
image_median = findpeaks.median_filter(img, win_size=winsize)


plt.figure(); plt.imshow(img_fastnl, cmap='gray'); plt.title('Fastnl')
plt.figure(); plt.imshow(img_bilateral, cmap='gray'); plt.title('Bilateral')
plt.figure(); plt.imshow(image_frost, cmap='gray'); plt.title('Frost')
plt.figure(); plt.imshow(image_kuan, cmap='gray'); plt.title('Kuan')
plt.figure(); plt.imshow(image_lee, cmap='gray'); plt.title('Lee')
plt.figure(); plt.imshow(image_lee_enhanced, cmap='gray'); plt.title('Lee Enhanced')
plt.figure(); plt.imshow(image_mean, cmap='gray'); plt.title('Mean')
plt.figure(); plt.imshow(image_median, cmap='gray'); plt.title('Median')

If you want to play around with the library, more examples can be found over here.

erdogant
  • 1,544
  • 14
  • 23