StarNetPyTorch/starnet_v1_pytorch.py

import pickle
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
import tifffile as tiff
import torch
import torch.nn.functional as F

from model import Discriminator, Generator


class StarNet:
    def __init__(

        self,

        mode,

        window_size=512,

        stride=256,

        lr=1e-4,

        train_folder="./train/",

        batch_size=1,

        device=None,

    ):
        if mode not in ("RGB", "Greyscale"):
            raise ValueError("Mode should be either RGB or Greyscale")
        if window_size % 256 != 0:
            raise ValueError("window_size must be divisible by 256")
        if stride <= 0 or stride > window_size:
            raise ValueError("stride must be greater than zero and no larger than window_size")
        if (window_size - stride) % 2 != 0:
            raise ValueError("window_size - stride must be even")

        self.mode = mode
        self.input_channels = 3 if mode == "RGB" else 1
        self.window_size = window_size
        self.stride = stride
        self.lr = lr
        self.train_folder = Path(train_folder)
        self.batch_size = batch_size
        self.device = torch.device(device or ("cuda" if torch.cuda.is_available() else "cpu"))

        self.history = {}
        self._ema = 0.9999
        self.weights = []
        self.original = []
        self.starless = []
        self.iters_per_epoch = 0

        self.G = None
        self.D = None
        self.gen_optimizer = None
        self.dis_optimizer = None

    def __str__(self):
        return "StarNet PyTorch instance"

    def load_training_dataset(self):
        original_dir = self.train_folder / "original"
        starless_dir = self.train_folder / "starless"
        original_files = sorted(path.name for path in original_dir.glob("*.tif"))
        starless_files = sorted(path.name for path in starless_dir.glob("*.tif"))

        if not original_files or not starless_files:
            raise ValueError("No training data found in %s" % self.train_folder)
        if original_files != starless_files:
            raise ValueError(
                "Corresponding names in the original and starless folders must be equal"
            )

        print("Total training images found: %d" % len(original_files))
        self.original = [tiff.imread(original_dir / name) for name in original_files]
        self.starless = [tiff.imread(starless_dir / name) for name in starless_files]
        self.weights = []
        total_pixels = 0

        for name, original, starless in zip(original_files, self.original, self.starless):
            if original.shape != starless.shape:
                raise ValueError("Image sizes are not equal for %s" % name)
            if original.shape[0] < self.window_size or original.shape[1] < self.window_size:
                raise ValueError(
                    "%s is smaller than the configured window size %d"
                    % (name, self.window_size)
                )
            pixels = original.shape[0] * original.shape[1]
            total_pixels += pixels
            self.weights.append(pixels)

        self.iters_per_epoch = max(1, total_pixels // (self.window_size * self.window_size))
        self.weights = np.asarray(self.weights, dtype=np.float64)
        self.weights /= self.weights.sum()

        print("Total size of training images: %.2f MP" % (total_pixels / 1e6))
        print("One epoch is set to %d iterations" % self.iters_per_epoch)
        print("Training dataset has been successfully loaded!")

    def _generator(self, m=64):
        if m != 64:
            raise ValueError("The translated architecture supports m=64")
        return Generator(input_channels=self.input_channels, batch_norm_eps=1e-3)

    def _discriminator(self):
        return Discriminator(input_channels=self.input_channels, batch_norm_eps=1e-3)

    def _weight_paths(self, weights):
        path = Path(weights)
        
        if path.suffix.lower() in (".pth", ".pt"):
            return path, None
        return (
            Path(str(path) + "_G_" + self.mode + ".pth"),
            Path(str(path) + "_D_" + self.mode + ".pth"),
        )

    @staticmethod
    def _load_state(module, path, device):
        try:
            checkpoint = torch.load(path, map_location=device, weights_only=True)
        except TypeError:  # PyTorch < 2.0 compatibility
            checkpoint = torch.load(path, map_location=device)
        state_dict = checkpoint.get("state_dict", checkpoint)
        module.load_state_dict(state_dict)

    def load_model(self, weights=None, history=None):
        self.G = self._generator().to(self.device)
        self.D = self._discriminator().to(self.device)
        self.gen_optimizer = torch.optim.Adam(self.G.parameters(), lr=self.lr)
        self.dis_optimizer = torch.optim.Adam(self.D.parameters(), lr=self.lr / 4)

        if weights:
            generator_path, discriminator_path = self._weight_paths(weights)
            if not generator_path.exists():
                raise FileNotFoundError("Generator weights were not found at %s" % generator_path)
            self._load_state(self.G, generator_path, self.device)

            if discriminator_path and discriminator_path.exists():
                self._load_state(self.D, discriminator_path, self.device)
            elif discriminator_path:
                print("Discriminator weights were not found; using initialized weights.")

        if history:
            history_path = Path(str(history) + "_" + self.mode + ".pkl")
            with history_path.open("rb") as handle:
                self.history = pickle.load(handle)

        self.G.eval()
        self.D.eval()
        return self.G

    def initialize_model(self):
        return self.load_model()

    @staticmethod
    def _ramp(x):
        if torch.is_tensor(x):
            return torch.clamp(x, 0.0, 1.0)
        return np.clip(x, 0.0, 1.0)

    @staticmethod
    def _image_scale(image):
        if image.dtype == np.uint8:
            return 255.0
        if image.dtype == np.uint16:
            return 65535.0
        if np.issubdtype(image.dtype, np.floating):
            return 1.0
        raise ValueError("Unsupported training image dtype: %s" % image.dtype)

    def _get_sample(self, r, h, w, type):
        if type not in ("original", "starless"):
            raise ValueError("type must be original or starless")
        source = self.original if type == "original" else self.starless
        sample = source[r][h:h + self.window_size, w:w + self.window_size]
        return sample.astype(np.float32) / self._image_scale(sample)

    def _augmentator(self, original, starless):
        original = np.array(original, copy=True)
        starless = np.array(starless, copy=True)

        if np.random.rand() < 0.5:
            original = np.flip(original, axis=1)
            starless = np.flip(starless, axis=1)
        if np.random.rand() < 0.5:
            original = np.flip(original, axis=0)
            starless = np.flip(starless, axis=0)
        if np.random.rand() < 0.5:
            rotations = np.random.randint(1, 4)
            original = np.rot90(original, rotations, axes=(1, 0))
            starless = np.rot90(starless, rotations, axes=(1, 0))

        if self.mode == "RGB":
            if original.ndim != 3 or original.shape[2] < 3:
                raise ValueError("RGB training images must contain at least three channels")
            original = original[:, :, :3]
            starless = starless[:, :, :3]
            if np.random.rand() < 0.7:
                channel = np.random.randint(3)
                minimum = min(float(original.min()), float(starless.min()))
                offset = np.random.rand() * 0.25 - np.random.rand() * minimum
                original[:, :, channel] += offset * (1.0 - original[:, :, channel])
                starless[:, :, channel] += offset * (1.0 - starless[:, :, channel])
            if np.random.rand() < 0.7:
                sequence = np.random.permutation(3)
                original = original[:, :, sequence]
                starless = starless[:, :, sequence]
        else:
            if np.random.rand() < 0.7:
                minimum = min(float(original.min()), float(starless.min()))
                offset = np.random.rand() * 0.25 - np.random.rand() * minimum
                original += offset * (1.0 - original)
                starless += offset * (1.0 - starless)
            if original.ndim == 3:
                channel = np.random.randint(min(3, original.shape[2]))
                original = original[:, :, channel]
                starless = starless[:, :, channel]
            original = original[:, :, None]
            starless = starless[:, :, None]

        return (
            np.ascontiguousarray(np.clip(original, 0.0, 1.0)),
            np.ascontiguousarray(np.clip(starless, 0.0, 1.0)),
        )

    def generate_input(self, iterations=1, augmentation=False):
        if not self.original:
            raise RuntimeError("Training dataset was not loaded")

        original_batch = None
        starless_batch = None
        for _ in range(iterations):
            original_batch = np.zeros(
                (
                    self.batch_size,
                    self.window_size,
                    self.window_size,
                    self.input_channels,
                ),
                dtype=np.float32,
            )
            starless_batch = np.zeros_like(original_batch)

            for index in range(self.batch_size):
                image_index = int(
                    np.random.choice(len(self.original), p=np.asarray(self.weights))
                )
                image_h, image_w = self.original[image_index].shape[:2]
                top = np.random.randint(0, image_h - self.window_size + 1)
                left = np.random.randint(0, image_w - self.window_size + 1)
                original = self._get_sample(image_index, top, left, "original")
                starless = self._get_sample(image_index, top, left, "starless")

                if augmentation:
                    original, starless = self._augmentator(original, starless)
                elif self.mode == "Greyscale":
                    if original.ndim == 3:
                        channel = np.random.randint(min(3, original.shape[2]))
                        original = original[:, :, channel]
                        starless = starless[:, :, channel]
                    original = original[:, :, None]
                    starless = starless[:, :, None]
                else:
                    original = original[:, :, :3]
                    starless = starless[:, :, :3]

                original_batch[index] = original
                starless_batch[index] = starless

        return original_batch, starless_batch

    @staticmethod
    def _set_requires_grad(module, enabled):
        for parameter in module.parameters():
            parameter.requires_grad_(enabled)

    def _losses(self, original, target):
        generated = self.G(original)
        real_outputs = self.D(target)
        fake_outputs = self.D(generated)

        real_features = real_outputs[:-1]
        fake_features = fake_outputs[:-1]
        predict_real = real_outputs[-1]
        predict_fake = fake_outputs[-1]

        dis_loss = -(
            torch.log(predict_real + 1e-8) + torch.log(1.0 - predict_fake + 1e-8)
        ).mean()
        gen_loss_gan = -torch.log(predict_fake + 1e-8).mean()
        perceptual = [
            torch.mean(torch.abs(fake - real))
            for fake, real in zip(fake_features, real_features)
        ]
        gen_l1 = F.l1_loss(generated, target)
        gen_loss = (
            gen_loss_gan * 0.1
            + perceptual[0] * 0.1
            + sum(loss * 10.0 for loss in perceptual[1:])
            + gen_l1 * 100.0
        )

        metrics = {
            "dis_loss": dis_loss,
            "gen_loss_GAN": gen_loss_gan,
            "gen_p1": perceptual[0],
            "gen_p2": perceptual[1],
            "gen_p3": perceptual[2],
            "gen_p4": perceptual[3],
            "gen_p5": perceptual[4],
            "gen_p6": perceptual[5],
            "gen_p7": perceptual[6],
            "gen_p8": perceptual[7],
            "gen_L1": gen_l1 * 100.0,
            "total": gen_loss,
        }
        return gen_loss, dis_loss, metrics

    def _train_step(self, original, target):
        self.G.train()
        self.D.train()
        self.gen_optimizer.zero_grad(set_to_none=True)
        self.dis_optimizer.zero_grad(set_to_none=True)

        gen_loss, dis_loss, metrics = self._losses(original, target)
        generator_gradients = torch.autograd.grad(
            gen_loss,
            tuple(self.G.parameters()),
            retain_graph=True,
        )
        discriminator_gradients = torch.autograd.grad(
            dis_loss,
            tuple(self.D.parameters()),
        )

        for parameter, gradient in zip(self.G.parameters(), generator_gradients):
            parameter.grad = gradient
        for parameter, gradient in zip(self.D.parameters(), discriminator_gradients):
            parameter.grad = gradient

        self.gen_optimizer.step()
        self.dis_optimizer.step()
        return {name: float(value.detach().cpu()) for name, value in metrics.items()}

    def _update_history(self, metrics):
        for name, value in metrics.items():
            if name in self.history and self.history[name]:
                value = value * (1.0 - self._ema) + self.history[name][-1] * self._ema
                self.history[name].append(value)
            else:
                self.history[name] = [value]

    def _plot_progress(self, original, target):
        self.G.eval()
        with torch.inference_mode():
            generated = self.G(original)

        original = ((original[0] + 1.0) / 2.0).detach().cpu().numpy()
        generated = ((generated[0] + 1.0) / 2.0).detach().cpu().numpy()
        target = ((target[0] + 1.0) / 2.0).detach().cpu().numpy()

        plt.close()
        _, axes = plt.subplots(1, 3, sharex=True, figsize=(16.5, 5.5))
        titles = ("Original", "Starless", "Target")
        for axis, title, image in zip(axes, titles, (original, generated, target)):
            if self.mode == "RGB":
                axis.imshow(np.transpose(image, (1, 2, 0)))
            else:
                axis.imshow(image[0], cmap="gray", vmin=0, vmax=1)
            axis.set_title(title)
        plt.pause(0.001)
        self.G.train()

    def train(

        self,

        epochs,

        augmentation=True,

        plot_progress=False,

        plot_interval=50,

        save_backups=True,

        warm_up=False,

    ):
        if not self.original:
            raise RuntimeError("Training dataset was not loaded; call load_training_dataset()")
        if self.G is None or self.D is None:
            self.initialize_model()

        for epoch in range(epochs):
            for iteration in range(self.iters_per_epoch):
                original, target = self.generate_input(augmentation=augmentation)
                original = torch.from_numpy(original.transpose(0, 3, 1, 2)).to(self.device)
                target = torch.from_numpy(target.transpose(0, 3, 1, 2)).to(self.device)
                original = original * 2.0 - 1.0
                target = target * 2.0 - 1.0
                if warm_up:
                    target = original

                if plot_progress and iteration % plot_interval == 0:
                    self._plot_progress(original, target)

                metrics = self._train_step(original, target)
                self._update_history(metrics)
                print(
                    "\rEpoch: %d. Iteration %d / %d Loss %f    "
                    % (epoch, iteration, self.iters_per_epoch, self.history["total"][-1]),
                    end="",
                )

            print()
            if save_backups:
                suffix = "even" if epoch % 2 == 0 else "odd"
                self._save_network(self.G, "./starnet_backup_G_%s.pth" % suffix)
                self._save_network(self.D, "./starnet_backup_D_%s.pth" % suffix)

        if plot_progress:
            plt.close()
        self.G.eval()
        self.D.eval()
        return self.history

    def plot_history(self, last=None):
        if not self.history:
            raise RuntimeError("Empty training history, nothing to plot")

        keys = list(self.history.keys())
        columns = 3
        rows = int(np.ceil(len(keys) / columns))
        figure, axes = plt.subplots(rows, columns, sharex=True, figsize=(16, 3.5 * rows))
        axes = np.atleast_1d(axes).reshape(-1)
        for axis, key in zip(axes, keys):
            values = self.history[key][-last:] if last else self.history[key]
            axis.plot(values)
            axis.set_title(key)
        for axis in axes[len(keys):]:
            axis.set_visible(False)
        figure.tight_layout()
        return figure, axes

    @staticmethod
    def _save_network(network, path):
        path = Path(path)
        path.parent.mkdir(parents=True, exist_ok=True)
        torch.save({"state_dict": network.state_dict()}, path)

    def save_model(self, weights_filename, history_filename=None):
        if self.G is None or self.D is None:
            raise RuntimeError("Model has not been initialized")

        generator_path = Path(str(weights_filename) + "_G_" + self.mode + ".pth")
        discriminator_path = Path(str(weights_filename) + "_D_" + self.mode + ".pth")
        self._save_network(self.G, generator_path)
        self._save_network(self.D, discriminator_path)

        if history_filename:
            history_path = Path(str(history_filename) + "_" + self.mode + ".pkl")
            history_path.parent.mkdir(parents=True, exist_ok=True)
            with history_path.open("wb") as handle:
                pickle.dump(self.history, handle)
        return generator_path, discriminator_path

    def _read_image(self, input_path):
        data = tiff.imread(input_path)
        if data.ndim > 3:
            layer = input("Tiff has %d layers, please enter layer to process: " % data.shape[0])
            data = data[int(layer)]

        if data.dtype == np.dtype("uint16"):
            scale = 65535.0
        elif data.dtype == np.dtype("uint8"):
            scale = 255.0
        else:
            raise ValueError("Unknown image dtype: %s" % data.dtype)

        image = (data / scale).astype(np.float32)
        if self.mode == "Greyscale":
            if image.ndim != 2:
                raise ValueError("You loaded Greyscale model, but the image is RGB")
            image = image[:, :, None]
        else:
            if image.ndim != 3:
                raise ValueError("You loaded RGB model, but the image is Greyscale")
            if image.shape[2] == 4:
                print("Input image has 4 channels. Removing Alpha-Channel")
                image = image[:, :, :3]
            if image.shape[2] != 3:
                raise ValueError("RGB input must have 3 or 4 channels")
        return image, data.dtype, scale

    @staticmethod
    def _imwrite(path, image):
        if hasattr(tiff, "imwrite"):
            tiff.imwrite(path, image)
        else:
            tiff.imsave(path, image)

    def transform(self, in_name, out_name):
        if self.G is None:
            raise RuntimeError("Call load_model() before transform()")

        self.G.eval()
        image, input_dtype, scale = self._read_image(in_name)
        offset = (self.window_size - self.stride) // 2
        h, w, _ = image.shape
        ith = h // self.stride + 1
        itw = w // self.stride + 1
        dh = ith * self.stride - h
        dw = itw * self.stride - w

        image = np.concatenate((image, image[(h - dh):, :, :]), axis=0)
        image = np.concatenate((image, image[:, (w - dw):, :]), axis=1)
        padded_h, padded_w, _ = image.shape
        image = np.concatenate((image, image[(padded_h - offset):, :, :]), axis=0)
        image = np.concatenate((image[:offset, :, :], image), axis=0)
        image = np.concatenate((image, image[:, (padded_w - offset):, :]), axis=1)
        image = np.concatenate((image[:, :offset, :], image), axis=1)
        image = image * 2.0 - 1.0
        output = np.empty_like(image)

        with torch.inference_mode():
            for i in range(ith):
                for j in range(itw):
                    x = self.stride * i
                    y = self.stride * j
                    tile = image[
                        x:x + self.window_size,
                        y:y + self.window_size,
                        :,
                    ]
                    tile = (
                        torch.from_numpy(tile.transpose(2, 0, 1))
                        .unsqueeze(0)
                        .to(self.device)
                    )
                    result = self.G(tile)
                    result = ((result[0] + 1.0) / 2.0).cpu().numpy().transpose(1, 2, 0)
                    output[
                        x + offset:x + offset + self.stride,
                        y + offset:y + offset + self.stride,
                        :,
                    ] = result[
                        offset:offset + self.stride,
                        offset:offset + self.stride,
                        :,
                    ]
                    progress = int((itw * i + j + 1) * 100 / (ith * itw))
                    print("Transforming input image... %d%%\r" % progress, end="")

        print("Transforming input image... Done!")
        output = np.clip(output, 0.0, 1.0)
        output = output[offset:offset + h, offset:offset + w, :]
        if self.mode == "Greyscale":
            output = output[:, :, 0]

        output_path = Path(out_name)
        output_path.parent.mkdir(parents=True, exist_ok=True)
        self._imwrite(output_path, (output * scale).astype(input_dtype))
        return output_path