apply spacing rescaling when computing centroid_weighted

[enricotagliavini]

Description

We noticed the numbers for the centroid_weighted returned by regionprops, when using anisotropic spacing (e.g. 2, 1 1), where not correct. Upon investigation it looks like the slice.start coordinates are not rescaled before adding to the, correctly rescaled, centroid_weighted_local. This commit fixes the issue and should now return the correct coordinates.

For reviewers

• Check that the PR title is short, concise, and will make sense 1 year
  later.
• Check that new functions are imported in corresponding __init__.py.
• Check that new features, API changes, and deprecations are mentioned in
  doc/release/release_dev.rst.
• There is a bot to help automate backporting a PR to an older branch. For
  example, to backport to v0.19.x after merging, add the following in a PR
  comment: @meeseeksdev backport to v0.19.x
• To run benchmarks on a PR, add the run-benchmark label. To rerun, the label
  can be removed and then added again. The benchmark output can be checked in
  the "Actions" tab.

[jni]

Thanks for the PR @enricotagliavini! And sorry about that oversight 😅

It looks like some tests are failing because of it:

https://github.com/scikit-image/scikit-image/actions/runs/4754448065/jobs/8447280174?pr=6900#step:5:1036

FAILED measure/tests/test_regionprops.py::test_multichannel_centroid_weighted_table - AssertionError: 
Arrays are not almost equal to 7 decimals

(shapes (2,), (2, 2) mismatch)
 x: array([5.5405405, 9.4459459])
 y: array([[5.5405405, 5.5405405],
       [9.4459459, 9.4459459]])
FAILED measure/tests/test_regionprops.py::test_multichannel - ValueError: operands could not be broadcast together with shapes (2,) (2,3)

I think probably there is a broadcasting/shape error in how you are multiplying the centroid and the spacing, particularly in the case of multichannel images? Could (a) double check that, and (b) also add a test to make sure that this fix is captured by tests in the future? Thank you!

[enricotagliavini]

Yeah I noticed. I had a quick look but didn't figure it out yet. I'll try to have a look in the coming weeks, hopefully I can figure it out.

[enricotagliavini]

Ok this passed the tests. It's a bit odd as the columns are the coordinate (e.g. (X Y Z)) for single channel images, but it becomes the channel index for multi channel images. I didn't know this when doing the first patch and it was tested only on single channel images. Although we always have multichannel images, we always analyze the channels separately.

As for a test I can think of this: an image with only one non zero point in the corner furthest from the origin. With a spacing greater than 1 the centroid would end up outside the image, with the unfixed version of the code. Would that be acceptable?

[enricotagliavini]

Actually you already have a test for this, but the input image being used is not general enough to cover all cases. In particular for that image the slice.start is always the image origin, so it's 0. Even if not correctly multiplied by the spacing the end result is always correct as any number multiplied by zero yields 0.

[lagru]

Thanks for working on this @enricotagliavini! If you merge the current main into this branch, the failing job on CircelCI should fix itself.

[lagru]

I'm probably not grasping all the bits and pieces properly yet but I'm confused by something. I checked how the spacing is taken into account in the non-weighted case in

scikit-image/skimage/measure/_regionprops.py

Lines 433 to 438 in 2f9b82c

	def coords_scaled(self):
	indices = np.argwhere(self.image)
	object_offset = np.array([
	self.slice[i].start for i in range(self._ndim)
	])
	return (object_offset + indices) * self._spacing + self._offset

(called by self.centroid). There it looks like the spacing is multiplied with the sum of index and offset / slice.start, while in this fix it's only multiplied with slice.start. Should it be

Suggested change

	return tuple(idx + slc.start * spc
	for idx, slc, spc in zip(ctr, self.slice, self._spacing))
	return tuple((idx + slc.start) * spc
	for idx, slc, spc in zip(ctr, self.slice, self._spacing))

instead?

[enricotagliavini]

The non-weighted centroid is computed in a much different way, as it's much simpler. It computes the centroid without taking into account spacing, and applies all the corrections at the very end.

When computing the weighted centroid we start from the local (to the prop bbox) centroid, which is already corrected for spacing, and we just need to apply a translation to the coordinate system of the original (full size) image. If we apply the correction again as suggested this would yield wrong coordinates, as the local centroid would be corrected twice.

More in details: the centroid weighted local uses the weighted moments

scikit-image/skimage/measure/_regionprops.py

Lines 638 to 660 in 2f9b82c

	def centroid_weighted_local(self):
	M = self.moments_weighted
	M0 = M[(0,) * self._ndim]
	def _get_element(axis):
	return (0,) * axis + (1,) + (0,) * (self._ndim - 1 - axis)
	return np.asarray(
	tuple(M[_get_element(axis)] / M0 for axis in range(self._ndim)))
	@property
	@_cached
	def moments_weighted(self):
	image = self._image_intensity_double()
	if self._multichannel:
	moments = np.stack(
	[_moments.moments(image[..., i], order=3, spacing=self._spacing)
	for i in range(image.shape[-1])],
	axis=-1,
	)
	else:
	moments = _moments.moments(image, order=3, spacing=self._spacing)
	return moments

and the weighted moments just use the moments which take spacing into account (see line 261 in particular):

scikit-image/skimage/measure/_moments.py

Lines 154 to 269 in 2f9b82c

	def moments(image, order=3, *, spacing=None):
	"""Calculate all raw image moments up to a certain order.
	The following properties can be calculated from raw image moments:
	* Area as: ``M[0, 0]``.
	* Centroid as: {``M[1, 0] / M[0, 0]``, ``M[0, 1] / M[0, 0]``}.
	Note that raw moments are neither translation, scale nor rotation
	invariant.
	Parameters
	----------
	image : nD double or uint8 array
	Rasterized shape as image.
	order : int, optional
	Maximum order of moments. Default is 3.
	spacing: tuple of float, shape (ndim, )
	The pixel spacing along each axis of the image.
	Returns
	-------
	m : (``order + 1``, ``order + 1``) array
	Raw image moments.
	References
	----------
	.. [1] Wilhelm Burger, Mark Burge. Principles of Digital Image Processing:
	Core Algorithms. Springer-Verlag, London, 2009.
	.. [2] B. Jähne. Digital Image Processing. Springer-Verlag,
	Berlin-Heidelberg, 6. edition, 2005.
	.. [3] T. H. Reiss. Recognizing Planar Objects Using Invariant Image
	Features, from Lecture notes in computer science, p. 676. Springer,
	Berlin, 1993.
	.. [4] https://en.wikipedia.org/wiki/Image_moment
	Examples
	--------
	>>> image = np.zeros((20, 20), dtype=np.float64)
	>>> image[13:17, 13:17] = 1
	>>> M = moments(image)
	>>> centroid = (M[1, 0] / M[0, 0], M[0, 1] / M[0, 0])
	>>> centroid
	(14.5, 14.5)
	"""
	return moments_central(image, (0,) * image.ndim, order=order, spacing=spacing)
	def moments_central(image, center=None, order=3, *, spacing=None, **kwargs):
	"""Calculate all central image moments up to a certain order.
	The center coordinates (cr, cc) can be calculated from the raw moments as:
	{``M[1, 0] / M[0, 0]``, ``M[0, 1] / M[0, 0]``}.
	Note that central moments are translation invariant but not scale and
	rotation invariant.
	Parameters
	----------
	image : nD double or uint8 array
	Rasterized shape as image.
	center : tuple of float, optional
	Coordinates of the image centroid. This will be computed if it
	is not provided.
	order : int, optional
	The maximum order of moments computed.
	spacing: tuple of float, shape (ndim, )
	The pixel spacing along each axis of the image.
	Returns
	-------
	mu : (``order + 1``, ``order + 1``) array
	Central image moments.
	References
	----------
	.. [1] Wilhelm Burger, Mark Burge. Principles of Digital Image Processing:
	Core Algorithms. Springer-Verlag, London, 2009.
	.. [2] B. Jähne. Digital Image Processing. Springer-Verlag,
	Berlin-Heidelberg, 6. edition, 2005.
	.. [3] T. H. Reiss. Recognizing Planar Objects Using Invariant Image
	Features, from Lecture notes in computer science, p. 676. Springer,
	Berlin, 1993.
	.. [4] https://en.wikipedia.org/wiki/Image_moment
	Examples
	--------
	>>> image = np.zeros((20, 20), dtype=np.float64)
	>>> image[13:17, 13:17] = 1
	>>> M = moments(image)
	>>> centroid = (M[1, 0] / M[0, 0], M[0, 1] / M[0, 0])
	>>> moments_central(image, centroid)
	array([[16., 0., 20., 0.],
	[ 0., 0., 0., 0.],
	[20., 0., 25., 0.],
	[ 0., 0., 0., 0.]])
	"""
	if center is None:
	# Note: No need for an explicit call to centroid.
	# The centroid will be obtained from the raw moments.
	moments_raw = moments(image, order=order, spacing=spacing)
	return moments_raw_to_central(moments_raw)
	if spacing is None:
	spacing = np.ones(image.ndim)
	float_dtype = _supported_float_type(image.dtype)
	calc = image.astype(float_dtype, copy=False)
	for dim, dim_length in enumerate(image.shape):
	delta = (
	np.arange(dim_length, dtype=float_dtype) * spacing[dim] - center[dim]
	)
	powers_of_delta = (
	delta[:, np.newaxis] ** np.arange(order + 1, dtype=float_dtype)
	)
	calc = np.rollaxis(calc, dim, image.ndim)
	calc = np.dot(calc, powers_of_delta)
	calc = np.rollaxis(calc, -1, dim)
	return calc

Does this clear the concern?

[lagru]

Yes, I think I'm following now. Thanks for the detailed explanation and the contribution in total! :D

[enricotagliavini]

Thanks for working on this @enricotagliavini! If you merge the current main into this branch, the failing job on CircelCI should fix itself.

Thanks @lagru for checking this out and confirming this was a CI issue. I was quite puzzled for a bit. I addressed your comment on the fix, let me know if this clarifies.

[stefanv]

Thank you for your contribution, @enricotagliavini.

[enricotagliavini]

Thank you so much for accepting the fix. Cheers.