Weird/Incorrect results from SimpleITK:ConnectedComponent()

jbi35 · April 10, 2020, 6:20pm

Hi there,
I’ve come across some unexpected behavior while using the procedural interface of the ConnectedComponentImageFilter in SimpleITK. When applying the filter in order to extract the largest connected components in a segmentation, I get only 5 connected components, while I expect hundreds or thousands. I’ve also doubled checked this by computing the connected components in Paraview.

In order to reproduce this, I’ve uploaded a version of the segmentation I am working on here, along with the results I get from Paraview and a minimal example of the code to reproduce this.

Blockquote
import SimpleITK as sitk
import numpy as np

# create some dummy checkerboard data 
data = np.kron([[1, 0] * 20, [0, 1] *20] * 20, np.ones((10, 10, 10)))
some_zeros = np.zeros((10, 400, 400))
data = np.vstack((data, some_zeros, data,
                  data, some_zeros, data,
                  data, some_zeros, data,
                  data, some_zeros, data,
                  data, some_zeros, data,
                  data, some_zeros, data,
                  data, some_zeros, data,
                  data, some_zeros, data,
                  data, some_zeros, data,
                  data, some_zeros, data,
                  data, some_zeros, data,
                  data, some_zeros, data))

conn_comps = sitk.GetImageFromArray(data)
conn_comps = sitk.Cast(conn_comps, sitk.sitkInt8)

#print("Dummy Image before {}".format(conn_compsconn_comps))
conn_comps = sitk.ConnectedComponent(conn_comps)
#print("Dummy Image after {}".format(conn_comps))

sitk.WriteImage(conn_comps, 'conn_comps.mhd')


segmentation =sitk.ReadImage('segmentation.mhd')
# some cleanup
segmentation = sitk.BinaryOpeningByReconstruction(segmentation, 5, sitk.sitkBall)

fully_connected_neighborhood = False  # face-connected neighbors only
# conn_comps = sitk.ScalarConnectedComponent(result, result, 10, fully_connected_neighborhood)
conn_comps = sitk.ConnectedComponent(segmentation, fully_connected_neighborhood)

sitk.WriteImage(conn_comps, 'conn_comps_2.mhd')

Blockquote

The upper dummy example creates a checkerboard pattern, for which the connected component filter produces the expected results. However, for the segmentation, it does not and I cannot figure out why.
I’ve tried v1.2.4 as well as SimpleITK 2.0 Release Candidate 1.

zivy · April 10, 2020, 9:10pm

Hello @jbi35,

I believe your code is doing what you want, you are making sure that you don’t have those hundreds/thousands of spurious connected components by calling the BinaryOpeningByReconstruction. If you don’t run this cleanup step you will get the large number of spurious components (in the thousands).

Please clarify, are you expecting to still get all these spurious components after you cleaned things up?

jbi35 · April 11, 2020, 6:09am

Hello @zivy,

sorry for being not precise enough. I probably also should have removed the call to the BinaryOpeningByReconstruction as this, I believe not the issue here.
IMHO the problem is that the ConnectComponentImageFilter does not correctly separate clearly non-connected components. When running the filter I get only 5 components while I should get many more.
When visually inspecting the result, it is apparent that the majority of the voxels within one component are not connected to each other. Nevertheless, the ConnectedComponentImageFilter labels them as such. I’ve attached a screenshot where the components are labeled by color.

Unfortunately, I was so far unable to reproduce this behavior within a simpler example and hence have tried to provide all the files necessary to reproduce this.

dzenanz · April 11, 2020, 1:41pm

Which version of SimpleITK are you using? A bug in connected components filter was recently fixed. But I don’t think that bug can quite explain the problem in your screenshot.

jbi35 · April 11, 2020, 1:46pm

Yes, I saw this and hence I tried the new release candidate of SimpleITK 2. Before that I worked with version 1.2.4

blowekamp · April 11, 2020, 9:03pm

This bug was introduced with the threading of ITK 5.0, and fix in 5.1. SimpleITK 1.2.x was build against ITK4.13, and SimpleITK 2.0rc1is build against ITK 5.1rc2. So this bug should have been missed with the SimpleITK binary releases.

zivy · April 12, 2020, 4:35pm

Hello @jbi35,

Still not sure about the problem. I ran the following code and got what I expected, without the cleanup line we get 12679 connected components. With the cleanup line we get 5 components and they do look like they are connected correctly (used volume rendering). Please take a look and confirm the problem, recreate the rendering for the saved volume after cleanup.

import SimpleITK as sitk
import numpy as np

segmentation =sitk.ReadImage('segmentation.mhd')

fully_connected_neighborhood = False # face-connected neighbors only
conn_comps = sitk.ConnectedComponent(segmentation, fully_connected_neighborhood)
unique_comps = np.unique(sitk.GetArrayViewFromImage(conn_comps))
print(len(unique_comps))
print(unique_comps)
sitk.WriteImage(conn_comps, 'segmentationWithoutCleanup.mhd')

# Same thing with cleanup
segmentation = sitk.BinaryOpeningByReconstruction(segmentation, 5, sitk.sitkBall)

fully_connected_neighborhood = False # face-connected neighbors only
conn_comps = sitk.ConnectedComponent(segmentation, fully_connected_neighborhood)
unique_comps = np.unique(sitk.GetArrayViewFromImage(conn_comps))
print(len(unique_comps))
print(unique_comps)
sitk.WriteImage(conn_comps, 'segmentationWithCleanup.mhd')

jbi35 · April 14, 2020, 4:09pm

Dear @zivy,

again, I am sorry for the confusion. The cleanup step is not the problem. I ran your code and as you said, without cleanup I get 12679 components. The problem is that these are too few and there are unconnected parts lumped together into one component. E.g. when looking at
segmentationWithoutCleanup.mhd in Paraview, I get this:

Then, looking only at the first component of this, see code below, I obtain this:

There are lots of unconnected voxels as can be seen when applying Paraviews Connectivity filter (just to this component) which yields this (I get additional 846 components in Paraview just for component 1):

So in my opinion, the ConnectedComponent filter lumps several unconnected components together, or am I wrong?

import SimpleITK as sitk

import numpy as np
segmentation =sitk.ReadImage('segmentation.mhd')
fully_connected_neighborhood = False # face-connected neighbors only
conn_comps = sitk.ConnectedComponent(segmentation, fully_connected_neighborhood)
unique_comps = np.unique(sitk.GetArrayViewFromImage(conn_comps))

print("Number of components SITK {}".format(len(unique_comps)))
print(unique_comps)

sitk.WriteImage(conn_comps, 'segmentationWithoutCleanup.mhd')
# get only component 1
result = sitk.BinaryThreshold(conn_comps, 1, 1.99, 1, 0)
conn_comps = sitk.ConnectedComponent(result, fully_connected_neighborhood)
unique_comps = np.unique(sitk.GetArrayViewFromImage(conn_comps))
sitk.WriteImage(result, 'segmentationWithoutCleanupComp1.mhd')

zivy · April 14, 2020, 5:40pm

Hello @jbi35,

This looks to me more of an issue with the visualization than with the segmentation. I tried to replicate what you were doing using ITK-SNAP (loaded segmentation.mhd as volume and the first component, sitk.WriteImage(conn_comps==1,'compOne.mhd'), as the segmentation. I zoomed onto a 3D component (bottom left view) which looks like it is separate. ITK-SNAP’s cursor is linked and you can see on the top right that the voxels are face connected. I suspect that what we are seeing is significant leakage when there are a small number of voxels with face connections bridging between large connected components. A single such voxel is sufficient and it is all but impossible to visually detect it in 3D. If a small binary morphology is enough to separate things then I trust the filter is working correctly. We can only verify it using artificial test images where it does work, at least with the ones used by ITK. If you can create a simple image that shows that the filter does not respect the face based connectivity then you have found a bug, otherwise it is more likely that the face connectivity is respected and we are seeing leakage.

jbi35 · May 3, 2020, 2:11pm

Dear @zivy,

thanks for your input and sorry for getting back to you so late. It was indeed a visualization issue, as you suspected. The problem was looking at the data in Paraview, which did not correctly display the data. was now able to fix the visualization problem and agree with you that the filter is working just fine. Anyways, thanks for your time and help.