diff --git a/examples/demo_dlmbl/convert-solution.py b/examples/demo_dlmbl/convert-solution.py new file mode 100644 index 000000000..279f78743 --- /dev/null +++ b/examples/demo_dlmbl/convert-solution.py @@ -0,0 +1,41 @@ +import argparse +from traitlets.config import Config +import nbformat as nbf +from nbconvert.preprocessors import TagRemovePreprocessor, ClearOutputPreprocessor +from nbconvert.exporters import NotebookExporter + + +def get_arg_parser(): + parser = argparse.ArgumentParser() + + parser.add_argument('input_file') + parser.add_argument('output_file') + + return parser + + +def convert(input_file, output_file): + c = Config() + c.TagRemovePreprocessor.remove_cell_tags = ("solution",) + c.TagRemovePreprocessor.enabled = True + c.ClearOutputPreprocesser.enabled = True + c.NotebookExporter.preprocessors = [ + "nbconvert.preprocessors.TagRemovePreprocessor", + "nbconvert.preprocessors.ClearOutputPreprocessor" + ] + + exporter = NotebookExporter(config=c) + exporter.register_preprocessor(TagRemovePreprocessor(config=c), True) + exporter.register_preprocessor(ClearOutputPreprocessor(), True) + + output = NotebookExporter(config=c).from_filename(input_file) + with open(output_file, 'w') as f: + f.write(output[0]) + + +if __name__ == "__main__": + parser = get_arg_parser() + args = parser.parse_args() + + convert(args.input_file, args.output_file) + print(f'Converted {args.input_file} to {args.output_file}') diff --git a/examples/demo_dlmbl/exercise.ipynb b/examples/demo_dlmbl/exercise.ipynb new file mode 100644 index 000000000..b6a3a256e --- /dev/null +++ b/examples/demo_dlmbl/exercise.ipynb @@ -0,0 +1,754 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "c9438eb5", + "metadata": { + "cell_marker": "\"\"\"" + }, + "source": [ + "# Image translation\n", + "---\n", + "\n", + "Written by Ziwen Liu and Shalin Mehta, CZ Biohub San Francisco.\n", + "---\n", + "\n", + "In this exercise, we will solve an image translation task to predict fluorescence images of nuclei and membrane markers from quantitative phase images of cells. In other words, we will _virtually stain_ the nuclei and membrane visible in the phase image. \n", + "\n", + "Here, the source domain is label-free microscopy (material density) and the target domain is fluorescence microscopy (fluorophore density). The goal is to learn a mapping from the source domain to the target domain. We will use a deep convolutional neural network (CNN), specifically, a U-Net model with residual connections to learn the mapping. The preprocessing, training, prediction, evaluation, and deployment steps are unified in a computer vision pipeline for single-cell analysis that we call [VisCy](https://github.com/mehta-lab/VisCy).\n", + "\n", + "VisCy evolved from our previous work on virtual staining of cellular components from their density and anisotropy.\n", + "![](https://iiif.elifesciences.org/lax/55502%2Felife-55502-fig1-v2.tif/full/1500,/0/default.jpg)\n", + "\n", + "[Guo et al. (2020) Revealing architectural order with quantitative label-free imaging and deep learning\n", + ". eLife](https://elifesciences.org/articles/55502).\n", + "\n", + "VisCy exploits recent advances in the data and metadata formats ([OME-zarr](https://www.nature.com/articles/s41592-021-01326-w)) and DL frameworks, [PyTorch Lightning](https://lightning.ai/) and [MONAI](https://monai.io/). Our previous pipeline, [microDL](https://github.com/mehta-lab/microDL), is deprecated and is now a public archive." + ] + }, + { + "cell_type": "markdown", + "id": "b36463af", + "metadata": { + "cell_marker": "\"\"\"", + "lines_to_next_cell": 0 + }, + "source": [ + "Today, we will train a 2D image translation model using a 2D U-Net with residual connections. We will use a dataset of 301 fields of view (FOVs) of Human Embryonic Kidney (HEK) cells, each FOV has 3 channels (phase, membrane, and nuclei). The cells were labeled with CRISPR editing. Intrestingly, not all cells during this experiment were labeled due to the stochastic nature of CRISPR editing. In such situations, virtual staining rescues missing labels.\n", + "![HEK](https://github.com/mehta-lab/VisCy/blob/dlmbl2023/docs/figures/phase_to_nuclei_membrane.svg?raw=true)\n", + "\n", + "
\n", + "The exercise is organized in 3 parts.\n", + "\n", + "* **Part 1** - Explore the data using tensorboard. Launch the training before lunch.\n", + "* Lunch break - The model will continue training during lunch.\n", + "* **Part 2** - Evaluate the training with tensorboard. Train another model.\n", + "* **Part 3** - Tune the models to improve performance.\n", + "
\n", + "\n", + "📖 As you work through parts 2 and 3, please share the layouts of the models you train and their performance with everyone via [this google doc](https://docs.google.com/document/d/1hZWSVRvt9KJEdYu7ib-vFBqAVQRYL8cWaP_vFznu7D8/edit#heading=h.n5u485pmzv2z) 📖.\n", + "\n", + "\n", + "Our guesstimate is that each of the three parts will take ~1.5 hours, but don't rush parts 1 and 2 if you need more time with them.\n", + "We will discuss your observations on google doc after checkpoints 2 and 3. The exercise is focused on understanding information contained in data, process of training and evaluating image translation models, and parameter exploration.\n", + "There are a few coding tasks sprinkled in.\n", + "\n", + "\n", + "Before you start,\n", + "\n", + "
\n", + "Set your python kernel to 04-image-translation\n", + "
" + ] + }, + { + "cell_type": "markdown", + "id": "ff42876f", + "metadata": { + "cell_marker": "\"\"\"", + "incorrectly_encoded_metadata": "id='1_phase2fluor'>", + "title": "\n", + "Task 1.1\n", + "Use \n", + "iohub.open_ome_zarr to read the dataset and explore several FOVs with matplotlib.\n", + "\n", + "\n", + "There should be 301 FOVs in the dataset (12 GB compressed).\n", + "\n", + "Each FOV consists of 3 channels of 2048x2048 images,\n", + "saved in the \n", + "High-Content Screening (HCS) layout\n", + "specified by the Open Microscopy Environment Next Generation File Format\n", + "(OME-NGFF).\n", + "\n", + "The layout on the disk is: row/col/field/pyramid_level/timepoint/channel/z/y/x.\n", + "Notice that labelling of nuclei channel is not complete - some cells are not expressing the fluorescent protein." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e0241037", + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "dataset = open_ome_zarr(data_path)\n", + "\n", + "print(f\"Number of positions: {len(list(dataset.positions()))}\")\n", + "\n", + "# Use the field and pyramid_level below to visualize data.\n", + "row = \"0\"\n", + "col = \"0\"\n", + "field = \"23\"\n", + "\n", + "# This dataset contains images at 3 resolutions.\n", + "# '0' is the highest resolution\n", + "# '1' is down-scaled 2x2,\n", + "# '2' is down-scaled 4x4.\n", + "# Such datasets are called image pyramids.\n", + "pyaramid_level = \"0\"\n", + "\n", + "# `channel_names` is the metadata that is stored with data according to the OME-NGFF spec.\n", + "n_channels = len(dataset.channel_names)\n", + "\n", + "image = dataset[f\"{row}/{col}/{field}/{pyaramid_level}\"].numpy()\n", + "print(f\"data shape: {image.shape}, FOV: {field}, pyramid level: {pyaramid_level}\")\n", + "\n", + "figure, axes = plt.subplots(1, n_channels, figsize=(9, 3))\n", + "\n", + "for i in range(n_channels):\n", + " for i in range(n_channels):\n", + " channel_image = image[0, i, 0]\n", + " # Adjust contrast to 0.5th and 99.5th percentile of pixel values.\n", + " p_low, p_high = np.percentile(channel_image, (0.5, 99.5))\n", + " channel_image = np.clip(channel_image, p_low, p_high)\n", + " axes[i].imshow(channel_image, cmap=\"gray\")\n", + " axes[i].axis(\"off\")\n", + " axes[i].set_title(dataset.channel_names[i])\n", + "plt.tight_layout()" + ] + }, + { + "cell_type": "markdown", + "id": "89e676bf", + "metadata": { + "cell_marker": "\"\"\"" + }, + "source": [ + "## Initialize data loaders and see the samples in tensorboard.\n", + "\n", + "
\n", + "Task 1.2\n", + "Setup the data loader and log several batches to tensorboard.\n", + "
`\n", + "\n", + "VisCy builds on top of PyTorch Lightning. PyTorch Lightning is a thin wrapper around PyTorch that allows rapid experimentation. It provides a [DataModule](https://lightning.ai/docs/pytorch/stable/data/datamodule.html) to handle loading and processing of data during training. VisCy provides a child class, `HCSDataModule` to make it intuitve to access data stored in the HCS layout.\n", + " \n", + "The dataloader in `HCSDataModule` returns a batch of samples. A `batch` is a list of dictionaries. The length of the list is equal to the batch size. Each dictionary consists of following key-value pairs.\n", + "- `source`: the input image, a tensor of size 1*1*Y*X\n", + "- `target`: the target image, a tensor of size 2*1*Y*X\n", + "- `index` : the tuple of (location of field in HCS layout, time, and z-slice) of the sample." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "05ce4dbb", + "metadata": {}, + "outputs": [], + "source": [ + "# Define a function to write a batch to tensorboard log.\n", + "\n", + "\n", + "def log_batch_tensorboard(batch, batchno, writer, card_name):\n", + " \"\"\"\n", + " Logs a batch of images to TensorBoard.\n", + "\n", + " Args:\n", + " batch (dict): A dictionary containing the batch of images to be logged.\n", + " writer (SummaryWriter): A TensorBoard SummaryWriter object.\n", + " card_name (str): The name of the card to be displayed in TensorBoard.\n", + "\n", + " Returns:\n", + " None\n", + " \"\"\"\n", + " batch_phase = batch[\"source\"][:, :, 0, :, :] # batch_size x z_size x Y x X tensor.\n", + " batch_membrane = batch[\"target\"][:, 1, 0, :, :].unsqueeze(\n", + " 1\n", + " ) # batch_size x 1 x Y x X tensor.\n", + " batch_nuclei = batch[\"target\"][:, 0, 0, :, :].unsqueeze(\n", + " 1\n", + " ) # batch_size x 1 x Y x X tensor.\n", + "\n", + " p1, p99 = np.percentile(batch_membrane, (0.1, 99.9))\n", + " batch_membrane = np.clip((batch_membrane - p1) / (p99 - p1), 0, 1)\n", + "\n", + " p1, p99 = np.percentile(batch_nuclei, (0.1, 99.9))\n", + " batch_nuclei = np.clip((batch_nuclei - p1) / (p99 - p1), 0, 1)\n", + "\n", + " p1, p99 = np.percentile(batch_phase, (0.1, 99.9))\n", + " batch_phase = np.clip((batch_phase - p1) / (p99 - p1), 0, 1)\n", + "\n", + " [N, C, H, W] = batch_phase.shape\n", + " interleaved_images = torch.zeros((3 * N, C, H, W), dtype=batch_phase.dtype)\n", + " interleaved_images[0::3, :] = batch_phase\n", + " interleaved_images[1::3, :] = batch_nuclei\n", + " interleaved_images[2::3, :] = batch_membrane\n", + "\n", + " grid = torchvision.utils.make_grid(interleaved_images, nrow=3)\n", + "\n", + " # add the grid to tensorboard\n", + " writer.add_image(card_name, grid, batchno)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f627e8e8", + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "# Initialize the data module.\n", + "\n", + "BATCH_SIZE = 42\n", + "# 42 is a perfectly reasonable batch size. After all, it is the answer to the ultimate question of life, the universe and everything.\n", + "# More seriously, batch size does not have to be a power of 2.\n", + "# See: https://sebastianraschka.com/blog/2022/batch-size-2.html\n", + "\n", + "data_module = HCSDataModule(\n", + " data_path,\n", + " source_channel=\"Phase\",\n", + " target_channel=[\"Nuclei\", \"Membrane\"],\n", + " z_window_size=1,\n", + " split_ratio=0.8,\n", + " batch_size=BATCH_SIZE,\n", + " num_workers=8,\n", + " architecture=\"2D\",\n", + " yx_patch_size=(512, 512), # larger patch size makes it easy to see augmentations.\n", + " augment=False, # Turn off augmentation for now.\n", + ")\n", + "data_module.setup(\"fit\")\n", + "\n", + "print(\n", + " f\"FOVs in training set: {len(data_module.train_dataset)}, FOVs in validation set:{len(data_module.val_dataset)}\"\n", + ")\n", + "train_dataloader = data_module.train_dataloader()\n", + "\n", + "# Instantiate the tensorboard SummaryWriter, logs the first batch and then iterates through all the batches and logs them to tensorboard.\n", + "\n", + "writer = SummaryWriter(log_dir=f\"{log_dir}/view_batch\")\n", + "# Draw a batch and write to tensorboard.\n", + "batch = next(iter(train_dataloader))\n", + "log_batch_tensorboard(batch, 0, writer, \"augmentation/none\")\n", + "\n", + "# Iterate through all the batches and log them to tensorboard.\n", + "for i, batch in enumerate(train_dataloader):\n", + " log_batch_tensorboard(batch, i, writer, \"augmentation/none\")\n", + "writer.close()" + ] + }, + { + "cell_type": "markdown", + "id": "f6ffc7d1", + "metadata": { + "cell_marker": "\"\"\"" + }, + "source": [ + "There are multiple ways of seeing the tensorboard.\n", + "1. Jupyter lab forwards the tensorboard port to the browser. Go to http://localhost:6006/ to see the tensorboard.\n", + "2. You likely have an open viewer in the first cell where you loaded tensorboard jupyter extension.\n", + "3. If you want to see tensorboard in a specific cell, use the following code.\n", + "```\n", + "notebook.list() # View open TensorBoard instances\n", + "notebook.display(port=6006, height=800) # Display the TensorBoard instance specified by the port.\n", + "```" + ] + }, + { + "cell_type": "markdown", + "id": "ed1448ff", + "metadata": { + "cell_marker": "\"\"\"", + "lines_to_next_cell": 0 + }, + "source": [ + "## View augmentations using tensorboard.\n", + "\n", + "
\n", + "Task 1.3\n", + "Turn on augmentation and view the batch in tensorboard." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "245ae4dd", + "metadata": { + "lines_to_next_cell": 2 + }, + "outputs": [], + "source": [ + "##########################\n", + "######## TODO ########\n", + "##########################\n", + "\n", + "# Write code to turn on augmentations, change batch sizes and log them to tensorboard.\n", + "# See how the training data changes as a function of these parameters.\n", + "# Remember to call `data_module.setup(\"fit\")` after changing the parameters." + ] + }, + { + "cell_type": "markdown", + "id": "85743be9", + "metadata": { + "cell_marker": "\"\"\"" + }, + "source": [ + "## Construct a 2D U-Net for image translation.\n", + "See ``viscy.unet.networks.Unet2D.Unet2d`` for configuration details.\n", + "We setup a fresh data module and instantiate the trainer class." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "44fddf02", + "metadata": { + "lines_to_next_cell": 2 + }, + "outputs": [], + "source": [ + "\n", + "# The entire training loop is contained in this cell.\n", + "\n", + "GPU_ID = 0\n", + "BATCH_SIZE = 10\n", + "YX_PATCH_SIZE = (512, 512)\n", + "\n", + "\n", + "# Dictionary that specifies key parameters of the model.\n", + "phase2fluor_config = {\n", + " \"architecture\": \"2D\",\n", + " \"num_filters\": [24, 48, 96, 192, 384],\n", + " \"in_channels\": 1,\n", + " \"out_channels\": 2,\n", + " \"residual\": True,\n", + " \"dropout\": 0.1, # dropout randomly turns off weights to avoid overfitting of the model to data.\n", + " \"task\": \"reg\", # reg = regression task.\n", + "}\n", + "\n", + "phase2fluor_model = VSUNet(\n", + " model_config=phase2fluor_config.copy(),\n", + " batch_size=BATCH_SIZE,\n", + " loss_function=torch.nn.functional.l1_loss,\n", + " schedule=\"WarmupCosine\",\n", + " log_num_samples=10, # Number of samples from each batch to log to tensorboard.\n", + " example_input_yx_shape=YX_PATCH_SIZE,\n", + ")\n", + "\n", + "# Reinitialize the data module.\n", + "phase2fluor_data = HCSDataModule(\n", + " data_path,\n", + " source_channel=\"Phase\",\n", + " target_channel=[\"Nuclei\", \"Membrane\"],\n", + " z_window_size=1,\n", + " split_ratio=0.8,\n", + " batch_size=BATCH_SIZE,\n", + " num_workers=8,\n", + " architecture=\"2D\",\n", + " yx_patch_size=YX_PATCH_SIZE,\n", + " augment=True,\n", + ")\n", + "phase2fluor_data.setup(\"fit\")" + ] + }, + { + "cell_type": "markdown", + "id": "292cb22a", + "metadata": { + "cell_marker": "\"\"\"" + }, + "source": [ + "
\n", + "Task 1.4\n", + "Setup the training for ~30 epochs\n", + "
\n", + "\n", + "Tips:\n", + "- Set ``default_root_dir`` to store the logs and checkpoints\n", + "in a specific directory." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b9cf215b", + "metadata": { + "lines_to_next_cell": 2, + "title": "Setup trainer and check for errors." + }, + "outputs": [], + "source": [ + "\n", + "# fast_dev_run runs a single batch of data through the model to check for errors.\n", + "trainer = VSTrainer(accelerator=\"gpu\", devices=[GPU_ID], fast_dev_run=True)\n", + "\n", + "# trainer class takes the model and the data module as inputs.\n", + "trainer.fit(phase2fluor_model, datamodule=phase2fluor_data)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5f11d604", + "metadata": { + "lines_to_next_cell": 2 + }, + "outputs": [], + "source": [ + "\n", + "GPU_ID = 0\n", + "n_samples = len(phase2fluor_data.train_dataset)\n", + "steps_per_epoch = n_samples // BATCH_SIZE # steps per epoch.\n", + "n_epochs = 30\n", + "\n", + "trainer = VSTrainer(\n", + " accelerator=\"gpu\",\n", + " devices=[GPU_ID],\n", + " max_epochs=n_epochs,\n", + " log_every_n_steps=steps_per_epoch // 2,\n", + " # log losses and image samples 2 times per epoch.\n", + " default_root_dir=Path(\n", + " log_dir, \"phase2fluor\"\n", + " ), # lightning trainer transparently saves logs and model checkpoints in this directory.\n", + ")\n", + "\n", + "# Log graph\n", + "trainer.logger.log_graph(phase2fluor_model, phase2fluor_data.train_dataset[0][\"source\"])\n", + "# Launch training.\n", + "trainer.fit(phase2fluor_model, datamodule=phase2fluor_data)" + ] + }, + { + "cell_type": "markdown", + "id": "157a989c", + "metadata": { + "cell_marker": "\"\"\"" + }, + "source": [ + "
\n", + "Checkpoint 1\n", + "\n", + "Now the training has started,\n", + "we can come back after a while and evaluate the performance!\n", + "
" + ] + }, + { + "cell_type": "markdown", + "id": "d68dab93", + "metadata": { + "cell_marker": "\"\"\"", + "incorrectly_encoded_metadata": "id='1_fluor2phase'>", + "title": "\n", + "Checkpoint 2\n", + "Please summarize hyperparameters and performance of your models in [this google doc](https://docs.google.com/document/d/1hZWSVRvt9KJEdYu7ib-vFBqAVQRYL8cWaP_vFznu7D8/edit#heading=h.n5u485pmzv2z)\n", + "\n", + "Now that you have trained two models, let's think about the following questions:\n", + "- What is the information content of each channel in the dataset?\n", + "- How would you use image translation models?\n", + "- What can you try to improve the performance of each model?\n", + "\n", + "\n", + "
" + ] + }, + { + "cell_type": "markdown", + "id": "7d0c4204", + "metadata": { + "cell_marker": "\"\"\"", + "incorrectly_encoded_metadata": "id='3_tuning'>", + "title": "\n", + "Checkpoint 3\n", + "\n", + "Congratulations! You have trained several image translation models now!\n", + "Please document hyperparameters, snapshots of predictioons on validation set, and loss curves for your models in [this google doc](https://docs.google.com/document/d/1hZWSVRvt9KJEdYu7ib-vFBqAVQRYL8cWaP_vFznu7D8/edit#heading=h.n5u485pmzv2z)\n", + "" + ] + } + ], + "metadata": { + "jupytext": { + "cell_metadata_filter": "all", + "main_language": "python" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/examples/demo_dlmbl/prepare-exercise.sh b/examples/demo_dlmbl/prepare-exercise.sh new file mode 100644 index 000000000..88f2ba1b1 --- /dev/null +++ b/examples/demo_dlmbl/prepare-exercise.sh @@ -0,0 +1,10 @@ +# Run black on .py files +# black solution.py + +# Convert .py to ipynb +# "cell_metadata_filter": "all" preserve cell tags including our solution tags +jupytext --to ipynb --update-metadata '{"jupytext": {"cell_metadata_filter":"all"}}' solution.py + +# Create the exercise notebook by removing cell outputs and deleting cells tagged with "solution" +# There is a bug in the nbconvert cli so we need to use the python API instead +python convert-solution.py solution.ipynb exercise.ipynb diff --git a/examples/demo_dlmbl/solution.ipynb b/examples/demo_dlmbl/solution.ipynb new file mode 100644 index 000000000..c97b85779 --- /dev/null +++ b/examples/demo_dlmbl/solution.ipynb @@ -0,0 +1,929 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "c9438eb5", + "metadata": { + "cell_marker": "\"\"\"" + }, + "source": [ + "# Image translation\n", + "---\n", + "\n", + "Written by Ziwen Liu and Shalin Mehta, CZ Biohub San Francisco.\n", + "---\n", + "\n", + "In this exercise, we will solve an image translation task to predict fluorescence images of nuclei and membrane markers from quantitative phase images of cells. In other words, we will _virtually stain_ the nuclei and membrane visible in the phase image. \n", + "\n", + "Here, the source domain is label-free microscopy (material density) and the target domain is fluorescence microscopy (fluorophore density). The goal is to learn a mapping from the source domain to the target domain. We will use a deep convolutional neural network (CNN), specifically, a U-Net model with residual connections to learn the mapping. The preprocessing, training, prediction, evaluation, and deployment steps are unified in a computer vision pipeline for single-cell analysis that we call [VisCy](https://github.com/mehta-lab/VisCy).\n", + "\n", + "VisCy evolved from our previous work on virtual staining of cellular components from their density and anisotropy.\n", + "![](https://iiif.elifesciences.org/lax/55502%2Felife-55502-fig1-v2.tif/full/1500,/0/default.jpg)\n", + "\n", + "[Guo et al. (2020) Revealing architectural order with quantitative label-free imaging and deep learning\n", + ". eLife](https://elifesciences.org/articles/55502).\n", + "\n", + "VisCy exploits recent advances in the data and metadata formats ([OME-zarr](https://www.nature.com/articles/s41592-021-01326-w)) and DL frameworks, [PyTorch Lightning](https://lightning.ai/) and [MONAI](https://monai.io/). Our previous pipeline, [microDL](https://github.com/mehta-lab/microDL), is deprecated and is now a public archive." + ] + }, + { + "cell_type": "markdown", + "id": "b36463af", + "metadata": { + "cell_marker": "\"\"\"", + "lines_to_next_cell": 0 + }, + "source": [ + "Today, we will train a 2D image translation model using a 2D U-Net with residual connections. We will use a dataset of 301 fields of view (FOVs) of Human Embryonic Kidney (HEK) cells, each FOV has 3 channels (phase, membrane, and nuclei). The cells were labeled with CRISPR editing. Intrestingly, not all cells during this experiment were labeled due to the stochastic nature of CRISPR editing. In such situations, virtual staining rescues missing labels.\n", + "![HEK](https://github.com/mehta-lab/VisCy/blob/dlmbl2023/docs/figures/phase_to_nuclei_membrane.svg?raw=true)\n", + "\n", + "
\n", + "The exercise is organized in 3 parts.\n", + "\n", + "* **Part 1** - Explore the data using tensorboard. Launch the training before lunch.\n", + "* Lunch break - The model will continue training during lunch.\n", + "* **Part 2** - Evaluate the training with tensorboard. Train another model.\n", + "* **Part 3** - Tune the models to improve performance.\n", + "
\n", + "\n", + "📖 As you work through parts 2 and 3, please share the layouts of the models you train and their performance with everyone via [this google doc](https://docs.google.com/document/d/1hZWSVRvt9KJEdYu7ib-vFBqAVQRYL8cWaP_vFznu7D8/edit#heading=h.n5u485pmzv2z) 📖.\n", + "\n", + "\n", + "Our guesstimate is that each of the three parts will take ~1.5 hours, but don't rush parts 1 and 2 if you need more time with them.\n", + "We will discuss your observations on google doc after checkpoints 2 and 3. The exercise is focused on understanding information contained in data, process of training and evaluating image translation models, and parameter exploration.\n", + "There are a few coding tasks sprinkled in.\n", + "\n", + "\n", + "Before you start,\n", + "\n", + "
\n", + "Set your python kernel to 04-image-translation\n", + "
" + ] + }, + { + "cell_type": "markdown", + "id": "ff42876f", + "metadata": { + "cell_marker": "\"\"\"", + "incorrectly_encoded_metadata": "id='1_phase2fluor'>", + "title": "\n", + "Task 1.1\n", + "Use \n", + "iohub.open_ome_zarr to read the dataset and explore several FOVs with matplotlib.\n", + "\n", + "\n", + "There should be 301 FOVs in the dataset (12 GB compressed).\n", + "\n", + "Each FOV consists of 3 channels of 2048x2048 images,\n", + "saved in the \n", + "High-Content Screening (HCS) layout\n", + "specified by the Open Microscopy Environment Next Generation File Format\n", + "(OME-NGFF).\n", + "\n", + "The layout on the disk is: row/col/field/pyramid_level/timepoint/channel/z/y/x.\n", + "Notice that labelling of nuclei channel is not complete - some cells are not expressing the fluorescent protein." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e0241037", + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "dataset = open_ome_zarr(data_path)\n", + "\n", + "print(f\"Number of positions: {len(list(dataset.positions()))}\")\n", + "\n", + "# Use the field and pyramid_level below to visualize data.\n", + "row = \"0\"\n", + "col = \"0\"\n", + "field = \"23\"\n", + "\n", + "# This dataset contains images at 3 resolutions.\n", + "# '0' is the highest resolution\n", + "# '1' is down-scaled 2x2,\n", + "# '2' is down-scaled 4x4.\n", + "# Such datasets are called image pyramids.\n", + "pyaramid_level = \"0\"\n", + "\n", + "# `channel_names` is the metadata that is stored with data according to the OME-NGFF spec.\n", + "n_channels = len(dataset.channel_names)\n", + "\n", + "image = dataset[f\"{row}/{col}/{field}/{pyaramid_level}\"].numpy()\n", + "print(f\"data shape: {image.shape}, FOV: {field}, pyramid level: {pyaramid_level}\")\n", + "\n", + "figure, axes = plt.subplots(1, n_channels, figsize=(9, 3))\n", + "\n", + "for i in range(n_channels):\n", + " for i in range(n_channels):\n", + " channel_image = image[0, i, 0]\n", + " # Adjust contrast to 0.5th and 99.5th percentile of pixel values.\n", + " p_low, p_high = np.percentile(channel_image, (0.5, 99.5))\n", + " channel_image = np.clip(channel_image, p_low, p_high)\n", + " axes[i].imshow(channel_image, cmap=\"gray\")\n", + " axes[i].axis(\"off\")\n", + " axes[i].set_title(dataset.channel_names[i])\n", + "plt.tight_layout()" + ] + }, + { + "cell_type": "markdown", + "id": "89e676bf", + "metadata": { + "cell_marker": "\"\"\"" + }, + "source": [ + "## Initialize data loaders and see the samples in tensorboard.\n", + "\n", + "
\n", + "Task 1.2\n", + "Setup the data loader and log several batches to tensorboard.\n", + "
`\n", + "\n", + "VisCy builds on top of PyTorch Lightning. PyTorch Lightning is a thin wrapper around PyTorch that allows rapid experimentation. It provides a [DataModule](https://lightning.ai/docs/pytorch/stable/data/datamodule.html) to handle loading and processing of data during training. VisCy provides a child class, `HCSDataModule` to make it intuitve to access data stored in the HCS layout.\n", + " \n", + "The dataloader in `HCSDataModule` returns a batch of samples. A `batch` is a list of dictionaries. The length of the list is equal to the batch size. Each dictionary consists of following key-value pairs.\n", + "- `source`: the input image, a tensor of size 1*1*Y*X\n", + "- `target`: the target image, a tensor of size 2*1*Y*X\n", + "- `index` : the tuple of (location of field in HCS layout, time, and z-slice) of the sample." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "05ce4dbb", + "metadata": {}, + "outputs": [], + "source": [ + "# Define a function to write a batch to tensorboard log.\n", + "\n", + "\n", + "def log_batch_tensorboard(batch, batchno, writer, card_name):\n", + " \"\"\"\n", + " Logs a batch of images to TensorBoard.\n", + "\n", + " Args:\n", + " batch (dict): A dictionary containing the batch of images to be logged.\n", + " writer (SummaryWriter): A TensorBoard SummaryWriter object.\n", + " card_name (str): The name of the card to be displayed in TensorBoard.\n", + "\n", + " Returns:\n", + " None\n", + " \"\"\"\n", + " batch_phase = batch[\"source\"][:, :, 0, :, :] # batch_size x z_size x Y x X tensor.\n", + " batch_membrane = batch[\"target\"][:, 1, 0, :, :].unsqueeze(\n", + " 1\n", + " ) # batch_size x 1 x Y x X tensor.\n", + " batch_nuclei = batch[\"target\"][:, 0, 0, :, :].unsqueeze(\n", + " 1\n", + " ) # batch_size x 1 x Y x X tensor.\n", + "\n", + " p1, p99 = np.percentile(batch_membrane, (0.1, 99.9))\n", + " batch_membrane = np.clip((batch_membrane - p1) / (p99 - p1), 0, 1)\n", + "\n", + " p1, p99 = np.percentile(batch_nuclei, (0.1, 99.9))\n", + " batch_nuclei = np.clip((batch_nuclei - p1) / (p99 - p1), 0, 1)\n", + "\n", + " p1, p99 = np.percentile(batch_phase, (0.1, 99.9))\n", + " batch_phase = np.clip((batch_phase - p1) / (p99 - p1), 0, 1)\n", + "\n", + " [N, C, H, W] = batch_phase.shape\n", + " interleaved_images = torch.zeros((3 * N, C, H, W), dtype=batch_phase.dtype)\n", + " interleaved_images[0::3, :] = batch_phase\n", + " interleaved_images[1::3, :] = batch_nuclei\n", + " interleaved_images[2::3, :] = batch_membrane\n", + "\n", + " grid = torchvision.utils.make_grid(interleaved_images, nrow=3)\n", + "\n", + " # add the grid to tensorboard\n", + " writer.add_image(card_name, grid, batchno)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f627e8e8", + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "# Initialize the data module.\n", + "\n", + "BATCH_SIZE = 42\n", + "# 42 is a perfectly reasonable batch size. After all, it is the answer to the ultimate question of life, the universe and everything.\n", + "# More seriously, batch size does not have to be a power of 2.\n", + "# See: https://sebastianraschka.com/blog/2022/batch-size-2.html\n", + "\n", + "data_module = HCSDataModule(\n", + " data_path,\n", + " source_channel=\"Phase\",\n", + " target_channel=[\"Nuclei\", \"Membrane\"],\n", + " z_window_size=1,\n", + " split_ratio=0.8,\n", + " batch_size=BATCH_SIZE,\n", + " num_workers=8,\n", + " architecture=\"2D\",\n", + " yx_patch_size=(512, 512), # larger patch size makes it easy to see augmentations.\n", + " augment=False, # Turn off augmentation for now.\n", + ")\n", + "data_module.setup(\"fit\")\n", + "\n", + "print(\n", + " f\"FOVs in training set: {len(data_module.train_dataset)}, FOVs in validation set:{len(data_module.val_dataset)}\"\n", + ")\n", + "train_dataloader = data_module.train_dataloader()\n", + "\n", + "# Instantiate the tensorboard SummaryWriter, logs the first batch and then iterates through all the batches and logs them to tensorboard.\n", + "\n", + "writer = SummaryWriter(log_dir=f\"{log_dir}/view_batch\")\n", + "# Draw a batch and write to tensorboard.\n", + "batch = next(iter(train_dataloader))\n", + "log_batch_tensorboard(batch, 0, writer, \"augmentation/none\")\n", + "\n", + "# Iterate through all the batches and log them to tensorboard.\n", + "for i, batch in enumerate(train_dataloader):\n", + " log_batch_tensorboard(batch, i, writer, \"augmentation/none\")\n", + "writer.close()" + ] + }, + { + "cell_type": "markdown", + "id": "f6ffc7d1", + "metadata": { + "cell_marker": "\"\"\"" + }, + "source": [ + "There are multiple ways of seeing the tensorboard.\n", + "1. Jupyter lab forwards the tensorboard port to the browser. Go to http://localhost:6006/ to see the tensorboard.\n", + "2. You likely have an open viewer in the first cell where you loaded tensorboard jupyter extension.\n", + "3. If you want to see tensorboard in a specific cell, use the following code.\n", + "```\n", + "notebook.list() # View open TensorBoard instances\n", + "notebook.display(port=6006, height=800) # Display the TensorBoard instance specified by the port.\n", + "```" + ] + }, + { + "cell_type": "markdown", + "id": "ed1448ff", + "metadata": { + "cell_marker": "\"\"\"", + "lines_to_next_cell": 0 + }, + "source": [ + "## View augmentations using tensorboard.\n", + "\n", + "
\n", + "Task 1.3\n", + "Turn on augmentation and view the batch in tensorboard." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "245ae4dd", + "metadata": { + "lines_to_next_cell": 2 + }, + "outputs": [], + "source": [ + "##########################\n", + "######## TODO ########\n", + "##########################\n", + "\n", + "# Write code to turn on augmentations, change batch sizes and log them to tensorboard.\n", + "# See how the training data changes as a function of these parameters.\n", + "# Remember to call `data_module.setup(\"fit\")` after changing the parameters." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e14f4faf", + "metadata": { + "tags": [ + "solution" + ] + }, + "outputs": [], + "source": [ + "##########################\n", + "######## Solution ########\n", + "##########################\n", + "\n", + "data_module.augment = True\n", + "data_module.batch_size = 21\n", + "data_module.split_ratio = 0.8\n", + "data_module.setup(\"fit\")\n", + "\n", + "train_dataloader = data_module.train_dataloader()\n", + "# Draw batches and write to tensorboard\n", + "writer = SummaryWriter(log_dir=f\"{log_dir}/view_batch\")\n", + "for i, batch in enumerate(train_dataloader):\n", + " log_batch_tensorboard(batch, i, writer, \"augmentation/some\")\n", + "writer.close()" + ] + }, + { + "cell_type": "markdown", + "id": "85743be9", + "metadata": { + "cell_marker": "\"\"\"" + }, + "source": [ + "## Construct a 2D U-Net for image translation.\n", + "See ``viscy.unet.networks.Unet2D.Unet2d`` for configuration details.\n", + "We setup a fresh data module and instantiate the trainer class." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "44fddf02", + "metadata": { + "lines_to_next_cell": 2 + }, + "outputs": [], + "source": [ + "\n", + "# The entire training loop is contained in this cell.\n", + "\n", + "GPU_ID = 0\n", + "BATCH_SIZE = 10\n", + "YX_PATCH_SIZE = (512, 512)\n", + "\n", + "\n", + "# Dictionary that specifies key parameters of the model.\n", + "phase2fluor_config = {\n", + " \"architecture\": \"2D\",\n", + " \"num_filters\": [24, 48, 96, 192, 384],\n", + " \"in_channels\": 1,\n", + " \"out_channels\": 2,\n", + " \"residual\": True,\n", + " \"dropout\": 0.1, # dropout randomly turns off weights to avoid overfitting of the model to data.\n", + " \"task\": \"reg\", # reg = regression task.\n", + "}\n", + "\n", + "phase2fluor_model = VSUNet(\n", + " model_config=phase2fluor_config.copy(),\n", + " batch_size=BATCH_SIZE,\n", + " loss_function=torch.nn.functional.l1_loss,\n", + " schedule=\"WarmupCosine\",\n", + " log_num_samples=10, # Number of samples from each batch to log to tensorboard.\n", + " example_input_yx_shape=YX_PATCH_SIZE,\n", + ")\n", + "\n", + "# Reinitialize the data module.\n", + "phase2fluor_data = HCSDataModule(\n", + " data_path,\n", + " source_channel=\"Phase\",\n", + " target_channel=[\"Nuclei\", \"Membrane\"],\n", + " z_window_size=1,\n", + " split_ratio=0.8,\n", + " batch_size=BATCH_SIZE,\n", + " num_workers=8,\n", + " architecture=\"2D\",\n", + " yx_patch_size=YX_PATCH_SIZE,\n", + " augment=True,\n", + ")\n", + "phase2fluor_data.setup(\"fit\")" + ] + }, + { + "cell_type": "markdown", + "id": "292cb22a", + "metadata": { + "cell_marker": "\"\"\"" + }, + "source": [ + "
\n", + "Task 1.4\n", + "Setup the training for ~30 epochs\n", + "
\n", + "\n", + "Tips:\n", + "- Set ``default_root_dir`` to store the logs and checkpoints\n", + "in a specific directory." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b9cf215b", + "metadata": { + "lines_to_next_cell": 2, + "title": "Setup trainer and check for errors." + }, + "outputs": [], + "source": [ + "\n", + "# fast_dev_run runs a single batch of data through the model to check for errors.\n", + "trainer = VSTrainer(accelerator=\"gpu\", devices=[GPU_ID], fast_dev_run=True)\n", + "\n", + "# trainer class takes the model and the data module as inputs.\n", + "trainer.fit(phase2fluor_model, datamodule=phase2fluor_data)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5f11d604", + "metadata": { + "lines_to_next_cell": 2 + }, + "outputs": [], + "source": [ + "\n", + "GPU_ID = 0\n", + "n_samples = len(phase2fluor_data.train_dataset)\n", + "steps_per_epoch = n_samples // BATCH_SIZE # steps per epoch.\n", + "n_epochs = 30\n", + "\n", + "trainer = VSTrainer(\n", + " accelerator=\"gpu\",\n", + " devices=[GPU_ID],\n", + " max_epochs=n_epochs,\n", + " log_every_n_steps=steps_per_epoch // 2,\n", + " # log losses and image samples 2 times per epoch.\n", + " default_root_dir=Path(\n", + " log_dir, \"phase2fluor\"\n", + " ), # lightning trainer transparently saves logs and model checkpoints in this directory.\n", + ")\n", + "\n", + "# Log graph\n", + "trainer.logger.log_graph(phase2fluor_model, phase2fluor_data.train_dataset[0][\"source\"])\n", + "# Launch training.\n", + "trainer.fit(phase2fluor_model, datamodule=phase2fluor_data)" + ] + }, + { + "cell_type": "markdown", + "id": "157a989c", + "metadata": { + "cell_marker": "\"\"\"" + }, + "source": [ + "
\n", + "Checkpoint 1\n", + "\n", + "Now the training has started,\n", + "we can come back after a while and evaluate the performance!\n", + "
" + ] + }, + { + "cell_type": "markdown", + "id": "d68dab93", + "metadata": { + "cell_marker": "\"\"\"", + "incorrectly_encoded_metadata": "id='1_fluor2phase'>", + "title": "\n", + "Checkpoint 2\n", + "Please summarize hyperparameters and performance of your models in [this google doc](https://docs.google.com/document/d/1hZWSVRvt9KJEdYu7ib-vFBqAVQRYL8cWaP_vFznu7D8/edit#heading=h.n5u485pmzv2z)\n", + "\n", + "Now that you have trained two models, let's think about the following questions:\n", + "- What is the information content of each channel in the dataset?\n", + "- How would you use image translation models?\n", + "- What can you try to improve the performance of each model?\n", + "\n", + "\n", + "
" + ] + }, + { + "cell_type": "markdown", + "id": "7d0c4204", + "metadata": { + "cell_marker": "\"\"\"", + "incorrectly_encoded_metadata": "id='3_tuning'>", + "title": "\n", + "Checkpoint 3\n", + "\n", + "Congratulations! You have trained several image translation models now!\n", + "Please document hyperparameters, snapshots of predictioons on validation set, and loss curves for your models in [this google doc](https://docs.google.com/document/d/1hZWSVRvt9KJEdYu7ib-vFBqAVQRYL8cWaP_vFznu7D8/edit#heading=h.n5u485pmzv2z)\n", + "" + ] + } + ], + "metadata": { + "jupytext": { + "cell_metadata_filter": "all", + "main_language": "python" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/tests/evaluation/test_evaluation_metrics.py b/tests/evaluation/test_evaluation_metrics.py index 78a7233f0..b0c36b430 100644 --- a/tests/evaluation/test_evaluation_metrics.py +++ b/tests/evaluation/test_evaluation_metrics.py @@ -135,6 +135,6 @@ def test_ms_ssim_25d(): img_pred = img + noise - 0.5 ssim_pred = ms_ssim_25d(img_pred, img) assert _is_within_unit(ssim_pred) - # Negative correlation should be zero when clamped + # clamped should be positive but very small ssim_inv = ms_ssim_25d(1 - img, img, clamp=True) - assert torch.allclose(ssim_inv, torch.tensor(0.0)) \ No newline at end of file + assert 0 < ssim_inv < 1e-3 diff --git a/viscy/evaluation/evaluation_metrics.py b/viscy/evaluation/evaluation_metrics.py index b1f8b4aed..7c6ec67a9 100644 --- a/viscy/evaluation/evaluation_metrics.py +++ b/viscy/evaluation/evaluation_metrics.py @@ -236,25 +236,26 @@ def ms_ssim_25d( :param torch.Tensor target: target images :param tuple[int, int] in_plane_window_size: kernel width and height, defaults to (11, 11) - :param bool clamp: clamp with ReLU for training stability when used in loss, + :param bool clamp: clamp to [1e-6, 1] for training stability when used in loss, defaults to False :param Sequence[float] betas: exponents of each resolution, defaults to (0.0448, 0.2856, 0.3001, 0.2363, 0.1333) :return torch.Tensor: multi-scale SSIM """ + base_min = 1e-4 mcs_list = [] for _ in range(len(betas)): ssim, contrast_sensitivity = ssim_25d( preds, target, in_plane_window_size, return_contrast_sensitivity=True ) if clamp: - contrast_sensitivity = torch.relu(contrast_sensitivity) + contrast_sensitivity = contrast_sensitivity.clamp(min=base_min) mcs_list.append(contrast_sensitivity) # do not downsample along depth preds = F.avg_pool3d(preds, (1, 2, 2)) target = F.avg_pool3d(target, (1, 2, 2)) if clamp: - ssim = torch.relu(ssim) + ssim = ssim.clamp(min=base_min) mcs_list[-1] = ssim mcs_stack = torch.stack(mcs_list) betas = torch.tensor(betas, device=mcs_stack.device).view(-1, 1) diff --git a/viscy/scripts/count_flops.py b/viscy/scripts/count_flops.py index c96f5a6b3..352da9c7a 100644 --- a/viscy/scripts/count_flops.py +++ b/viscy/scripts/count_flops.py @@ -13,6 +13,8 @@ "in_stack_depth": 5, "backbone": "convnextv2_tiny", "stem_kernel_size": (5, 4, 4), + "decoder_mode": "pixelshuffle", + "head_expansion_ratio": 4, }, ) diff --git a/viscy/scripts/network_diagram.py b/viscy/scripts/network_diagram.py index 94b56f1bf..419d69a45 100644 --- a/viscy/scripts/network_diagram.py +++ b/viscy/scripts/network_diagram.py @@ -6,18 +6,15 @@ # %% 2D UNet model = VSUNet( architecture="2D", - model_config={ - "in_channels": 2, - "out_channels": 1, - }, + model_config={"in_channels": 1, "out_channels": 2}, ) - +# %% model_graph = draw_graph( model, model.example_input_array, graph_name="2D UNet", roll=True, - depth=2, + depth=4, # graph_dir="LR", # save_graph=True, ) @@ -75,12 +72,12 @@ model = VSUNet( architecture="2.2D", model_config={ - "in_channels": 2, - "out_channels": 1, + "in_channels": 1, + "out_channels": 2, "in_stack_depth": 9, "backbone": "convnextv2_tiny", - "stem_kernel_size": (3, 2, 2), - "decoder_mode": "deconv", + "decoder_mode": "pixelshuffle", + "stem_kernel_size": (3, 4, 4), }, ) diff --git a/viscy/unet/networks/Unet21D.py b/viscy/unet/networks/Unet21D.py index 4cdbf4e8f..7c32e34b8 100644 --- a/viscy/unet/networks/Unet21D.py +++ b/viscy/unet/networks/Unet21D.py @@ -1,23 +1,67 @@ -from typing import Literal, Sequence +from typing import Callable, Literal, Optional, Sequence, Union import timm import torch -from monai.networks.blocks import ResidualUnit, UpSample +from monai.networks.blocks import Convolution, ResidualUnit, UpSample from monai.networks.blocks.dynunet_block import get_conv_layer +from monai.networks.utils import normal_init from torch import nn -def _get_convnext_stage(in_channels: int, out_channels: int, depth: int): - return timm.models.convnext.ConvNeXtStage( +def icnr_init( + conv: nn.Module, + upsample_factor: int, + upsample_dims: int, + init=nn.init.kaiming_normal_, +): + """ + ICNR initialization for 2D/3D kernels adapted from Aitken et al.,2017 , + "Checkerboard artifact free sub-pixel convolution". + + Adapted from MONAI v1.2.0, added support for upsampling dimensions + that are not the same as the kernel dimension. + + :param conv: convolution layer + :param upsample_factor: upsample factor + :param upsample_dims: upsample dimensions, 2 or 3 + :param init: initialization function + """ + out_channels, in_channels, *dims = conv.weight.shape + scale_factor = upsample_factor**upsample_dims + + oc2 = int(out_channels / scale_factor) + + kernel = torch.zeros([oc2, in_channels] + dims) + kernel = init(kernel) + kernel = kernel.transpose(0, 1) + kernel = kernel.reshape(oc2, in_channels, -1) + kernel = kernel.repeat(1, 1, scale_factor) + kernel = kernel.reshape([in_channels, out_channels] + dims) + kernel = kernel.transpose(0, 1) + conv.weight.data.copy_(kernel) + + +def _get_convnext_stage( + in_channels: int, + out_channels: int, + depth: int, + upsample_factor: Optional[int] = None, +) -> nn.Module: + stage = timm.models.convnext.ConvNeXtStage( in_chs=in_channels, out_chs=out_channels, stride=1, depth=depth, ls_init_value=None, + conv_mlp=True, use_grn=True, norm_layer=timm.layers.LayerNorm2d, norm_layer_cl=timm.layers.LayerNorm, ) + stage.apply(timm.models.convnext._init_weights) + if upsample_factor: + icnr_init(stage.blocks[-1].mlp.fc2, upsample_factor, upsample_dims=2) + return stage class Conv21dStem(nn.Module): @@ -51,11 +95,13 @@ class Unet2dUpStage(nn.Module): def __init__( self, in_channels: int, + skip_channels: int, out_channels: int, scale_factor: int, mode: Literal["deconv", "pixelshuffle"], conv_blocks: int, norm_name: str, + upsample_pre_conv: Optional[Union[Literal["default"], Callable]], ) -> None: super().__init__() spatial_dims = 2 @@ -81,19 +127,23 @@ def __init__( nn.Conv2d(in_channels, out_channels, kernel_size=(1, 1)), ) elif mode == "pixelshuffle": + mid_channels = in_channels // scale_factor**2 self.upsample = UpSample( spatial_dims=spatial_dims, in_channels=in_channels, - out_channels=out_channels, + out_channels=mid_channels, scale_factor=scale_factor, mode=mode, - pre_conv="default", - apply_pad_pool=True, + pre_conv=upsample_pre_conv, + apply_pad_pool=False, ) + conv_weight_init_factor = None if upsample_pre_conv else scale_factor self.conv = _get_convnext_stage( - out_channels + out_channels, out_channels, conv_blocks + mid_channels + skip_channels, + out_channels, + conv_blocks, + upsample_factor=conv_weight_init_factor, ) - self.conv.apply(timm.models.convnext._init_weights) def forward(self, inp: torch.Tensor, skip: torch.Tensor) -> torch.Tensor: """ @@ -112,26 +162,45 @@ def __init__( in_channels: int, out_channels: int, out_stack_depth: int, + expansion_ratio: int, + pool: bool, ) -> None: super().__init__() - self.norm = timm.layers.LayerNorm2d(num_channels=in_channels) - self.gelu = nn.GELU() - self.conv = nn.Conv3d( - in_channels // out_stack_depth, - out_channels, - kernel_size=3, - stride=1, - padding=1, + first_scale = 2 + self.upsample = UpSample( + spatial_dims=2, + in_channels=in_channels, + out_channels=in_channels // first_scale**2, + scale_factor=first_scale, + mode="pixelshuffle", + pre_conv=None, + apply_pad_pool=pool, + ) + mid_channels = out_channels * expansion_ratio * 2**2 + self.conv = nn.Sequential( + Convolution( + spatial_dims=3, + in_channels=in_channels // first_scale**2 // (out_stack_depth + 2), + out_channels=mid_channels, + kernel_size=3, + padding=(0, 1, 1), + ), + nn.Conv3d(mid_channels, out_channels * 2**2, 1), ) + normal_init(self.conv[0]) + icnr_init(self.conv[-1], 2, upsample_dims=2) + self.out = nn.PixelShuffle(2) self.out_stack_depth = out_stack_depth def forward(self, x: torch.Tensor) -> torch.Tensor: - x = self.norm(x) - x = self.gelu(x) + x = self.upsample(x) + d = self.out_stack_depth + 2 b, c, h, w = x.shape - x = x.reshape((b, c // self.out_stack_depth, self.out_stack_depth, h, w)) + x = x.reshape((b, c // d, d, h, w)) x = self.conv(x) - return x + x = x.transpose(1, 2) + x = self.out(x) + return x.transpose(1, 2) class UnsqueezeHead(nn.Module): @@ -149,44 +218,35 @@ class Unet2dDecoder(nn.Module): def __init__( self, num_channels: list[int], - out_channels: int, norm_name: str, mode: Literal["deconv", "pixelshuffle"], conv_blocks: int, strides: list[int], + upsample_pre_conv: Optional[Union[Literal["default"], Callable]], ) -> None: super().__init__() self.decoder_stages = nn.ModuleList([]) - stages = len(num_channels) - num_channels.append(num_channels[-1]) + stages = len(num_channels) - 1 for i in range(stages): - stride = strides[i] stage = Unet2dUpStage( in_channels=num_channels[i], + skip_channels=num_channels[i] // 2, out_channels=num_channels[i + 1], - scale_factor=stride, + scale_factor=strides[i], mode=mode, conv_blocks=conv_blocks, norm_name=norm_name, + upsample_pre_conv=upsample_pre_conv, ) self.decoder_stages.append(stage) - self.head = UpSample( - spatial_dims=2, - in_channels=num_channels[-1], - out_channels=out_channels, - scale_factor=strides[-1], - mode=mode, - pre_conv="default", - apply_pad_pool=False, - ) def forward(self, features: Sequence[torch.Tensor]) -> torch.Tensor: feat = features[0] # padding features.append(None) - for skip, stage in zip(features[1:], self.decoder_stages[:-1]): + for skip, stage in zip(features[1:], self.decoder_stages): feat = stage(feat, skip) - return self.head(feat) + return feat class Unet21d(nn.Module): @@ -202,6 +262,9 @@ def __init__( decoder_mode: Literal["deconv", "pixelshuffle"] = "pixelshuffle", decoder_conv_blocks: int = 2, decoder_norm_layer: str = "instance", + decoder_upsample_pre_conv: bool = False, + head_pool: bool = False, + head_expansion_ratio: int = 4, drop_path_rate: float = 0.0, ) -> None: super().__init__() @@ -216,6 +279,12 @@ def __init__( f"the same as `input_stack_depth` ({in_stack_depth}), " f"but got {out_stack_depth}." ) + if not (in_stack_depth == out_stack_depth or out_stack_depth == 1): + raise ValueError( + "`out_stack_depth` must be either 1 or " + f"the same as `input_stack_depth` ({in_stack_depth}), " + f"but got {out_stack_depth}." + ) multi_scale_encoder = timm.create_model( backbone, pretrained=pretrained, @@ -231,24 +300,27 @@ def __init__( ) decoder_channels = num_channels decoder_channels.reverse() - if out_stack_depth == 1: - decoder_out_channels = out_channels - self.head = UnsqueezeHead() - else: - decoder_out_channels = ( - out_stack_depth * decoder_channels[-1] // stem_kernel_size[-1] ** 2 - ) - self.head = PixelToVoxelHead( - decoder_out_channels, out_channels, out_stack_depth - ) + decoder_channels[-1] = ( + (out_stack_depth + 2) * out_channels * 2**2 * head_expansion_ratio + ) self.decoder = Unet2dDecoder( decoder_channels, - decoder_out_channels, norm_name=decoder_norm_layer, mode=decoder_mode, conv_blocks=decoder_conv_blocks, - strides=[2] * len(num_channels) + [stem_kernel_size[-1]], + strides=[2] * (len(num_channels) - 1) + [stem_kernel_size[-1]], + upsample_pre_conv="default" if decoder_upsample_pre_conv else None, ) + if out_stack_depth == 1: + self.head = UnsqueezeHead() + else: + self.head = PixelToVoxelHead( + decoder_channels[-1], + out_channels, + out_stack_depth, + head_expansion_ratio, + pool=head_pool, + ) self.out_stack_depth = out_stack_depth @property