Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,6 @@
import numpy as np
import os
import torch
import torch.nn as nn
import torch.optim as optim
from skimage import io

Expand Down
227 changes: 183 additions & 44 deletions src/aind_exaspim_image_compression/machine_learning/unet3d.py
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,14 @@ class UNet(nn.Module):
Final 1x1x1 convolution mapping features to the output channel.
"""

def __init__(self, width_multiplier=1, trilinear=True, residual=True):
def __init__(
self,
width_multiplier=1,
trilinear=True,
residual=True,
maxblurpool=False,
remove_top_skip=False,

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The default residual=True would defeat the purpose of remove_top_skip=True. If a user enables remove_top_skip=True to reproduce N2V2's hot-pixel/checkerboard suppression but leaves residual at its default True; pixel noise still passes straight through via the residual add, so the architectural change does nothing and the experiment silently evaluates the wrong hypothesis. At minimum document the interaction; consider warning or disallowing remove_top_skip=True with residual=True.

):
"""
Instantiates a UNet object.

Expand All @@ -56,38 +63,87 @@ def __init__(self, width_multiplier=1, trilinear=True, residual=True):
Default is True.
"""
# Call parent class
super(UNet, self).__init__()
super().__init__()

if (
isinstance(width_multiplier, bool)
or not isinstance(width_multiplier, Real)
isinstance(width_multiplier, Real)

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why change this condition? isinstance(width_multiplier, Real) is True for every ordinary int or float so attempting to construct this class will always raise a ValueError.

or width_multiplier < 1
or not float(width_multiplier).is_integer()
):
raise ValueError("width_multiplier must be a positive integer")

# Initializations
_channels = (32, 64, 128, 256, 512)
base_channels = (32, 64, 128, 256, 512)
factor = 2 if trilinear else 1

# Instance attributes
self.width_multiplier = int(width_multiplier)
self.channels = [c * self.width_multiplier for c in _channels]
self.channels = [
c * self.width_multiplier
for c in base_channels
]

self.trilinear = trilinear
self.residual = residual
self.maxblurpool = maxblurpool
self.remove_top_skip = remove_top_skip

# Contracting layers
# Encoder
self.inc = DoubleConv(1, self.channels[0])
self.down1 = Down(self.channels[0], self.channels[1])
self.down2 = Down(self.channels[1], self.channels[2])
self.down3 = Down(self.channels[2], self.channels[3])
self.down4 = Down(self.channels[3], self.channels[4] // factor)

# Expanding layers
self.up1 = Up(self.channels[4], self.channels[3] // factor, trilinear)
self.up2 = Up(self.channels[3], self.channels[2] // factor, trilinear)
self.up3 = Up(self.channels[2], self.channels[1] // factor, trilinear)
self.up4 = Up(self.channels[1], self.channels[0], trilinear)

self.down1 = Down(
self.channels[0],
self.channels[1],
maxblurpool=maxblurpool,
)

self.down2 = Down(
self.channels[1],
self.channels[2],
maxblurpool=maxblurpool,
)

self.down3 = Down(
self.channels[2],
self.channels[3],
maxblurpool=maxblurpool,
)

self.down4 = Down(
self.channels[3],
self.channels[4] // factor,
maxblurpool=maxblurpool,
)

# Decoder
self.up1 = Up(
self.channels[4],
self.channels[3] // factor,
trilinear=trilinear,
use_skip=True,
)

self.up2 = Up(
self.channels[3],
self.channels[2] // factor,
trilinear=trilinear,
use_skip=True,
)

self.up3 = Up(
self.channels[2],
self.channels[1] // factor,
trilinear=trilinear,
use_skip=True,
)

self.up4 = Up(
self.channels[1],
self.channels[0],
trilinear=trilinear,
use_skip=not remove_top_skip,
)

self.outc = OutConv(self.channels[0], 1)

@property
Expand All @@ -97,6 +153,8 @@ def config(self):
"width_multiplier": self.width_multiplier,
"trilinear": self.trilinear,
"residual": self.residual,
"maxblurpool": self.maxblurpool,
"remove_top_skip": self.remove_top_skip,
}

def forward(self, x):
Expand All @@ -114,23 +172,29 @@ def forward(self, x):
Output tensor with shape (B, 1, D, H, W), representing the
denoised image.
"""
# Contracting layers
# Encoder
x1 = self.inc(x)
x2 = self.down1(x1)
x3 = self.down2(x2)
x4 = self.down3(x3)
x5 = self.down4(x4)

# Expanding layers
# Decoder
d = self.up1(x5, x4)
d = self.up2(d, x3)
d = self.up3(d, x2)
d = self.up4(d, x1)

if self.remove_top_skip:
d = self.up4(d)
else:
d = self.up4(d, x1)

logits = self.outc(d)

# Residual denoising: predict the correction added to the input
if self.residual:
return x + logits

return logits


Expand Down Expand Up @@ -211,15 +275,14 @@ def forward(self, x):
class Down(nn.Module):
"""
A downsampling module for a 3D U-Net.

Attributes
----------
maxpool_conv : nn.Sequential
Sequential module containing a MaxPool3d layer followed by a
DoubleConv block.
"""

def __init__(self, in_channels, out_channels):
def __init__(
self,
in_channels,
out_channels,
maxblurpool=False,
):
"""
Instantiates a Down object.

Expand All @@ -229,13 +292,23 @@ def __init__(self, in_channels, out_channels):
Number of input channels to this module.
out_channels : int
Number of output channels produced by this module.
maxblurpool : bool, optional
True if max-blur pooling should be used to downsample. Default is
False.
"""
# Call parent class
super().__init__()

# Initializations
if maxblurpool:
downsample = MaxBlurPool3D(in_channels)
else:
downsample = nn.MaxPool3d(2)

# Instance attributes
self.maxpool_conv = nn.Sequential(
nn.MaxPool3d(2), DoubleConv(in_channels, out_channels)
downsample,
DoubleConv(in_channels, out_channels),
)

def forward(self, x):
Expand Down Expand Up @@ -269,7 +342,13 @@ class Up(nn.Module):
connection.
"""

def __init__(self, in_channels, out_channels, trilinear=True):
def __init__(
self,
in_channels,
out_channels,
trilinear=True,
use_skip=True,
):
"""
Instantiates an Up object.

Expand All @@ -287,20 +366,32 @@ def __init__(self, in_channels, out_channels, trilinear=True):
super().__init__()

# Instance attributes
self.use_skip = use_skip

if trilinear:
self.up = nn.Upsample(
scale_factor=2, mode="trilinear", align_corners=True
)
self.conv = DoubleConv(
in_channels, out_channels, mid_channels=in_channels // 2

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

When loading a checkpoint created before this patch with trilinear=True, omitting
mid_channels=in_channels // 2 changes parameter shapes throughout up1–up3. Thus, load_state_dict in load_model will reject existing checkpoints. This architecture should probably be a separate class, e.g., UNetN2V2 so we can preserve backwards compatibility and do A/B testing more easily.

scale_factor=2,
mode="trilinear",
align_corners=True,
)
else:
self.up = nn.ConvTranspose3d(
in_channels, in_channels // 2, kernel_size=2, stride=2
in_channels,
in_channels // 2,
kernel_size=2,
stride=2,
)
self.conv = DoubleConv(in_channels, out_channels)

def forward(self, x1, x2):
conv_in_channels = (
in_channels if use_skip else in_channels // 2
)

self.conv = DoubleConv(
conv_in_channels,
out_channels,
)

def forward(self, x1, x2=None):
"""
Forward pass of the upsampling block in a 3D U-Net.

Expand All @@ -309,7 +400,7 @@ def forward(self, x1, x2):
x1 : torch.Tensor
Input tensor from the previous decoder layer with shape
(B, C1, D, H1, W1).
x2 : torch.Tensor
x2 : torch.Tensor, optional
Skip connection tensor from the encoder path with shape
(B, C2, D, H2, W2).

Expand All @@ -321,17 +412,65 @@ def forward(self, x1, x2):
(B, out_channels, D, H2, W2).
"""
x1 = self.up(x1)
diffY = x2.size()[2] - x1.size()[2]
diffX = x2.size()[3] - x1.size()[3]

x1 = F.pad(
x1,
[diffX // 2, diffX - diffX // 2, diffY // 2, diffY - diffY // 2],
)
x = torch.cat([x2, x1], dim=1)
if self.use_skip:
diff_z = x2.size(2) - x1.size(2)
diff_y = x2.size(3) - x1.size(3)
diff_x = x2.size(4) - x1.size(4)

x1 = F.pad(
x1,
[
diff_x // 2,
diff_x - diff_x // 2,
diff_y // 2,
diff_y - diff_y // 2,
diff_z // 2,
diff_z - diff_z // 2,
],
)

x = torch.cat([x2, x1], dim=1)

else:

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

When remove_top_skip=True and any input spatial dimension is odd, the final upsampling doubles floor(input_size / 2) without aligning it to x1, leaving the output one voxel smaller in that dimension. With the default residual mode, x + logits then raises a size-mismatch error; without residual mode, the output tensor still has the wrong shape.

Example:
Using an odd-sized input (1, 1, 17, 19, 21) with remove_top_skip=True:

residual=False:
input = (1, 1, 17, 19, 21)
output = (1, 1, 16, 18, 20)

residual=True:
RuntimeError: The size of tensor a (21) must match the size of tensor b (20)

The encoder’s first pooling produces (8, 9, 10). The final decoder upsampling doubles that to (16, 18, 20), losing one voxel along each odd dimension.

x = x1

return self.conv(x)


class MaxBlurPool3D(nn.Module):

def __init__(self, channels):
super().__init__()

self.pool = nn.MaxPool3d(2, stride=1)

kernel = torch.tensor([1., 2., 1.])
kernel = (
kernel[:, None, None]
* kernel[None, :, None]
* kernel[None, None, :]
)
kernel /= kernel.sum()

self.register_buffer(
"kernel",
kernel[None, None].repeat(channels, 1, 1, 1, 1),
)
self.channels = channels

def forward(self, x):
x = self.pool(x)
x = F.conv3d(
x,
self.kernel,
stride=2,
padding=1,

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

MaxBlurPool3D's blur conv uses zero padding, so the normalized [1,2,1] binomial kernel systematically attenuates border voxels at every encoder scale (reference BlurPool implementations use reflect/replicate padding e.g., https://github.com/seominseok0429/Impact-of-Aliasing-on-Generalization-in-Deep-Convolutional-Networks-pytorch/blob/main/blurpool.py). For example, constant input 1.0 yields interior 1.0 but faces 0.75, edges 0.56, corners 0.42, vs 1.0 everywhere for the nn.MaxPool3d(2) it replaces.

groups=self.channels,
)
return x


class OutConv(nn.Module):
"""
Final output convolution layer for a 3D U-Net.
Expand Down
Loading