Hi Matteo,

Can you try putting this notebook up as a gist and pasting a link to the notebook? It's hard for me to follow all of the steps (and the polarity of the image) without the images inline. Is it just the inverse of what you want? And anyway why aren't you just using ndimage's binary_fill_holes?

https://docs.scipy.org/doc/scipy-0.15.1/reference/generated/scipy.ndimage.morphology.binary_fill_holes.html

Juan.

On Fri, Mar 27, 2015 at 9:09 AM, Matteo <matteo....@gmail.com> wrote:

Hello Juan

Thanks so much for your suggestions.

Once I convertedthe image as you suggested:

# import back image
cfthdr=io.imread('filled_contour_THDR.png')
cfthdr = color.rgb2gray(cfthdr) > 0.5

I get good results with opening:

# clean it up with opening
selem17 = disk(17)
opened_thdr = opening(cfthdr, selem17)/255
# plot it
fig = plt.figure(figsize=(5, 5))
ax = fig.add_subplot(1, 1, 1)
ax.set_xticks([])
ax.set_yticks([])
plt.imshow(opened_thdr,cmap='bone')
plt.show()
# not bad

With remove_small_objects the advantage is that it does not join blobs in the original:

cfthdr_inv = ~cfthdr
test=remove_small_objects(cfthdr,10000)
# plot it
fig = plt.figure(figsize=(5, 5))
ax = fig.add_subplot(1, 1, 1)
ax.set_xticks([])
ax.set_yticks([])
plt.imshow(test,cmap='bone')
plt.show()

but with reconstruction done as this:

# filling holes with morphological reconstruction
seed = np.copy(cfthdr_inv)
seed[1:-1, 1:-1] = cfthdr_inv.max()
mask = cfthdr_inv
filled = reconstruction(seed, mask, method='erosion')
# plot it
fig = plt.figure(figsize=(5, 5))
ax = fig.add_subplot(1, 1, 1)
ax.set_xticks([])
ax.set_yticks([])
plt.imshow(filled,cmap='bone',vmin=cfthdr_inv.min(), vmax=cfthdr_inv.max())
plt.show()

I get a completely white image. Do you have any suggestions as to why?

Thank again. Cheers,

Matteo

On Wednesday, March 25, 2015 at 6:29:43 PM UTC-6, Juan Nunez-Iglesias wrote:

Hi Matteo,

My guess is that even though you are looking at a "black and white" image, the png is actually an RGB png. Just check the output of "print(cfthdr.shape)". Should be straightforward to make it a binary image:

from skimage import color

cfthdr = color.rgb2gray(cfthdr) > 0.5

Then you should have a binary image. (And inverting should be as simple as "cfthdr_inv = ~cfthdr") We have morphology.binary_fill_holes to do what you want.

btw, there's also morphology.remove_small_objects, which does exactly what you did but as a function call. Finally, it looks like you are not using the latest version of scikit-image (0.11), so you might want to upgrade.

Hope that helps!

Juan.

On Thu, Mar 26, 2015 at 8:48 AM, Matteo <matteo....@gmail.com> wrote:

Issues with morphological filters when trying to remove white holes in black objects in a binary images. Using opening or filling holes on inverted (or complement) of the original binary.

Hi there

I have a series of derivatives calculated on geophysical data.

Many of these derivatives have nice continuous maxima, so I treat them as images on which I do some cleanup with morphological filter.

Here's one example of operations that I do routinely, and successfully:

# threshold theta map  using Otsu method

thresh_th = threshold_otsu(theta)

binary_th = theta > thresh_th

# clean up small objects

label_objects_th, nb_labels_th = sp.ndimage.label(binary_th)

sizes_th = np.bincount(label_objects_th.ravel())

mask_sizes_th = sizes_th > 175

mask_sizes_th[0] = 0

binary_cleaned_th = mask_sizes_th[label_objects_th]

# further enhance with morphological closing (dilation followed by an erosion) to remove small dark spots and connect small bright cracks

# followed by an extra erosion

selem = disk(1)

closed_th = closing(binary_cleaned_th, selem)/255

eroded_th = erosion(closed_th, selem)/255

# Finally, extract lienaments using skeletonization

skeleton_th=skeletonize(binary_th)

skeleton_cleaned_th=skeletonize(binary_cleaned_th)

# plot to compare

fig = plt.figure(figsize=(20, 7))

ax = fig.add_subplot(1, 2, 1)

imshow(skeleton_th, cmap='bone_r', interpolation='none')

ax2 = fig.add_subplot(1, 3, 2)

imshow(skeleton_cleaned_th, cmap='bone_r', interpolation='none')

ax.set_xticks([])

ax.set_yticks([])

ax2.set_xticks([])

ax2.set_yticks([])

Unfortunately I cannot share the data as it is proprietary, but I will for the next example, which is the one that does not work.

There's one derivative that shows lots of detail but not continuous maxima. As a workaround I created filled contours in Matplotlib

exported as an image. The image is attached.

Now I want to import back the image and plot it to test:

# import back image

cfthdr=io.imread('filled_contour.png')

# threshold using using Otsu method

thresh_thdr = threshold_otsu(cfthdr)

binary_thdr = cfthdr > thresh_thdr

# plot it

fig = plt.figure(figsize=(5, 5))

ax = fig.add_subplot(1, 1, 1)

ax.set_xticks([])

ax.set_yticks([])

plt.imshow(binary_thdr, cmap='bone')

plt.show()

The above works without issues.

 

Next I want to fill the white holes inside the black blobs. I thought of 2 strategies.

The first would be to use opening; the second to invert the image, and then fill the holes as in here:

http://scikit-image.org/docs/dev/auto_examples/plot_holes_and_peaks.html

By the way, I found a similar example for opencv here

http://stackoverflow.com/questions/10316057/filling-holes-inside-a-binary-object

 

Let's start with opening. When I try:

selem = disk(1)

opened_thdr = opening(binary_thdr, selem)

or:

selem = disk(1)

opened_thdr = opening(cfthdr, selem)

I get an error message like this:

---------------------------------------------------------------------------

ValueError                                Traceback (most recent call last)

<ipython-input-49-edc0d01ba327> in <module>()

      1 #binary_thdr=img_as_float(binary_thdr,force_copy=False)

----> 2 opened_thdr = opening(binary_thdr, selem)/255

      3

...