问在噪声图像中计数物体

EN

from scipy import ndimage
from skimage import filters, feature, measure, color
import numpy as np
import tifffile
import matplotlib.pyplot as plt
from scipy.ndimage import median_filter, morphology
from scipy.fft import fftfreq, fft2, fftshift, ifft2, ifftshift
from skimage.segmentation import watershed

plt.rcParams["image.cmap"] = "gray"
plt.rcParams["text.color"] = "white"
plt.rcParams["xtick.color"] = "gray"
plt.rcParams["ytick.color"] = "gray"

def fourier_spectrum_to_log_spectrum(im): # reduce max and min values with logarithm
    return np.log10(1 + np.abs(im) ** 2)

def make_notch(first_shape, locations, sigma): # implemented from exercise.ipynb from moodle and slightly modified
    u = fftfreq(first_shape, 1 / first_shape)
    u, v = np.meshgrid(u, u)
    res = np.zeros((first_shape, first_shape))
    for (v_0, u_0) in locations:
        diff = [(v-v_0), (u-u_0)]
        sums = [(v+v_0), (u+u_0)]
        res += np.exp(-(diff[0]  ** 2 + diff[1] ** 2) / (2 * sigma ** 2))
        res += np.exp(-(sums[0]  ** 2 + sums[1] ** 2) / (2 * sigma ** 2))
        res += np.exp(-(sums[0]  ** 2 + diff[1] ** 2) / (2 * sigma ** 2))
        res += np.exp(-(diff[0]  ** 2 + sums[1] ** 2) / (2 * sigma ** 2))
    return res

orig_img = tifffile.imread('input.tif')
orig_img_filtered = median_filter(orig_img, 6)

orig_fft = fftshift(fft2(orig_img_filtered))
orig_fft_log = fourier_spectrum_to_log_spectrum(orig_fft)

brightness_threshold = 10
sigma = 20
highpass_filter = 1 - make_notch(orig_fft.shape[0], [np.array(orig_fft.shape) / 2], sigma)
bright_pixels = morphology.binary_dilation((fourier_spectrum_to_log_spectrum(orig_fft) * highpass_filter) > brightness_threshold, iterations=3)

orig_fft_no_grid = (1 - bright_pixels) * orig_fft # use the inverted Bright Pixels on the original Image-FT
orig_img_no_grid = ifft2(ifftshift(orig_fft_no_grid)).real # revert the FT

thresh = filters.threshold_otsu(orig_img_no_grid)
print("Threshold: %f" % thresh)
holes = morphology.binary_fill_holes(orig_img_no_grid < thresh) # mark corns in b/w-Picture
img_opened = morphology.binary_opening(holes, iterations=2) # remove Noise-Corns
img_eroded = morphology.binary_erosion(img_opened, iterations=4) # needed cause some peppercorns touch the Image-Borders

dist_trans = ndimage.distance_transform_edt(img_eroded)
local_max = feature.peak_local_max(dist_trans, min_distance=7)
local_max_mask = np.zeros(dist_trans.shape, dtype=bool)
local_max_mask[tuple(local_max.T)] = True
labels = watershed(-dist_trans, measure.label(local_max_mask), mask=img_eroded) # separate merged corns
fig, ax = plt.subplots(ncols=4, nrows=3, figsize=(30,20))
ax[0][0].set_title("Original")
ax[0][0].imshow(orig_img)
ax[0][1].set_title("Original Log-Fourier")
ax[0][1].imshow(orig_fft_log)
ax[0][2].set_title("Noise-Filtered")
ax[0][2].imshow(orig_img_filtered)
ax[0][3].set_title("Filtered Log-Fourier")
ax[0][3].imshow(fourier_spectrum_to_log_spectrum(fftshift(fft2(orig_img_filtered))))
ax[1][0].set_title("Highpass-Filter inverted")
ax[1][0].imshow(highpass_filter)
ax[1][1].set_title("Bright Pixels w/o Center")
ax[1][1].imshow(bright_pixels)
ax[1][2].set_title("No-Grid Fourier")
ax[1][2].imshow(fourier_spectrum_to_log_spectrum(orig_fft_no_grid))
ax[1][3].set_title("No-Grid")
ax[1][3].imshow(orig_img_no_grid)
ax[2][0].set_title("Filled Holes")
ax[2][0].imshow(holes)
ax[2][1].set_title("Opened Image after Holes")
ax[2][1].imshow(img_opened)
ax[2][2].set_title("Eroded after Opening")
ax[2][2].imshow(img_eroded)
ax[2][3].set_title("Watershed Transform")
ax[2][3].imshow(color.label2rgb(labels, bg_label=0))
print("Number of Peppercorns: %d" % labels.max())

0. invert (objects white)
1. gaussian blur
2a. threshold for absolute response, just high enough to not pick up any background
2b. NMS via `dilate(blurred) == blurred`, dilation radius = NMS radius, roundish kernel is better than square
3. logical `and`
4. cosmetics (dilate the peaks to make them fat)