Convallaria
Find me on Github Convallaria (c. majalis) acquired by Britta Schroth-Diez of the MPI-CBG Light Microscopy Facility, denoised using structN2V. The data consists of 100 images of the same field of view, in order to show structN2V denoising capabilities, we only select a single image as to not overestimate the results by providing too much information to the network. The ground truth is obtained by averaging all the images together.
In [1]:
Copied!
# Imports necessary to execute the code
from pathlib import Path
import matplotlib.pyplot as plt
import numpy as np
import pooch
import tifffile
from careamics import CAREamist
from careamics.config import create_n2v_configuration
from careamics.config.transformations import XYFlipModel
from careamics.utils import autocorrelation
from careamics.utils.metrics import scale_invariant_psnr
from PIL import Image
# Imports necessary to execute the code from pathlib import Path import matplotlib.pyplot as plt import numpy as np import pooch import tifffile from careamics import CAREamist from careamics.config import create_n2v_configuration from careamics.config.transformations import XYFlipModel from careamics.utils import autocorrelation from careamics.utils.metrics import scale_invariant_psnr from PIL import Image
Import the dataset¶
The dataset can be directly downloaded using the careamics-portfolio package, which uses pooch to download the data.
In [ ]:
Copied!
# # instantiate data portfolio manage
# portfolio = PortfolioManager()
# # and download the data
# root_path = Path("./data")
# files = portfolio.denoising.Flower.download(root_path)
# create a folder for our data.
# TODO add it to the portfolio
path = Path("data")
if not path.exists():
path.mkdir()
# download file
path_to_file = pooch.retrieve(
"https://download.fht.org/jug/n2v/flower.tif",
known_hash="08c2c77ae63d43d22631d69ae905ee695e8225cf9872183faf49bd0771e5b1a2",
path=path,
)
# # instantiate data portfolio manage # portfolio = PortfolioManager() # # and download the data # root_path = Path("./data") # files = portfolio.denoising.Flower.download(root_path) # create a folder for our data. # TODO add it to the portfolio path = Path("data") if not path.exists(): path.mkdir() # download file path_to_file = pooch.retrieve( "https://download.fht.org/jug/n2v/flower.tif", known_hash="08c2c77ae63d43d22631d69ae905ee695e8225cf9872183faf49bd0771e5b1a2", path=path, )
Visualize data and autocorrelation¶
In [3]:
Copied!
x_min, x_max = 580, 630
y_min, y_max = 380, 430
# load training and validation image and show them side by side
data = tifffile.imread(path_to_file)
print(f"Shape of the data {data.shape}")
train_image = data[50]
print(f"Shape of the training data {train_image.shape}")
gt_image = data.mean(axis=0)
print(f"Shape of the ground truth {gt_image.shape}")
# plot images
fig, ax = plt.subplots(2, 2, figsize=(10, 10))
ax[0, 0].imshow(gt_image, cmap="gray")
ax[0, 0].set_title("Ground truth")
ax[0, 1].imshow(train_image, cmap="gray")
ax[0, 1].set_title("Training image")
ax[1, 0].imshow(gt_image[x_min:x_max, y_min:y_max], cmap="gray")
ax[1, 1].imshow(train_image[x_min:x_max, y_min:y_max], cmap="gray")
x_min, x_max = 580, 630 y_min, y_max = 380, 430 # load training and validation image and show them side by side data = tifffile.imread(path_to_file) print(f"Shape of the data {data.shape}") train_image = data[50] print(f"Shape of the training data {train_image.shape}") gt_image = data.mean(axis=0) print(f"Shape of the ground truth {gt_image.shape}") # plot images fig, ax = plt.subplots(2, 2, figsize=(10, 10)) ax[0, 0].imshow(gt_image, cmap="gray") ax[0, 0].set_title("Ground truth") ax[0, 1].imshow(train_image, cmap="gray") ax[0, 1].set_title("Training image") ax[1, 0].imshow(gt_image[x_min:x_max, y_min:y_max], cmap="gray") ax[1, 1].imshow(train_image[x_min:x_max, y_min:y_max], cmap="gray")
Shape of the data (100, 1024, 1024) Shape of the training data (1024, 1024) Shape of the ground truth (1024, 1024)
Out[3]:
<matplotlib.image.AxesImage at 0x7fb0d795c370>
Visualize autocorrelation¶
The autocorrelation shows a slight horizontal correlation.
In [4]:
Copied!
# compute autocorrelation
autocorr = autocorrelation(train_image)
# crop the correlation around (0, 0)
midpoint = train_image.shape[0] // 2
crop_size = 10
slices = (
slice(midpoint - crop_size, midpoint + crop_size),
slice(midpoint - crop_size, midpoint + crop_size),
)
# plot image and autocorrelation
fig, ax = plt.subplots(1, 2, figsize=(10, 5))
ax[0].imshow(train_image, cmap="gray")
ax[0].set_title("Training image")
ax[1].imshow(autocorr[slices], cmap="gray")
ax[1].set_title("Autocorrelation")
# compute autocorrelation autocorr = autocorrelation(train_image) # crop the correlation around (0, 0) midpoint = train_image.shape[0] // 2 crop_size = 10 slices = ( slice(midpoint - crop_size, midpoint + crop_size), slice(midpoint - crop_size, midpoint + crop_size), ) # plot image and autocorrelation fig, ax = plt.subplots(1, 2, figsize=(10, 5)) ax[0].imshow(train_image, cmap="gray") ax[0].set_title("Training image") ax[1].imshow(autocorr[slices], cmap="gray") ax[1].set_title("Autocorrelation")
Out[4]:
Text(0.5, 1.0, 'Autocorrelation')
Train with CAREamics¶
The easiest way to use CAREamics is to create a configuration and a CAREamist.
Create configuration¶
The configuration can be built from scratch, giving the user full control over the various parameters available in CAREamics. However, a straightforward way to create a configuration for a particular algorithm is to use one of the convenience functions.
To use structN2V, we simply add the relevant parameters to the configuration.
In [5]:
Copied!
config = create_n2v_configuration(
experiment_name="convallaria_structn2v",
data_type="array",
axes="YX",
patch_size=(64, 64),
batch_size=16,
num_epochs=50,
# structN2V parameters
struct_n2v_axis="horizontal",
struct_n2v_span=11,
# disable augmentations because of the noise correlations
augmentations=[],
)
# only use x-axis flip in augmentations
config.data_config.transforms.insert(
0,
XYFlipModel(flip_x=True, flip_y=False),
)
print(config)
config = create_n2v_configuration( experiment_name="convallaria_structn2v", data_type="array", axes="YX", patch_size=(64, 64), batch_size=16, num_epochs=50, # structN2V parameters struct_n2v_axis="horizontal", struct_n2v_span=11, # disable augmentations because of the noise correlations augmentations=[], ) # only use x-axis flip in augmentations config.data_config.transforms.insert( 0, XYFlipModel(flip_x=True, flip_y=False), ) print(config)
{'algorithm_config': {'algorithm': 'n2v',
'loss': 'n2v',
'lr_scheduler': {'name': 'ReduceLROnPlateau',
'parameters': {}},
'model': {'architecture': 'UNet',
'conv_dims': 2,
'depth': 2,
'final_activation': 'None',
'in_channels': 1,
'independent_channels': True,
'n2v2': False,
'num_channels_init': 32,
'num_classes': 1},
'optimizer': {'name': 'Adam',
'parameters': {'lr': 0.0001}}},
'data_config': {'axes': 'YX',
'batch_size': 16,
'data_type': 'array',
'patch_size': [64, 64],
'transforms': [{'flip_x': True,
'flip_y': False,
'name': 'XYFlip',
'p': 0.5},
{'masked_pixel_percentage': 0.2,
'name': 'N2VManipulate',
'roi_size': 11,
'strategy': 'uniform',
'struct_mask_axis': 'horizontal',
'struct_mask_span': 11}]},
'experiment_name': 'convallaria_structn2v',
'training_config': {'accumulate_grad_batches': 1,
'check_val_every_n_epoch': 1,
'checkpoint_callback': {'auto_insert_metric_name': False,
'mode': 'min',
'monitor': 'val_loss',
'save_last': True,
'save_top_k': 3,
'save_weights_only': False,
'verbose': False},
'enable_progress_bar': True,
'gradient_clip_algorithm': 'norm',
'max_steps': -1,
'num_epochs': 50,
'precision': '32'},
'version': '0.1.0'}
Train¶
A CAREamist can be created using a configuration alone, and then be trained by using the data already loaded in memory.
In [ ]:
Copied!
# instantiate a CAREamist
careamist = CAREamist(source=config)
# train
careamist.train(
train_source=train_image,
)
# instantiate a CAREamist careamist = CAREamist(source=config) # train careamist.train( train_source=train_image, )
Predict with CAREamics¶
Prediction is done with the same CAREamist used for training. Because the image is large we predict using tiling.
In [ ]:
Copied!
prediction = careamist.predict(source=train_image, tile_size=(256, 256))
prediction = careamist.predict(source=train_image, tile_size=(256, 256))
Visualize the prediction¶
In [8]:
Copied!
x_min, x_max = 580, 630
y_min, y_max = 380, 430
psnr_noisy = scale_invariant_psnr(gt_image, train_image)
pnsr_pred = scale_invariant_psnr(gt_image, prediction[0].squeeze())
# plot images
fig, ax = plt.subplots(2, 3, figsize=(8, 5))
ax[0, 0].imshow(train_image, cmap="gray")
ax[0, 0].set_title(f"Training image\npsnr={psnr_noisy:.2f}")
ax[0, 0].set_xticks([])
ax[0, 0].set_yticks([])
ax[0, 1].imshow(prediction[0].squeeze(), cmap="gray")
ax[0, 1].set_title(f"Prediction\npsnr={pnsr_pred:.2f}")
ax[0, 1].set_xticks([])
ax[0, 1].set_yticks([])
ax[0, 2].imshow(gt_image, cmap="gray")
ax[0, 2].set_title("Ground truth")
ax[0, 2].set_xticks([])
ax[0, 2].set_yticks([])
ax[1, 0].imshow(train_image[x_min:x_max, y_min:y_max], cmap="gray")
ax[1, 0].set_xticks([])
ax[1, 0].set_yticks([])
ax[1, 1].imshow(prediction[0].squeeze()[x_min:x_max, y_min:y_max], cmap="gray")
ax[1, 1].set_xticks([])
ax[1, 1].set_yticks([])
ax[1, 2].imshow(gt_image[x_min:x_max, y_min:y_max], cmap="gray")
ax[1, 2].set_xticks([])
ax[1, 2].set_yticks([])
fig.tight_layout()
plt.show()
x_min, x_max = 580, 630 y_min, y_max = 380, 430 psnr_noisy = scale_invariant_psnr(gt_image, train_image) pnsr_pred = scale_invariant_psnr(gt_image, prediction[0].squeeze()) # plot images fig, ax = plt.subplots(2, 3, figsize=(8, 5)) ax[0, 0].imshow(train_image, cmap="gray") ax[0, 0].set_title(f"Training image\npsnr={psnr_noisy:.2f}") ax[0, 0].set_xticks([]) ax[0, 0].set_yticks([]) ax[0, 1].imshow(prediction[0].squeeze(), cmap="gray") ax[0, 1].set_title(f"Prediction\npsnr={pnsr_pred:.2f}") ax[0, 1].set_xticks([]) ax[0, 1].set_yticks([]) ax[0, 2].imshow(gt_image, cmap="gray") ax[0, 2].set_title("Ground truth") ax[0, 2].set_xticks([]) ax[0, 2].set_yticks([]) ax[1, 0].imshow(train_image[x_min:x_max, y_min:y_max], cmap="gray") ax[1, 0].set_xticks([]) ax[1, 0].set_yticks([]) ax[1, 1].imshow(prediction[0].squeeze()[x_min:x_max, y_min:y_max], cmap="gray") ax[1, 1].set_xticks([]) ax[1, 1].set_yticks([]) ax[1, 2].imshow(gt_image[x_min:x_max, y_min:y_max], cmap="gray") ax[1, 2].set_xticks([]) ax[1, 2].set_yticks([]) fig.tight_layout() plt.show()
Compute metrics over the whole stack¶
In [ ]:
Copied!
predictions = careamist.predict(source=data, tile_size=(256, 256), axes="SYX")
predictions = careamist.predict(source=data, tile_size=(256, 256), axes="SYX")
In [10]:
Copied!
psnrs = np.zeros((len(predictions), 1))
for i, image in enumerate(predictions):
psnrs[i] = scale_invariant_psnr(gt_image, image.squeeze())
print(f"PSNR: {psnrs.mean():.2f} +/- {psnrs.std():.2f}")
psnrs = np.zeros((len(predictions), 1)) for i, image in enumerate(predictions): psnrs[i] = scale_invariant_psnr(gt_image, image.squeeze()) print(f"PSNR: {psnrs.mean():.2f} +/- {psnrs.std():.2f}")
PSNR: 31.02 +/- 0.08
Export the model¶
The model is automatically saved during training (the so-called checkpoints) and can be loaded back easily, but you can also export the model to the BioImage Model Zoo format.
In [11]:
Copied!
# create a cover image
x_start, width = 650, 256
y_start, height = 450, 256
# rotate images
rot_train = np.rot90(train_image, 3)
rot_pred = np.rot90(prediction[0].squeeze(), 3)
# create image
cover = np.zeros((height, width))
# normalize train and prediction
norm_train = (rot_train - rot_train.min()) / (rot_train.max() - rot_train.min())
norm_pred = (rot_pred - rot_pred.min()) / (rot_pred.max() - rot_pred.min())
# fill in halves
cover[:, : width // 2] = norm_train[
y_start : y_start + height, x_start : x_start + width // 2
]
cover[:, width // 2 :] = norm_pred[
y_start : y_start + height, x_start + width // 2 : x_start + width
]
# plot the single image
plt.imshow(cover, cmap="gray")
# save the image
im = Image.fromarray(cover * 255)
im = im.convert("L")
im.save("Convallaria_structN2V.jpeg")
# create a cover image x_start, width = 650, 256 y_start, height = 450, 256 # rotate images rot_train = np.rot90(train_image, 3) rot_pred = np.rot90(prediction[0].squeeze(), 3) # create image cover = np.zeros((height, width)) # normalize train and prediction norm_train = (rot_train - rot_train.min()) / (rot_train.max() - rot_train.min()) norm_pred = (rot_pred - rot_pred.min()) / (rot_pred.max() - rot_pred.min()) # fill in halves cover[:, : width // 2] = norm_train[ y_start : y_start + height, x_start : x_start + width // 2 ] cover[:, width // 2 :] = norm_pred[ y_start : y_start + height, x_start + width // 2 : x_start + width ] # plot the single image plt.imshow(cover, cmap="gray") # save the image im = Image.fromarray(cover * 255) im = im.convert("L") im.save("Convallaria_structN2V.jpeg")
In [12]:
Copied!
# Export the model
careamist.export_to_bmz(
path_to_archive="convallaria_structn2v_model.zip",
friendly_model_name="Convallaria_StructN2V",
input_array=train_image.astype(np.float32),
authors=[{"name": "CAREamics authors", "affiliation": "Human Technopole"}],
)
# Export the model careamist.export_to_bmz( path_to_archive="convallaria_structn2v_model.zip", friendly_model_name="Convallaria_StructN2V", input_array=train_image.astype(np.float32), authors=[{"name": "CAREamics authors", "affiliation": "Human Technopole"}], )
LOCAL_RANK: 0 - CUDA_VISIBLE_DEVICES: [0]
Predicting: | | 0/? [00:00<?, ?it/s]
computing SHA256 of inputs.npy (result: f107951e7b7860117ff4f40501aee64acc15f6cdba7bf83876f7be31a7d50274): 100%|██████████| 33/33 [00:00<00:00, 7912.88it/s] computing SHA256 of outputs.npy (result: b8257b57c728d5953888f43ea9d4d9ac1dc50e91e78f4617bab8564b24374d75): 100%|██████████| 33/33 [00:00<00:00, 8033.67it/s] computing SHA256 of environment.yml (result: d743faa28b30f054d6caa699271006e27d3cda48e07c02cf092b4e0732830574): 100%|██████████| 1/1 [00:00<00:00, 1284.23it/s] computing SHA256 of weights.pth (result: d696af5cafd4a6c578ec22cf5643f31c5768cc3bdd81eba8f997ae03297ea80a): 100%|██████████| 16/16 [00:00<00:00, 6931.30it/s] computing SHA256 of config.yml (result: 9fdcc6e487e30e5c29abdbcbaf3d3422a1fe748366bd4350348dcd09fc04da3b): 100%|██████████| 1/1 [00:00<00:00, 1279.14it/s] 2024-09-19 18:27:04.222 | Level 30 | bioimageio.spec._internal.field_warning:issue_warning:149 - documentation: No '# Validation' (sub)section found in /home/melisande.croft/.careamics/README.md. 2024-09-19 18:27:04.748 | INFO | bioimageio.core._resource_tests:_test_model_inference:147 - starting 'Reproduce test outputs from test inputs (pytorch_state_dict)' /localscratch/miniforge3/envs/careamics/lib/python3.10/site-packages/bioimageio/core/model_adapters/_pytorch_model_adapter.py:44: FutureWarning: You are using `torch.load` with `weights_only=False` (the current default value), which uses the default pickle module implicitly. It is possible to construct malicious pickle data which will execute arbitrary code during unpickling (See https://github.com/pytorch/pytorch/blob/main/SECURITY.md#untrusted-models for more details). In a future release, the default value for `weights_only` will be flipped to `True`. This limits the functions that could be executed during unpickling. Arbitrary objects will no longer be allowed to be loaded via this mode unless they are explicitly allowlisted by the user via `torch.serialization.add_safe_globals`. We recommend you start setting `weights_only=True` for any use case where you don't have full control of the loaded file. Please open an issue on GitHub for any issues related to this experimental feature. state: Any = torch.load( 2024-09-19 18:27:05.004 | INFO | bioimageio.core._resource_tests:_test_model_inference_parametrized:242 - Testing inference with 2 different input tensor sizes /localscratch/miniforge3/envs/careamics/lib/python3.10/site-packages/bioimageio/core/model_adapters/_pytorch_model_adapter.py:44: FutureWarning: You are using `torch.load` with `weights_only=False` (the current default value), which uses the default pickle module implicitly. It is possible to construct malicious pickle data which will execute arbitrary code during unpickling (See https://github.com/pytorch/pytorch/blob/main/SECURITY.md#untrusted-models for more details). In a future release, the default value for `weights_only` will be flipped to `True`. This limits the functions that could be executed during unpickling. Arbitrary objects will no longer be allowed to be loaded via this mode unless they are explicitly allowlisted by the user via `torch.serialization.add_safe_globals`. We recommend you start setting `weights_only=True` for any use case where you don't have full control of the loaded file. Please open an issue on GitHub for any issues related to this experimental feature. state: Any = torch.load( computing SHA256 of config.yml (result: 9fdcc6e487e30e5c29abdbcbaf3d3422a1fe748366bd4350348dcd09fc04da3b): 100%|██████████| 1/1 [00:00<00:00, 262.31it/s] computing SHA256 of inputs.npy (result: f107951e7b7860117ff4f40501aee64acc15f6cdba7bf83876f7be31a7d50274): 100%|██████████| 33/33 [00:00<00:00, 1838.85it/s] computing SHA256 of outputs.npy (result: b8257b57c728d5953888f43ea9d4d9ac1dc50e91e78f4617bab8564b24374d75): 100%|██████████| 33/33 [00:00<00:00, 1827.39it/s] computing SHA256 of environment.yml (result: d743faa28b30f054d6caa699271006e27d3cda48e07c02cf092b4e0732830574): 100%|██████████| 1/1 [00:00<00:00, 198.02it/s] computing SHA256 of weights.pth (result: d696af5cafd4a6c578ec22cf5643f31c5768cc3bdd81eba8f997ae03297ea80a): 100%|██████████| 16/16 [00:00<00:00, 1137.55it/s]
In [ ]:
Copied!