Convolution with segmentation.make_2dgaussian_kernel for image segmentation

[astqx]

Hi,

I have been working on a pipeline for photometry and recently came across this issue (reproducible code has been attached at the end):

convolution enhancement = convolved data - data
convolved data residual = convolved data [i] - convolved data [i-1]
colorbar: +- 5 * lowest amplitude
green pixels: above threshold

For the Gaussian kernel returned by segmentation.make_2dgaussian_kernel, the convolved data loses intensity in the central region as a result of convolution with a volume normalized kernel, though further surrounding pixels are enhanced. This could still be helpful in increasing the number of connected pixels above a low enough threshold, but may not always provide the required enhancement for image segmentation. Like in the first two cases the faintest object loses 8-connectedness as a result of convolution.

I believe it would be good to have an option to create (return) unnormalized kernels, with something like astropy.convolution.discretize_model to allow more control in cases where appropriate. I understand that using unnormalized kernels could also result in modifying the threshold for the convolved data accordingly, if significant background residual remains after background subtraction.

Aside, an article (Convolving with Unnormalized Kernels) in astropy.convolution documentation discusses:

There are some tasks, such as source finding, where you want to apply a filter with a kernel that is not normalized.
For data that are well-behaved (contain no missing or infinite values), this can be done in one step:

convolve(image, kernel)

The text makes sense but their aforementioned "one step" contradicts default implementation:

astropy.convolution.convolve(array, kernel, ..., normalize_kernel=True, ...)

Interestingly, the reproducible example provided in that article calls:

astropy_conv = convolve(img, kernel, normalize_kernel=False, ...)

to achieve reasonable results:

Reproducible code

import numpy as np
import matplotlib.cm as cm
import matplotlib.colors as mcolors
from astropy.modeling import models
from astropy.visualization import MinMaxInterval
from photutils.segmentation import make_2dgaussian_kernel
from astropy.convolution import convolve
from astropy.convolution.kernels import CustomKernel
import matplotlib.pyplot as plt
from astropy.stats import gaussian_fwhm_to_sigma
from astropy.convolution import discretize_model
from copy import deepcopy

def to_shape(a, shape):
    y_, x_ = shape
    y, x = a.shape
    y_pad = (y_-y)
    x_pad = (x_-x)
    return np.pad(a,((y_pad//2 + y_pad%2, y_pad//2), (x_pad//2 + x_pad%2, x_pad//2)), mode = 'constant')    

amp_high = 1
amp_low = amp_high/10
# model = models.Moffat2D(x_0=0,y_0=0,amplitude=amp,gamma=2,alpha=2) 
all_models = [
    models.Gaussian2D(amplitude=amp_high, x_mean=3, y_mean=3, x_stddev=1, y_stddev=1, theta=0),
    models.Gaussian2D(amplitude=amp_high, x_mean=-0, y_mean=-0, x_stddev=1, y_stddev=1, theta=0),
    models.Gaussian2D(amplitude=amp_low, x_mean=4, y_mean=-4, x_stddev=1, y_stddev=1, theta=0),
    models.Gaussian2D(amplitude=amp_high, x_mean=-5, y_mean=-5, x_stddev=1, y_stddev=2, theta=45),
    models.Moffat2D(amplitude=amp_high, x_0=-5, y_0=5, alpha=2, gamma=1),
    models.Moffat2D(amplitude=amp_high, x_0=12, y_0=12, alpha=2, gamma=2),
]

size=25
grid_size=np.round(size/2)
x,y = np.mgrid[-grid_size:grid_size+1,-grid_size:grid_size+1]

data = all_models[0](x,y)
for model in all_models[1:]:
    data_add = model(x,y)
    data += data_add

thresh_factor = 1
min_factor = 0

add_bkg = True
if add_bkg:
    # bkg = amp_low/10
    bkg=0
    data+=np.ones(data.shape)*bkg
    low_model = (7,17)
    min_val = data[low_model[1]][low_model[0]]
    print('min',min_val)
    bkg_rms = (min_val + bkg)*(1+min_factor/100)

cmap_max_factor = 5
cmap_max = amp_low*cmap_max_factor

SEX_ALL_GROUND_CONV = CustomKernel(np.array([
    [1,2,1],
    [2,4,2],
    [1,2,1]
]))

SEX_ALL_GROUND_CONV_NORM = CustomKernel(np.array([
    [1,2,1],
    [2,4,2],
    [1,2,1]
]))
SEX_ALL_GROUND_CONV_NORM.normalize(mode='integral')

kernel_fwhm = 2
# Gaussian2DKernel(**kernel_param) for kernel_param in kernel_params]
gauss_model = models.Gaussian2D(
    amplitude=1,
    x_mean=0,
    y_mean=0,
    x_stddev=kernel_fwhm * gaussian_fwhm_to_sigma,
    y_stddev=kernel_fwhm * gaussian_fwhm_to_sigma,
    theta=0,
) 

gauss_model_gal = models.Gaussian2D(
    amplitude=1,
    x_mean=0,
    y_mean=0,
    x_stddev=2,
    y_stddev=1,
    theta=45,
) 

moffat_model = models.Moffat2D(
    amplitude=1,
    x_0=0,
    y_0=0,
    alpha=1,
    gamma=1
)

kernels = {
    'SExtractor\n\"all-ground\"\n(normalized)': [SEX_ALL_GROUND_CONV_NORM,False],
    'make_2d-\ngaussian_kernel\nFWHM=2': [make_2dgaussian_kernel(kernel_fwhm,3),False],
    'SE \"all-ground\"\nGaussian amp=4\nFWHM=2': [SEX_ALL_GROUND_CONV,False],
    # 'Gaussian2D\nFWHM=2 (3,3)': [make_2dgaussian_kernel_unnorm(kernel_fwhm,3),False],
    # 'Gaussian\nFWHM=2 (5,5)': [make_2dgaussian_kernel(kernel_fwhm,5),False]
    'discretize_model\nGaussian amp=1\nFWHM=2 (3)': [CustomKernel(discretize_model(gauss_model,x_range=(-1,2),y_range=(-1,2))),False],
    'discretize_model\nMoffat amp=1\nFWHM=2 'r'$\alpha$=1': [CustomKernel(discretize_model(moffat_model,x_range=(-1,2),y_range=(-1,2))),False],
    'discretize_model\nGaussian amp=1\nFWHM=2 (5)': [CustomKernel(discretize_model(gauss_model,x_range=(-2,3),y_range=(-2,3))),False],
    'discretize_model\nGaussian amp=1\n'r'$\sigma_y/\sigma_x=2,\theta=45^\circ$': [CustomKernel(discretize_model(gauss_model_gal,x_range=(-2,3),y_range=(-2,3))),False]
    
}

conv_data = np.array([convolve(data,kernel[0],normalize_kernel=kernel[1]) for kernel in kernels.values()])

kernels_data = [kernel[0].array for kernel in kernels.values()]
kernels_data_padded = [to_shape(kernel_data,data.shape) for kernel_data in kernels_data]

conv_ehancement = conv_data-data

plot_x = np.arange(0,size,1)
center_x = size//2

mul = 1.5
step = 1
nrows = len(kernels)*step

# ----------------------------- #
# cmap='viridis'
cmap = cm.get_cmap('seismic')
cmap.set_bad(color='tab:green')

interval = MinMaxInterval()
_,vmax = interval.get_limits(data)
vmin,_ = interval.get_limits(conv_ehancement[-1])
norm = mcolors.CenteredNorm(halfrange=cmap_max)
if add_bkg:
    norm = mcolors.CenteredNorm(halfrange=cmap_max)

# ----------------------------- #

fig, axs = plt.subplots(nrows=nrows,ncols=5,figsize=(5*mul,nrows*mul),sharex=True,constrained_layout=True)
for i in range(0,nrows,step):
    data_idx = i//step
    
    axs[i][0].imshow(data,origin='lower',cmap=cmap,norm=norm)
    axs[i][0].set_ylabel(list(kernels.keys())[data_idx])
    axs[i][1].imshow(kernels_data_padded[data_idx],origin='lower',cmap=cmap,norm=norm)
    axs[i][1].text(1,1,f"{kernels_data[data_idx].shape}")
    
    if add_bkg:
        conv_thresh = bkg_rms*thresh_factor
        conv_data_copy = deepcopy(conv_data[data_idx])
        conv_data_copy[conv_data_copy >= conv_thresh] = np.nan
        axs[i][2].imshow(conv_data_copy,origin='lower',cmap=cmap,norm=norm)
    else:
        axs[i][2].imshow(conv_data[data_idx],origin='lower',cmap=cmap,norm=norm)
    axs[i][3].imshow(conv_ehancement[data_idx],origin='lower',cmap=cmap,norm=norm)
    if data_idx > 0:
        axs[i][4].imshow(conv_data[data_idx]-conv_data[data_idx-1],origin='lower',cmap=cmap,norm=norm)

    if i == len(kernels)-1:
        axs[i][0].set_xlabel('Collection of\ndistributions')
        axs[i][1].set_xlabel('convolution\nkernel')
        axs[i][2].set_xlabel('convolved\ndata')
        axs[i][3].set_xlabel('convolution\nenhancement')
        axs[i][4].set_xlabel('convolved data\nresidual')

    for ax in axs[i]:
        ax.set_xticks([])
        ax.set_yticks([])
        
[fig.delaxes(ax) for ax in axs.flatten() if not ax.has_data()]
im = cm.ScalarMappable(cmap=cmap,norm=norm)
fig.colorbar(im, ax=axs.ravel().tolist(), shrink=0.3)
fig.suptitle(f'Initial convolution kernel for source detection'f'\n'fr'(threshold: {cmap_max_factor} times the min. amplitude)')
if add_bkg:
    fig.suptitle(f'Initial convolution kernel for source detection'f'\n'fr'(threshold: $\geq$ {1+min_factor/100} times min. value for 8-connectedness of faintest source)')

[larrybradley]

@astqx You can input your convolved image into any of the source detection tools (detect_sources ,deblend_sources, SourceFinder). It's completely up to you how you want to make that convolved image. The make_2dgaussian_kernel is a simple convenience function for making a normalized 2D Gaussian kernel based on FWHM. This is by far the most common choice. If you want to use something else, there is a large list of pre-defined kernels in astropy.convolution: https://docs.astropy.org/en/stable/convolution/kernels.html#available-kernels. You can also use your own custom kernel. Just convolve the data with your kernel (convolved_data = convolve(data, kernel)) and pass the convolved image into the source detection tools. As you've noted, you may need adjust the detection threshold based on your convolution kernel.