Crop a (DASK) array of (tomographic) data to the biggest blob present in it.

cropper.py

def cropper(image,
            despeckle=True,
            verbose=False,
            threshold=50):
    '''Crop array to biggest item in it'''
    dimensions = len(image.shape)
    if verbose:
        print('Cropping %s-dimensional image' % dimensions)
    # Threshold
    thresholded = image > threshold
    if despeckle:
        if verbose:
            print('Removing small objects')
        despeckled = skimage.util.apply_parallel(skimage.morphology.remove_small_objects,
                                                 thresholded,
                                                 extra_keywords={'min_size': 10**dimensions},
                                                 chunks=[15,15,15])
    # Find biggest object
    # https://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.find_objects.html
    if verbose:
        print('Finding objects')
    if despeckle:
        cropdimensions = scipy.ndimage.find_objects(despeckled)[0]
    else:
        cropdimensions = scipy.ndimage.find_objects(thresholded)[0]
    if verbose:
        print('Cutting the image down to')
        for c, cd in enumerate(cropdimensions):
            print('\t- %s: %s' % (c, cd))
        if dimensions > 2:
            # Calculate cardinal direction MIPs
            MIPs = [image.max(axis=d) for d, direction in notebook.tqdm(enumerate(directions),
                                                                                  desc='Calculating MIPs',
                                                                                  total=len(directions))]
            for d, direction in enumerate(directions):
                plt.subplot(1, dimensions, d + 1)
                plt.imshow(MIPs[d], alpha=0.618)
                plt.imshow(dask.array.ma.masked_less(MIPs[d]>threshold, 1),
                           alpha=0.309,
                           cmap='viridis_r')
                if d == 0:
                    plt.axhline(cropdimensions[1].start,
                                c=seaborn.color_palette()[0],
                                label=cropdimensions[1].start)
                    plt.axhline(cropdimensions[1].stop,
                                c=seaborn.color_palette()[1],
                                label=cropdimensions[1].stop)
                    plt.axvline(cropdimensions[2].start,
                                c=seaborn.color_palette()[2],
                                label=cropdimensions[2].start)
                    plt.axvline(cropdimensions[2].stop,
                                c=seaborn.color_palette()[3],
                                label=cropdimensions[2].stop)
                elif d == 1:
                    plt.axhline(cropdimensions[0].start,
                                c=seaborn.color_palette()[0],
                                label=cropdimensions[0].start)
                    plt.axhline(cropdimensions[0].stop,
                                c=seaborn.color_palette()[1],
                                label=cropdimensions[0].stop)
                    plt.axvline(cropdimensions[2].start,
                                c=seaborn.color_palette()[2],
                                label=cropdimensions[2].start)
                    plt.axvline(cropdimensions[2].stop,
                                c=seaborn.color_palette()[3],
                                label=cropdimensions[2].stop)                    
                elif d == 2:
                    plt.axhline(cropdimensions[0].start,
                                c=seaborn.color_palette()[0],
                                label=cropdimensions[0].start)
                    plt.axhline(cropdimensions[0].stop,
                                c=seaborn.color_palette()[1],
                                label=cropdimensions[0].stop)
                    plt.axvline(cropdimensions[1].start,
                                c=seaborn.color_palette()[2],
                                label=cropdimensions[1].start)
                    plt.axvline(cropdimensions[1].stop,
                                c=seaborn.color_palette()[3],
                                label=cropdimensions[1].stop)                                        
                plt.title('%s MIP' % direction)
                plt.legend(loc='lower right')
        else:
            plt.subplot(121)
            plt.imshow(image, alpha=0.618)
            plt.imshow(dask.array.ma.masked_less(thresholded, 1), alpha=0.618, cmap='viridis_r')
            plt.axhline(cropdimensions[0].start, c=seaborn.color_palette()[0], label=cropdimensions[0].start)
            plt.axhline(cropdimensions[0].stop, c=seaborn.color_palette()[1], label=cropdimensions[0].stop)
            plt.axvline(cropdimensions[1].start, c=seaborn.color_palette()[2], label=cropdimensions[1].start)
            plt.axvline(cropdimensions[1].stop, c=seaborn.color_palette()[3], label=cropdimensions[1].stop)
            plt.title('Original image with thresholded overlay')
            plt.legend(loc='lower right')
            plt.subplot(122)
            plt.imshow(image[cropdimensions])
            plt.title('Output')
        plt.show()
    if verbose:
        print('Cropped image from %s to %s' % (image.shape, image[cropdimensions].shape))
        print(cropdimensions)
    return(image[cropdimensions])