Add support for using center of intensity as cell centers

[matham]

Description

What is this PR

• Bug fix
• Addition of a new feature
• Other

Why is this PR needed?

This is the final feature I had to add to make cellfinding work well for our images. In our FosTRAP images we encountered an issue where the location of some detected cells were significantly offset from the actual cell center. Looking closely with the help of the benchmark scripts it seemed to be a problem with cells that are a bit comet shaped. Specifically, these are bright cells that have other less bright neurites connected to it.

Because we have a range of intensities of our cells, the edge detection threshold is on the lower side. So if we have a cell with a very bright sphere, but connected to it is tissue that is much less bright but still above threshold, then all the bright voxels are considered equally when calculating the center of the cell. This will place the center to the side of the sphere. The solution is to either set the threshold much higher, or for us, it's to use the center of intensity when calculating the center so it'll be pulled towards the brightest voxels.

The following is an example cell from a test I added that tests the issue. We have a sphere with radius 7 whose center is located in z, y, x at (10, 25, 15). The center intensity is 1000 with a dropoff to 100 on the sphere's surface. But on the x-axis there's a bright region connected to it that drops off further to 10 and stops at x=37. Given we set the threshold low, all the non-zero voxels on the edges are considered bright, even though those two orders of magnitude lower than the sphere center. So the cell center using all the voxels would just be at (10, 25, 22). But ideally the location would be close to the sphere center.

The following image is the volume after the 3d filtering step. You can see there's a hole in the center where the gradient is very low outside the sphere, but everywhere else the gradient is above the threshold and so all those white voxels would be used to calculate the center of the cell.

What does this PR do?

This adds an option, detect_centre_of_intensity that is off by default. But when enabled, it weighs the voxels according to the intensity such that for the example above it correctly locates the center close to the sphere center. I also added full test coverage and tests that specifically test such comet shapes with the option ON and OFF.

To be able to do that I had to modify the filtering code to pass along the original input planes so we can get the original voxel intensity for this purpose. Similarly for structure splitting we pass along the intensity of the voxels in the structure being split. The additional memory and time cost should be minimal relative to everything else.

References

None.

How has this PR been tested?

I added thorough test coverage and also have been using with all of my samples.

I used the example above in the test. But additionally, the above example drops the intensity along the tail from 100 to 10, resulting in low gradients in the middle and the "hole" at those low gradients (i.e. below threshold) as shown above. To test the structure splitting code, I couldn't use this comet because the hole makes it split into two cells. I wanted an example comet that after splitting looks like just before so that we can test the center calculation of the splitting code after "splitting". The following example has similar linear gradients throughout because the center intensity is 1000 with a dropoff to 700 on the sphere's surface and 50 at the end of the tail where it stops. Now all gradients are above threshold and there's no whole so after 1 splitting iteration we end up with the same volume as the input except some dilation on the edges. And this lets us verify that splitting also correctly works with the feature ON or OFF:

Is this a breaking change?

No, by default it's off. Although I think we should consider turning this ON by default because this is probably a desired feature for everyone?

Does this PR require an update to the documentation?

There are new parameters.

Checklist:

• The code has been tested locally
• Tests have been added to cover all new functionality (unit & integration)
• The documentation has been updated to reflect any changes
• The code has been formatted with pre-commit

[matham]

The failed test is just a network issue and needs to be rerun.

I did forget to mention that I think

cellfinder/cellfinder/core/detect/filters/volume/structure_splitting.py

Lines 199 to 203 in 13274b0

	if (relative_coords > max_coords).all():
	logger.info(
	'Relative coordinates "{}" exceed maximum volume '
	'dimension of "{}"'.format(relative_coords, max_coords)
	)

is a bug. It should be any, not all. And also I don't think the logging is useful - it should either raise an error or handle it quietly except if we want to log to debugging. But I'm not quite sure how the point can ever be outside the valid cube if a valid position is passed in.

For now I added a test for that code that checks that it does the current behavior so if (when) we fix the behavior the test will also need to be updated - as it should be.

[matham]

As per my last comment, I did discover a bug that I fixed in this branch in the last commit.

Basically, after splitting, we checked each resulting structure to see if the center falls within the cuboid. I observed a bug in this check in that comment, but also wondered why it would ever be outside (other than rounding errors, which can just be forced to the edges). In fact there's this parameter about whether to allow these outliers or drop them.

Looking more closely at this because I was testing splitting, turns out there is a bug that resulted in these outliers. But once fixed, there's no possibility of outliers so we don't have to check for it. Specifically, given a cuboid of the points to be split, previously we padded the cuboid to accommodate the full ball filter. But, that padding and resulting offset due to this padding was incorrectly calculated and partially ignored. So the resulting split cells were not only shifted (which I had noticed visually, but attributed to "stuff"), but could sometimes fall outside the original unpadded cuboid.

The last commit fixes it all to use the correct padding and to compensate for the offset. There's even a test where in previous PRs where I added this comment # for some reason, same with pytorch, it's shifted by 1. Probably rounding about split cells. Now this if fixed and not shifted anymore.

I stuck it on to this branch, because I was working on this branch when I realized the issue. And I didn't want to make another PR and then have to later modify this PR on top of that, because they both work on the same code. But I can separate this fix if wanted.

Also, I didn't remove the outlier_keep parameter to detection, even though now it doesn't do anything, since there are no outliers anymore. Perhaps they can be removed later 🤷