Allow passing one value per axis of the input image as the min_distance parameter of skimage.feature.peak_local_max

[ElpadoCan]

Description:

As discussed in #7317, it would be cool to have the possibility to detect peaks in anisotropic images with a min_distance made of one value per axis of the input image.

I will start to look into the code to see whether I can implement it myself, but if someone wants to help or has ideas feel free to step in.

[jni]

If I remember correctly, we use scipy.spatial.cKDTree to implement min_distance, and that doesn't support anisotropic distance checks. However, you could hack it by doing something like:

• diving all coordinates by the spacing
• using 1 as the min-distance
• multiplying the final coordinates back by the spacing

I think the API is a bit tricky... It probably doesn't make sense to use different min_distance per axis because these are not manhattan distances. Rather, we could add a spacing= keyword argument, and the (single) distance is computed relative to that spacing, using the method above.

But either way, I think this would indeed be a nice feature to add. 😊

[ElpadoCan]

Hi @jni,

Interesting idea! I like the spacing= approach. I saw the use of scipy.spatial.cKDTree, but I implemented something different for my library where I use masked arrays (where you don't check the distance from/to already dropped points) and it checks whether multiple points are lying within the same ellipsoid. The implementation is quite short:

masked_points = np.ma.masked_array(
    data=points, 
    mask=np.zeros(points.shape, bool), 
    fill_value=-1
)
valid_points_mask = np.ones(num_points, dtype=bool)
for i, point in enumerate(points):
    if not valid_points_mask[i]:
        # Skip points that have already been dropped
        continue
    points_ellipsoid = np.square((masked_points - point)/min_distance)
    points_too_close_mask = np.sum(points_ellipsoid, axis=1) < 1
    valid_points_mask[i+1:] = np.invert(points_too_close_mask[i+1:])
    
    # Mask dropped points to avoid computing the distance to them in future iterations
    masked_points.mask[points_too_close_mask] = True
    
valid_points = points[valid_points_mask]

The idea is that we could sort the points coordinates by intensities and my implementation would keep only the brightest peak within the same ellipsoid. What do you think?

The problem I see is that I'm using Euclidean distances but I think there is a way to use all distances like the API does.

Would your implementation keep only the brightest peak if multiple peaks are within the minimum distance?

[lagru]

Note that there is also #4165 which has been in the work for quite some time. It also uses cKDTree so the same limitations to how distance is calculated applies. It has been already working and should also work on peaks as returned by peak_local_max. However, the approach ignores anisotropic images for now.

I think it's mainly been stalled due to discussion around definitions of priority and performance optimization and reviews in general. But I can and want to revisit it, after some other priority items are addressed.

[lagru]

I added support for anisotropic data to #4165. That function can be run on the labeled output of peak_local_max to ensure a correct distance.

Additionally, we could try to implement a similar approach for ensure_spacing which is used by peak_local_max.

[ElpadoCan]

Fantastic @lagru, that sounds great thanks!

[ElpadoCan]

@lagru I have a question about ensure_spacing. If two points are too close, which one is kept? Because in peak_local_max it would be nice to keep the coordinate of the higher peak, right?

[lagru]

I think higher peaks take precedence. The coordinates passed to ensure_spacing are sorted by intensity in

scikit-image/skimage/feature/peak.py

Lines 16 to 19 in 9b91d3a

	intensities = image[coord]
	# Highest peak first
	idx_maxsort = np.argsort(-intensities, kind="stable")
	coord = np.transpose(coord)[idx_maxsort]

The "-" makes highest value go first. These are then passed to the kd-tree in batches and iterated in the same order:

scikit-image/skimage/_shared/coord.py

Lines 31 to 37 in 9b91d3a

	tree = cKDTree(coord)
	indices = tree.query_ball_point(coord, r=spacing, p=p_norm)
	rejected_peaks_indices = set()
	naccepted = 0
	for idx, candidates in enumerate(indices):
	if idx not in rejected_peaks_indices:

--

remove_near_objects in #4165, allows sorting by an optional priority parameter.