After some experimentation, it seems that binary_closing is the best helper. And it is important to vary the structure, to avoid artifacts with "square corners"

import scipy.ndimage as nd

bones = img_data > args.threshold
n = 5
for t in range(2):
    nd.binary_closing(bones, iterations=1, output=bones, structure=np.ones((n,n,n)))
    nd.binary_closing(bones, iterations=1, output=bones)
nd.binary_opening(bones, iterations=1, output=bones)

bone_labels, num_feat = nd.label(bones)
vals, counts = np.unique(bone_labels, return_counts=True)
vals = vals[1:]
counts = counts[1:]
print('begin erase small labels')
for i,value in enumerate(vals):
    if counts[i] < args.min_size:
        bones[bone_labels == value] = False

The bones is a binary mask. The last part labels all components and then aggressively deletes the small ones. The min_size is 1500 in my case.