Last active
September 28, 2020 15:03
-
-
Save habi/16ce36a93f66b503ddd0034166d598cc to your computer and use it in GitHub Desktop.
Function to crop a dataset to the smallest cuboid that contains (denoised) data. Displays the data if requested.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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)) | |
return(image[cropdimensions]) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment