Filtering and restoration

TOP HAT FILTER

%matplotlib inline

import numpy as np

import matplotlib.pyplot as plt

from skimage import data

from skimage import color, morphology

image = color.rgb2gray(data.hubble_deep_field())[:500, :500]

selem = morphology.disk(1)

res = morphology.white_tophat(image, selem)

fig, ax = plt.subplots(ncols=3, figsize=(20, 8))

ax[0].set_title('Original')

ax[0].imshow(image, cmap='gray')

ax[1].set_title('White tophat')

ax[1].imshow(res, cmap='gray')

ax[2].set_title('Complementary')

ax[2].imshow(image - res, cmap='gray')

plt.show()

apply_hysteresis_threshold

%matplotlib inline

import matplotlib.pyplot as plt

from skimage import data, filters

fig, ax = plt.subplots(nrows=2, ncols=2)

image = data.coins()

edges = filters.sobel(image)

low = 0.1

high = 0.35

lowt = (edges > low).astype(int)

hight = (edges > high).astype(int)

hyst = filters.apply_hysteresis_threshold(edges, low, high)

ax[0, 0].imshow(image, cmap='gray')

ax[0, 0].set_title('Original image')

ax[0, 1].imshow(edges, cmap='magma')

ax[0, 1].set_title('Sobel edges')

ax[1, 0].imshow(lowt, cmap='magma')

ax[1, 0].set_title('Low threshold')

ax[1, 1].imshow(hight + hyst, cmap='magma')

ax[1, 1].set_title('Hysteresis threshold')

for a in ax.ravel():

a.axis('off')

plt.tight_layout()

plt.show()

Image Deconvolution

%matplotlib inline

import numpy as np

import matplotlib.pyplot as plt

from skimage import color, data, restoration

astro = color.rgb2gray(data.astronaut())

from scipy.signal import convolve2d as conv2

psf = np.ones((5, 5)) / 25

astro = conv2(astro, psf, 'same')

astro += 0.1 * astro.std() * np.random.standard_normal(astro.shape)

deconvolved, _ = restoration.unsupervised_wiener(astro, psf)

fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(8, 5),

sharex=True, sharey=True)

plt.gray()

ax[0].imshow(astro, vmin=deconvolved.min(), vmax=deconvolved.max())

ax[0].axis('off')

ax[0].set_title('Data')

ax[1].imshow(deconvolved)

ax[1].axis('off')

ax[1].set_title('Self tuned restoration')

fig.tight_layout()

plt.show()

Window + FFT (frekuensi)

%matplotlib inline

import matplotlib.pyplot as plt

import numpy as np

from scipy.fftpack import fft2, fftshift

from skimage import img_as_float

from skimage.color import rgb2gray

from skimage.data import astronaut

from skimage.filters import window

image = img_as_float(rgb2gray(astronaut()))

wimage = image * window('hann', image.shape)

image_f = np.abs(fftshift(fft2(image)))

wimage_f = np.abs(fftshift(fft2(wimage)))

fig, axes = plt.subplots(2, 2, figsize=(8, 8))

ax = axes.ravel()

ax[0].set_title("Original image")

ax[0].imshow(image, cmap='gray')

ax[1].set_title("Windowed image")

ax[1].imshow(wimage, cmap='gray')

ax[2].set_title("Original FFT (frequency)")

ax[2].imshow(np.log(image_f), cmap='magma')

ax[3].set_title("Window + FFT (frequency)")

ax[3].imshow(np.log(wimage_f), cmap='magma')

plt.show()

Mean Filter

%matplotlib inline

import matplotlib.pyplot as plt

from skimage import data

from skimage.morphology import disk

from skimage.filters import rank

image = data.coins()

selem = disk(20)

percentile_result = rank.mean_percentile(image, selem=selem, p0=.1, p1=.9)

bilateral_result = rank.mean_bilateral(image, selem=selem, s0=500, s1=500)

normal_result = rank.mean(image, selem=selem)

fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(10, 10),

sharex=True, sharey=True)

ax = axes.ravel()

titles = ['Original', 'Percentile mean', 'Bilateral mean', 'Local mean']

imgs = [image, percentile_result, bilateral_result, normal_result]

for n in range(0, len(imgs)):

ax[n].imshow(imgs[n], cmap=plt.cm.gray)

ax[n].set_title(titles[n])

ax[n].axis('off')

plt.tight_layout()

plt.show()

unsharp_mask

from skimage import data

from skimage.filters import unsharp_mask

import matplotlib.pyplot as plt

image = data.moon()

result_1 = unsharp_mask(image, radius=1, amount=1)

result_2 = unsharp_mask(image, radius=5, amount=2)

result_3 = unsharp_mask(image, radius=20, amount=1)

fig, axes = plt.subplots(nrows=2, ncols=2,

sharex=True, sharey=True, figsize=(10, 10))

ax = axes.ravel()

ax[0].imshow(image, cmap=plt.cm.gray)

ax[0].set_title('Original image')

ax[1].imshow(result_1, cmap=plt.cm.gray)

ax[1].set_title('Enhanced image, radius=1, amount=1.0')

ax[2].imshow(result_2, cmap=plt.cm.gray)

ax[2].set_title('Enhanced image, radius=5, amount=2.0')

ax[3].imshow(result_3, cmap=plt.cm.gray)

ax[3].set_title('Enhanced image, radius=20, amount=1.0')

for a in ax:

a.axis('off')

fig.tight_layout()

plt.show()

Image Deconvolution

import numpy as np

import matplotlib.pyplot as plt

from scipy.signal import convolve2d as conv2

from skimage import color, data, restoration

astro = color.rgb2gray(data.astronaut())

psf = np.ones((5, 5)) / 25

astro = conv2(astro, psf, 'same')

# Add Noise to Image

astro_noisy = astro.copy()

astro_noisy += (np.random.poisson(lam=25, size=astro.shape) - 10) / 255.

# Restore Image using Richardson-Lucy algorithm

deconvolved_RL = restoration.richardson_lucy(astro_noisy, psf, iterations=30)

fig, ax = plt.subplots(nrows=1, ncols=3, figsize=(8, 5))

plt.gray()

for a in (ax[0], ax[1], ax[2]):

a.axis('off')

ax[0].imshow(astro)

ax[0].set_title('Original Data')

ax[1].imshow(astro_noisy)

ax[1].set_title('Noisy data')

ax[2].imshow(deconvolved_RL, vmin=astro_noisy.min(), vmax=astro_noisy.max())

ax[2].set_title('Restoration using\nRichardson-Lucy')

fig.subplots_adjust(wspace=0.02, hspace=0.2,

top=0.9, bottom=0.05, left=0, right=1)

plt.show()

Inpainting

import numpy as np

import matplotlib.pyplot as plt

from skimage import data

from skimage.restoration import inpaint

image_orig = data.astronaut()[0:200, 0:200]

# Create mask with three defect regions: left, middle, right respectively

mask = np.zeros(image_orig.shape[:-1])

mask[20:60, 0:20] = 1

mask[160:180, 70:155] = 1

mask[30:60, 170:195] = 1

# Defect image over the same region in each color channel

image_defect = image_orig.copy()

for layer in range(image_defect.shape[-1]):

image_defect[np.where(mask)] = 0

image_result = inpaint.inpaint_biharmonic(image_defect, mask,

multichannel=True)

fig, axes = plt.subplots(ncols=2, nrows=2)

ax = axes.ravel()

ax[0].set_title('Original image')

ax[0].imshow(image_orig)

ax[1].set_title('Mask')

ax[1].imshow(mask, cmap=plt.cm.gray)

ax[2].set_title('Defected image')

ax[2].imshow(image_defect)

ax[3].set_title('Inpainted image')

ax[3].imshow(image_result)

for a in ax:

a.axis('off')

fig.tight_layout()

plt.show()

Calibrating Denoisers dengsn J-Invariance

import numpy as np

from matplotlib import pyplot as plt

from skimage.data import chelsea

from skimage.restoration import calibrate_denoiser, denoise_wavelet

from skimage.util import img_as_float, random_noise

from functools import partial

# rescale_sigma=True required to silence deprecation warnings

_denoise_wavelet = partial(denoise_wavelet, rescale_sigma=True)

image = img_as_float(chelsea())

sigma = 0.3

noisy = random_noise(image, var=sigma ** 2)

# Parameters to test when calibrating the denoising algorithm

parameter_ranges = {'sigma': np.arange(0.1, 0.3, 0.02),

'wavelet': ['db1', 'db2'],

'convert2ycbcr': [True, False],

'multichannel': [True]}

# Denoised image using default parameters of `denoise_wavelet`

default_output = denoise_wavelet(noisy, multichannel=True, rescale_sigma=True)

# Calibrate denoiser

calibrated_denoiser = calibrate_denoiser(noisy,

_denoise_wavelet,

denoise_parameters=parameter_ranges)

# Denoised image using calibrated denoiser

calibrated_output = calibrated_denoiser(noisy)

fig, axes = plt.subplots(1, 3, sharex=True, sharey=True, figsize=(15, 5))

for ax, img, title in zip(

axes,

[noisy, default_output, calibrated_output],

['Noisy Image', 'Denoised (Default)', 'Denoised (Calibrated)']

ax.imshow(img)

ax.set_title(title)

ax.set_yticks([])

ax.set_xticks([])

plt.show()

Denoising

import matplotlib.pyplot as plt

from skimage.restoration import (denoise_tv_chambolle, denoise_bilateral,

denoise_wavelet, estimate_sigma)

from skimage import data, img_as_float

from skimage.util import random_noise

original = img_as_float(data.chelsea()[100:250, 50:300])

sigma = 0.155

noisy = random_noise(original, var=sigma**2)

fig, ax = plt.subplots(nrows=2, ncols=4, figsize=(8, 5),

sharex=True, sharey=True)

plt.gray()

# Estimate the average noise standard deviation across color channels.

sigma_est = estimate_sigma(noisy, multichannel=True, average_sigmas=True)

# Due to clipping in random_noise, the estimate will be a bit smaller than the

# specified sigma.

print(f"Estimated Gaussian noise standard deviation = {sigma_est}")

ax[0, 0].imshow(noisy)

ax[0, 0].axis('off')

ax[0, 0].set_title('Noisy')

ax[0, 1].imshow(denoise_tv_chambolle(noisy, weight=0.1, multichannel=True))

ax[0, 1].axis('off')

ax[0, 1].set_title('TV')

ax[0, 2].imshow(denoise_bilateral(noisy, sigma_color=0.05, sigma_spatial=15,

multichannel=True))

ax[0, 2].axis('off')

ax[0, 2].set_title('Bilateral')

ax[0, 3].imshow(denoise_wavelet(noisy, multichannel=True, rescale_sigma=True))

ax[0, 3].axis('off')

ax[0, 3].set_title('Wavelet denoising')

ax[1, 1].imshow(denoise_tv_chambolle(noisy, weight=0.2, multichannel=True))

ax[1, 1].axis('off')

ax[1, 1].set_title('(more) TV')

ax[1, 2].imshow(denoise_bilateral(noisy, sigma_color=0.1, sigma_spatial=15,

multichannel=True))

ax[1, 2].axis('off')

ax[1, 2].set_title('(more) Bilateral')

ax[1, 3].imshow(denoise_wavelet(noisy, multichannel=True, convert2ycbcr=True,

rescale_sigma=True))

ax[1, 3].axis('off')

ax[1, 3].set_title('Wavelet denoising\nin YCbCr colorspace')

ax[1, 0].imshow(original)

ax[1, 0].axis('off')

ax[1, 0].set_title('Original')

fig.tight_layout()

plt.show()

Denoising

import numpy as np

import matplotlib.pyplot as plt

from skimage.morphology import diameter_closing

from skimage import data

from skimage.morphology import closing

from skimage.morphology import square

datasets = {

'retina': {'image': data.microaneurysms(),

'figsize': (15, 9),

'diameter': 10,

'vis_factor': 3,

'title': 'Detection of microaneurysm'},

'page': {'image': data.page(),

'figsize': (15, 7),

'diameter': 23,

'vis_factor': 1,

'title': 'Text detection'}

}

for dataset in datasets.values():

# image with printed letters

image = dataset['image']

figsize = dataset['figsize']

diameter = dataset['diameter']

fig, ax = plt.subplots(2, 3, figsize=figsize)

# Original image

ax[0, 0].imshow(image, cmap='gray', aspect='equal',

vmin=0, vmax=255)

ax[0, 0].set_title('Original', fontsize=16)

ax[0, 0].axis('off')

ax[1, 0].imshow(image, cmap='gray', aspect='equal',

vmin=0, vmax=255)

ax[1, 0].set_title('Original', fontsize=16)

ax[1, 0].axis('off')

# Diameter closing : we remove all dark structures with a maximal

# extension of less than <diameter> (12 or 23). I.e. in closed_attr, all

# local minima have at least a maximal extension of <diameter>.

closed_attr = diameter_closing(image, diameter, connectivity=2)

# We then calculate the difference to the original image.

tophat_attr = closed_attr - image

ax[0, 1].imshow(closed_attr, cmap='gray', aspect='equal',

vmin=0, vmax=255)

ax[0, 1].set_title('Diameter Closing', fontsize=16)

ax[0, 1].axis('off')

ax[0, 2].imshow(dataset['vis_factor'] * tophat_attr, cmap='gray',

aspect='equal', vmin=0, vmax=255)

ax[0, 2].set_title('Tophat (Difference)', fontsize=16)

ax[0, 2].axis('off')

# A morphological closing removes all dark structures that cannot

# contain a structuring element of a certain size.

closed = closing(image, square(diameter))

# Again we calculate the difference to the original image.

tophat = closed - image

ax[1, 1].imshow(closed, cmap='gray', aspect='equal',

vmin=0, vmax=255)

ax[1, 1].set_title('Morphological Closing', fontsize=16)

ax[1, 1].axis('off')

ax[1, 2].imshow(dataset['vis_factor'] * tophat, cmap='gray',

aspect='equal', vmin=0, vmax=255)

ax[1, 2].set_title('Tophat (Difference)', fontsize=16)

ax[1, 2].axis('off')

fig.suptitle(dataset['title'], fontsize=18)

fig.tight_layout(rect=(0, 0, 1, 0.88))

plt.show()

Denoisers Using J-Invariance

import numpy as np

from matplotlib import pyplot as plt

from matplotlib import gridspec

from skimage.data import chelsea, hubble_deep_field

from skimage.metrics import mean_squared_error as mse

from skimage.metrics import peak_signal_noise_ratio as psnr

from skimage.restoration import (calibrate_denoiser,

denoise_wavelet,

denoise_tv_chambolle, denoise_nl_means,

estimate_sigma)

from skimage.util import img_as_float, random_noise

from skimage.color import rgb2gray

from functools import partial

_denoise_wavelet = partial(denoise_wavelet, rescale_sigma=True)

image = img_as_float(chelsea())

sigma = 0.2

noisy = random_noise(image, var=sigma ** 2)

# Parameters to test when calibrating the denoising algorithm

parameter_ranges = {'sigma': np.arange(0.1, 0.3, 0.02),

'wavelet': ['db1', 'db2'],

'convert2ycbcr': [True, False],

'multichannel': [True]}

# Denoised image using default parameters of `denoise_wavelet`

default_output = denoise_wavelet(noisy, multichannel=True, rescale_sigma=True)

# Calibrate denoiser

calibrated_denoiser = calibrate_denoiser(noisy,

_denoise_wavelet,

denoise_parameters=parameter_ranges

)

# Denoised image using calibrated denoiser

calibrated_output = calibrated_denoiser(noisy)

fig, axes = plt.subplots(1, 3, sharex=True, sharey=True, figsize=(15, 5))

for ax, img, title in zip(axes,

[noisy, default_output, calibrated_output],

['Noisy Image', 'Denoised (Default)',

'Denoised (Calibrated)']):

ax.imshow(img)

ax.set_title(title)

ax.set_yticks([])

ax.set_xticks([])

Share on :

Filtering and restoration

5 Responses to "Filtering and restoration "