0

After reading an interesting topic on scipy.ndimage.label (Variable area threshold for identifying objects - python), I'd like to include an 'error margin' in the labelling.

In the above linked discussion: How can the blue dot on top be included, too (let's say it is wrongly disconnected from the orange, biggest, object)?

I found the structure attribute, which should be able to include that dot by changing the array (from np.ones(3,3,3) to anything more than that (I'd like it to be 3D). However, adjusting the 'structure' attribute to a larger array does not seem to work, unfortunately. It either gives an error of dimensions (RuntimeError: structure and input must have equal rank ) or it does not change anything..

Thanks!

this is the code:

labels, nshapes = ndimage.label(a, structure=np.ones((3,3,3)))

in which a is a 3D array.

Jules
  • 7
  • 5
  • The number of dimensions of `structure` must match the number of dimensions of `a`. The length of each dimension of `structure` must be exactly 3. So it looks like you can't use `structure` to do what you want. (This is a comment instead of an answer because the fundamental question is "How can the blue dot on top be included, too?") – Warren Weckesser Jul 18 '17 at 13:35
  • By they way, any simple approach to including an error margin would probably also merge the two shapes that are against the left edge. – Warren Weckesser Jul 18 '17 at 13:46

1 Answers1

0

Here's a possible approach that uses scipy.ndimage.binary_dilation. It is easier to see what is going on in a 2D example, but I'll show how to generalize to 3D at the end.

In [103]: a
Out[103]: 
array([[0, 0, 0, 1, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0],
       [1, 1, 0, 0, 1, 0, 0],
       [1, 1, 0, 0, 0, 1, 1],
       [0, 0, 0, 0, 0, 1, 1],
       [1, 1, 1, 0, 0, 0, 0]])

In [104]: from scipy.ndimage import label, binary_dilation

Extend each "shape" by one pixel down and to the right:

In [105]: b = binary_dilation(a, structure=np.array([[0, 0, 0], [0, 1, 1], [0, 1, 1]])).astype(int)

In [106]: b
Out[106]: 
array([[0, 0, 0, 1, 1, 0, 0],
       [0, 0, 0, 1, 1, 0, 0],
       [1, 1, 1, 0, 1, 1, 0],
       [1, 1, 1, 0, 1, 1, 1],
       [1, 1, 1, 0, 0, 1, 1],
       [1, 1, 1, 1, 0, 1, 1]])

Apply label to the padded array:

In [107]: labels, numlabels = label(b)

In [108]: numlabels
Out[108]: 2

In [109]: labels
Out[109]: 
array([[0, 0, 0, 1, 1, 0, 0],
       [0, 0, 0, 1, 1, 0, 0],
       [2, 2, 2, 0, 1, 1, 0],
       [2, 2, 2, 0, 1, 1, 1],
       [2, 2, 2, 0, 0, 1, 1],
       [2, 2, 2, 2, 0, 1, 1]], dtype=int32)

By multiplying a by labels, we get the desired array of labels of a:

In [110]: alab = labels*a

In [111]: alab
Out[111]: 
array([[0, 0, 0, 1, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0],
       [2, 2, 0, 0, 1, 0, 0],
       [2, 2, 0, 0, 0, 1, 1],
       [0, 0, 0, 0, 0, 1, 1],
       [2, 2, 2, 0, 0, 0, 0]])

(This assumes that the values in a are 0 or 1. If they are not, you can use alab = labels * (a > 0).)

For a 3D input, you have to change the structure argument to binary_dilation:

struct = np.zeros((3, 3, 3), dtype=int)
struct[1:, 1:, 1:] = 1
b = binary_dilation(a, structure=struct).astype(int)
Warren Weckesser
  • 110,654
  • 19
  • 194
  • 214
  • That is a clever approach, thank you! However, is there a reason you choose to not dilate [1, 1, 0] (above and left?) Would that merge the two shapes at the left edge, as you said already? Why? – Jules Jul 18 '17 at 14:15
  • By expanding just one side of each dimension, the error margin is one pixel. If you expand in all directions, you'll get an error margin of two pixels. – Warren Weckesser Jul 18 '17 at 14:25
  • Brilliant. And just for clarification: what is happening in the second line of the 3D add-on? – Jules Jul 18 '17 at 14:30
  • That fills the 2x2x2 "cube" in `struct` with 1s. It's a concise way of writing `struct[1,1,1] = 1; struct[2,1,1] = 1; struct[1, 2, 1] = 1; etc.`. – Warren Weckesser Jul 18 '17 at 14:37
  • Would you maybe have any ideas on how to threshold the minimal object size in 3D? Say I only wanted objects that are a minimal of np.ones((4,4,4)). – Jules Jul 19 '17 at 11:39
  • That should probably be a new SO question. I have ideas, but I don't have time at the moment to check that they work. If you make it a question, *someone* could probably help you sooner. – Warren Weckesser Jul 19 '17 at 14:09