I'm trying to modify the sorrounding values around clusters in an image so that the neighbours pixels decrease to zero following a exponential decay. I need to find a way to control the decaying rate with a parameter. I also want to keep the original non-zero values in the image.

I have created an example were I used convolution (with average kernel) to create the haloe of decaying values. Unfortunately, it does not create what I really want, since the values decay too quickly, and the non zero values are not kept. I tried as well using a gaussian kernel and varying the sigma value so I can control the spread, the problem is that it reduces the value too much near the clusters. Maybe using a custom filter could help?

import matplotlib.pyplot as plt
import numpy as np
from scipy.ndimage import convolve

matrix = np.zeros((100, 100))

# Create clusters of values
matrix[45:55, 40:50] = 0.5
for i, value in enumerate(np.linspace(0.5, 0.2, 20)):
    matrix[70:90, 20 + i] = value  # Decreasing along the columns

# Create an average kernel (15x15 window)
size = 3
kernel = np.ones((size, size)) / (
    size * size
)  # Normalize the kernel to get the average

# Apply the average filter
filtered_matrix = convolve(matrix, kernel, mode="nearest")

# Combine the original values with the filtered (halo) values
result = np.where(matrix > 0, matrix, filtered_matrix)

# Plot the original, filtered, and combined matrices
fig, ax = plt.subplots(1, 3, figsize=(18, 6))

# Original matrix
im1 = ax[0].imshow(matrix, cmap="viridis")
ax[0].set_title("Original Matrix")
plt.colorbar(im1, ax=ax[0], label="Increment [m]", fraction=0.046, pad=0.04)

# Filtered matrix (after applying the average kernel)
im2 = ax[1].imshow(filtered_matrix, cmap="viridis")
ax[1].set_title("Filtered Matrix with Average Kernel")
plt.colorbar(im2, ax=ax[1], label="Increment [m]", fraction=0.046, pad=0.04)

# Final result with original values preserved and haloes included
im3 = ax[2].imshow(result, cmap="viridis")
ax[2].set_title("Result with Haloes Included")
plt.colorbar(im3, ax=ax[2], label="Increment [m]", fraction=0.046, pad=0.04)

plt.tight_layout()
plt.show()

解决方案

Note that Convolution with Gaussian kernel represents an isotropic diffusion process that will decrease the intensity of objects isotropically, so pasting the original image on top will destroy the halo.

Something like the following, using mathematical morphology should work. Iteratively expand the contour of the objects (e.g., with dilation). Everytime you expand the contour of objects , ensure that the added pixel intensities are decayed exponentially (with a decay factor with which you can control the spread of the halo).

The decay_factor, num_iter (the number of iterations) and the structuring element (footprint) type (e.g., disk or square) and size (se_sz) are the parameters you can vary to get the desired effect.

from skimage.morphology import dilation, disk, square

decay_factor = 0.05
se_sz = 3
num_iter = 7

out_matrix = matrix.copy()
prev_filtered_matrix = matrix.copy()

for i in range(num_iter):
    filtered_matrix = dilation(prev_filtered_matrix, square(se_sz))
    decayed_expansion = (filtered_matrix - prev_filtered_matrix)*np.exp(-decay_factor*i)
    out_matrix = out_matrix + decayed_expansion # add exp-decayed expansion 
    prev_filtered_matrix = filtered_matrix

# Plot the original, filtered, and combined matrices
fig, ax = plt.subplots(1, 2, figsize=(18, 6))

plt.gray()
# Original matrix
im1 = ax[0].imshow(matrix)
ax[0].set_title("Original Matrix")
plt.colorbar(im1, ax=ax[0], label="Increment [m]", fraction=0.046, pad=0.04)

# Final result with original values preserved and haloes included
im3 = ax[1].imshow(out_matrix)
ax[1].set_title("Result with Haloes Included")
plt.colorbar(im3, ax=ax[1], label="Increment [m]", fraction=0.046, pad=0.04)

plt.tight_layout()
plt.show() 

The output obtained with square SE is the one shown below:

Note that here we are not changing the original image (so no change in original image intensity), only adding a fringe layer (with decayed intensity) at every iteration, eliminating the need to paste the original image to the output in the end.

You can optionally keep the images normalized (e.g., scaling the maximum intensity value to 1) at each iteration.

The next animation shows halo expansion: