Decomposition-Based Variational Network for Multi-Contrast MRI Super-Resolution and Reconstruction (ICCV2023)

import logging
import numpy as np
import cv2
import lmdb
import torch
import torch.utils.data as data
import data.util as util
import os
import random


def roll(x, shift, dim):
    """
    Similar to np.roll but applies to PyTorch Tensors.
    Args:
        x (torch.Tensor): A PyTorch tensor.
        shift (int): Amount to roll.
        dim (int): Which dimension to roll.
    Returns:
        torch.Tensor: Rolled version of x.
    """
    if isinstance(shift, (tuple, list)):
        assert len(shift) == len(dim)
        for s, d in zip(shift, dim):
            x = roll(x, s, d)
        return x
    shift = shift % x.size(dim)
    if shift == 0:
        return x
    left = x.narrow(dim, 0, x.size(dim) - shift)
    right = x.narrow(dim, x.size(dim) - shift, shift)
    return torch.cat((right, left), dim=dim)


def fftshift(x, dim=(-2,-1)):
    """
    Similar to np.fft.fftshift but applies to PyTorch Tensors
    Args:
        x (torch.Tensor): A PyTorch tensor.
        dim (int): Which dimension to fftshift.
    Returns:
        torch.Tensor: fftshifted version of x.
    """
    if dim is None:
        dim = tuple(range(x.dim()))
        shift = [dim // 2 for dim in x.shape]
    elif isinstance(dim, int):
        shift = x.shape[dim] // 2
    else:
        shift = [x.shape[i] // 2 for i in dim]

    return roll(x, shift, dim)

def real_to_complex(img):
    if len(img.shape)==3:
        data = img.unsqueeze(0)
    else:
        data = img
    y = torch.fft.fft2(data)
    y = fftshift(y, dim=(-2,-1))  ## (1,1,h,w)
    y_complex = torch.cat([y.real, y.imag], 1)  ## (1,2,h,w)
    if len(img.shape)==3:
        y_complex = y_complex[0]
    return y_complex

def complex_to_real(data):
    if len(data.shape)==3:
        data1 = data.unsqueeze(0)
    else:
        data1 = data
    h, w = data.shape[-2], data.shape[-1]
    y_real, y_imag = torch.chunk(data1, 2, dim=1)
    y = torch.complex(y_real, y_imag)
    y = fftshift(y, dim=(-2,-1))  ## (1,1,h,w)
    y = torch.fft.irfft2(y,s=(h,w))
    # y = torch.fft.ifft2(y, s=(h, w)).abs()
    if len(data.shape)==3:
        y = y[0]
    return y

def crop_k_data(data, scale):
    _,h,w = data.shape
    lr_h = h//scale
    lr_w = w//scale
    top_left_h = h//2-lr_h//2
    top_left_w = w//2-lr_w//2
    croped_data = data[:, top_left_h:(top_left_h+lr_h), top_left_w:(top_left_w+lr_w)]
    return croped_data

def gen_mask(size, scale):
    h,w = size
    mask = torch.zeros(size)
    lr_h = h//scale
    lr_w = w//scale
    top_left_h = h//2-lr_h//2
    top_left_w = w//2-lr_w//2
    mask[top_left_h:(top_left_h+lr_h), top_left_w:(top_left_w+lr_w)] = torch.ones(lr_h,lr_w)
    return mask

class IXI_train(data.Dataset):
    def __init__(self, opt, train):
        super(IXI_train, self).__init__()
        path = opt['dataroot_GT']
        GT_list = sorted(os.listdir(path))
        self.GT_paths = [os.path.join(path, i) for i in GT_list]
        self.train = train
        if self.train:
            self.crop_size = opt['crop_size']
        else:
            self.crop_size = opt['test_size']
        self.task = opt['task']
        self.scale = int(opt['scale'])
        self.hr_in = opt['hr_in']

    def __len__(self):
        if self.train:
            return len(self.GT_paths)
        else:
            return 2000  # val时只测试200张图片

    def __getitem__(self, idx):

        GT_img_path = self.GT_paths[idx]
        # ref_GT_img_path = self.GT_paths[idx].replace('T2', 'PD')
        ref_GT_img_path = self.GT_paths[idx].replace('T2', 'PD_noalign')
        if self.task == 'rec':
            if self.scale == 4:
                mask_path = '/home/lpc/dataset/IXI/MC_MRI/mask_x4.png'
            elif self.scale == 8:
                mask_path = '/home/lpc/dataset/IXI/MC_MRI/mask_x8.png'
            else:
                print('Wrong scale for reconstruction!')
        elif self.task == 'sr':
            if self.scale == 2:
                mask_path = '/home/lpc/dataset/IXI/MCSR/mask_sr_x2.png'
            elif self.scale == 4:
                mask_path = '/home/lpc/dataset/IXI/MCSR/mask_sr_x4.png'
            else:
                print('Wrong scale for SR!')

        # read image file
        im1_GT = cv2.imread(GT_img_path, cv2.IMREAD_UNCHANGED)
        im2_GT = cv2.imread(ref_GT_img_path, cv2.IMREAD_UNCHANGED)
        mask = cv2.imread(mask_path, cv2.IMREAD_UNCHANGED)

        im1_GT = torch.tensor(im1_GT).unsqueeze(0).float()/255.
        im2_GT = torch.tensor(im2_GT).unsqueeze(0).float()/255. 
        mask = torch.tensor(mask).unsqueeze(0).float().repeat(2,1,1)/255

        # FFT
        im1_GT_k = real_to_complex(im1_GT)
        im2_GT_k = real_to_complex(im2_GT)
        # apply mask
        # zero-padding
        if self.hr_in:
            im1_LQ_k = im1_GT_k * mask
            im2_LQ_k = im2_GT_k * mask
        # center_crop
        else:
            im1_LQ_k = crop_k_data(im1_GT_k, self.scale)
            im2_LQ_k = crop_k_data(im2_GT_k, self.scale)
        # IFFT
        im1_LQ = complex_to_real(im1_LQ_k)
        im2_LQ = complex_to_real(im2_LQ_k)

        # crop
        if self.train:
            _, H, W = im1_GT.shape 
            if self.hr_in:
                rnd_h = random.randint(0, max(0, H - self.crop_size))
                rnd_w = random.randint(0, max(0, W - self.crop_size))
                im1_GT = im1_GT[:, rnd_h:rnd_h + self.crop_size, rnd_w:rnd_w + self.crop_size]
                im2_GT = im2_GT[:, rnd_h:rnd_h + self.crop_size, rnd_w:rnd_w + self.crop_size]
                im1_LQ = im1_LQ[:, rnd_h:rnd_h + self.crop_size, rnd_w:rnd_w + self.crop_size]
                im2_LQ = im2_LQ[:, rnd_h:rnd_h + self.crop_size, rnd_w:rnd_w + self.crop_size]
            else:
                rnd_h = random.randint(0, max(0, H//self.scale - self.crop_size))
                rnd_w = random.randint(0, max(0, W//self.scale - self.crop_size))
                rnd_h_HR, rnd_w_HR = int(rnd_h * self.scale), int(rnd_w * self.scale)
                im1_LQ = im1_LQ[:, rnd_h:rnd_h + self.crop_size, rnd_w:rnd_w + self.crop_size]
                im2_LQ = im2_LQ[:, rnd_h:rnd_h + self.crop_size, rnd_w:rnd_w + self.crop_size]
                im1_GT = im1_GT[:, rnd_h_HR:rnd_h_HR + self.crop_size*self.scale, rnd_w_HR:rnd_w_HR + self.crop_size*self.scale]
                im2_GT = im2_GT[:, rnd_h_HR:rnd_h_HR + self.crop_size*self.scale, rnd_w_HR:rnd_w_HR + self.crop_size*self.scale]
                

        return {'im1_LQ':im1_LQ, 'im1_GT':im1_GT, 'im2_LQ':im2_LQ, 'im2_GT':im2_GT, 'mask':mask}


        

data/IXI_dataset.py
lpcccc-cv/MC-VarNet
'''create dataset and dataloader'''
import logging
import torch
import torch.utils.data


def create_dataloader(dataset, dataset_opt, opt, sampler):
    phase = dataset_opt['phase']
    if phase == 'train':
        if opt['dist']:
            world_size = torch.distributed.get_world_size()
            num_workers = dataset_opt['n_workers']
            assert dataset_opt['batch_size'] % world_size == 0
            batch_size = dataset_opt['batch_size'] // world_size
            shuffle = False
        else:
            num_workers = dataset_opt['n_workers'] * len(opt['gpu_ids'])
            batch_size = dataset_opt['batch_size']
            shuffle = True
        return torch.utils.data.DataLoader(dataset, batch_size=batch_size, shuffle=shuffle,
                                           num_workers=num_workers, sampler=sampler, drop_last=True,
                                           pin_memory=False)
    else:
        return torch.utils.data.DataLoader(dataset, batch_size=8, shuffle=False, num_workers=2,
                                           pin_memory=True)


def create_dataset(dataset_opt, train):
    mode = dataset_opt['mode']
    if mode == 'IXI':
        from data.IXI_dataset import IXI_train as D
    elif mode == 'brain':
        from data.brain_dataset import brain_train as D
    else:
        raise NotImplementedError('Dataset [{:s}] is not recognized.'.format(mode))
    dataset = D(dataset_opt, train)

    return dataset

data/__init__.py
lpcccc-cv/MC-VarNet
import logging
import numpy as np
import cv2
import lmdb
import torch
import torch.utils.data as data
import data.util as util
import os
import random


def roll(x, shift, dim):
    """
    Similar to np.roll but applies to PyTorch Tensors.
    Args:
        x (torch.Tensor): A PyTorch tensor.
        shift (int): Amount to roll.
        dim (int): Which dimension to roll.
    Returns:
        torch.Tensor: Rolled version of x.
    """
    if isinstance(shift, (tuple, list)):
        assert len(shift) == len(dim)
        for s, d in zip(shift, dim):
            x = roll(x, s, d)
        return x
    shift = shift % x.size(dim)
    if shift == 0:
        return x
    left = x.narrow(dim, 0, x.size(dim) - shift)
    right = x.narrow(dim, x.size(dim) - shift, shift)
    return torch.cat((right, left), dim=dim)


def fftshift(x, dim=(-2,-1)):
    """
    Similar to np.fft.fftshift but applies to PyTorch Tensors
    Args:
        x (torch.Tensor): A PyTorch tensor.
        dim (int): Which dimension to fftshift.
    Returns:
        torch.Tensor: fftshifted version of x.
    """
    if dim is None:
        dim = tuple(range(x.dim()))
        shift = [dim // 2 for dim in x.shape]
    elif isinstance(dim, int):
        shift = x.shape[dim] // 2
    else:
        shift = [x.shape[i] // 2 for i in dim]

    return roll(x, shift, dim)

def real_to_complex(img):
    if len(img.shape)==3:
        data = img.unsqueeze(0)
    else:
        data = img
    y = torch.fft.fft2(data)
    y = fftshift(y, dim=(-2,-1))  ## (1,1,h,w)
    y_complex = torch.cat([y.real, y.imag], 1)  ## (1,2,h,w)
    if len(img.shape)==3:
        y_complex = y_complex[0]
    return y_complex

def complex_to_real(data):
    if len(data.shape)==3:
        data1 = data.unsqueeze(0)
    else:
        data1 = data
    h, w = data.shape[-2], data.shape[-1]
    y_real, y_imag = torch.chunk(data1, 2, dim=1)
    y = torch.complex(y_real, y_imag)
    y = fftshift(y, dim=(-2,-1))  ## (1,1,h,w)
    y = torch.fft.irfft2(y, s=(h, w))
    # y = torch.fft.ifft2(y,s=(h,w)).abs()
    if len(data.shape)==3:
        y = y[0]
    return y

def crop_k_data(data, scale):
    _,h,w = data.shape
    lr_h = h//scale
    lr_w = w//scale
    top_left_h = h//2-lr_h//2
    top_left_w = w//2-lr_w//2
    croped_data = data[:, top_left_h:(top_left_h+lr_h), top_left_w:(top_left_w+lr_w)]
    return croped_data

class brain_train(data.Dataset):
    def __init__(self, opt, train):
        super(brain_train, self).__init__()
        path = opt['dataroot_GT']
        GT_list = sorted(os.listdir(path))
        self.GT_paths = [os.path.join(path, i) for i in GT_list]
        self.train = train
        self.crop_size = opt['crop_size']
        self.task = opt['task']
        self.scale = int(opt['scale'])
        self.hr_in = opt['hr_in']

    def __len__(self):
        if self.train:
            return len(self.GT_paths)
        else:
            return len(self.GT_paths)
            # return 200

    def __getitem__(self, idx):

        GT_img_path = self.GT_paths[idx]
        ref_GT_img_path = self.GT_paths[idx].replace('T1', 'T2')
        if self.task == 'rec':
            if self.scale == 4:
                mask_path = '/home/lpc/dataset/IXI/MC_MRI/mask_x4.png'
            elif self.scale == 8:
                mask_path = '/home/lpc/dataset/IXI/MC_MRI/mask_x8.png'
            else:
                print('Wrong scale for reconstruction!')
        elif self.task == 'sr':
            if self.scale == 2:
                mask_path = '/home/lpc/dataset/BrainTS/mask_sr_x2.png'
            elif self.scale == 4:
                mask_path = '/home/lpc/dataset/BrainTS/mask_sr_x4.png'
            else:
                print('Wrong scale for SR!')

        # read image file
        im1_GT = cv2.imread(GT_img_path, cv2.IMREAD_UNCHANGED)
        im2_GT = cv2.imread(ref_GT_img_path, cv2.IMREAD_UNCHANGED)
        mask = cv2.imread(mask_path, cv2.IMREAD_UNCHANGED)

        im1_GT = torch.tensor(im1_GT).unsqueeze(0).float()/255.
        im2_GT = torch.tensor(im2_GT).unsqueeze(0).float()/255. 
        mask = torch.tensor(mask).unsqueeze(0).float().repeat(2,1,1)/255

        # FFT
        im1_GT_k = real_to_complex(im1_GT)
        im2_GT_k = real_to_complex(im2_GT)
        # apply mask
        # zero-padding
        if self.hr_in:
            im1_LQ_k = im1_GT_k * mask
            im2_LQ_k = im2_GT_k * mask
        # center_crop
        else:
            im1_LQ_k = crop_k_data(im1_GT_k, self.scale)
            im2_LQ_k = crop_k_data(im2_GT_k, self.scale)
        # IFFT
        im1_LQ = complex_to_real(im1_LQ_k)
        im2_LQ = complex_to_real(im2_LQ_k)

        # crop
        if self.train:
            _, H, W = im1_GT.shape 
            if self.hr_in:
                rnd_h = random.randint(0, max(0, H - self.crop_size))
                rnd_w = random.randint(0, max(0, W - self.crop_size))
                im1_GT = im1_GT[:, rnd_h:rnd_h + self.crop_size, rnd_w:rnd_w + self.crop_size]
                im2_GT = im2_GT[:, rnd_h:rnd_h + self.crop_size, rnd_w:rnd_w + self.crop_size]
                im1_LQ = im1_LQ[:, rnd_h:rnd_h + self.crop_size, rnd_w:rnd_w + self.crop_size]
                im2_LQ = im2_LQ[:, rnd_h:rnd_h + self.crop_size, rnd_w:rnd_w + self.crop_size]
            else:
                rnd_h = random.randint(0, max(0, H//self.scale - self.crop_size))
                rnd_w = random.randint(0, max(0, W//self.scale - self.crop_size))
                rnd_h_HR, rnd_w_HR = int(rnd_h * self.scale), int(rnd_w * self.scale)
                im1_LQ = im1_LQ[:, rnd_h:rnd_h + self.crop_size, rnd_w:rnd_w + self.crop_size]
                im2_LQ = im2_LQ[:, rnd_h:rnd_h + self.crop_size, rnd_w:rnd_w + self.crop_size]
                im1_GT = im1_GT[:, rnd_h_HR:rnd_h_HR + self.crop_size*self.scale, rnd_w_HR:rnd_w_HR + self.crop_size*self.scale]
                im2_GT = im2_GT[:, rnd_h_HR:rnd_h_HR + self.crop_size*self.scale, rnd_w_HR:rnd_w_HR + self.crop_size*self.scale]
                

        return {'im1_LQ':im1_LQ, 'im1_GT':im1_GT, 'im2_LQ':im2_LQ, 'im2_GT':im2_GT, 'mask':mask}
        

data/brain_dataset.py
lpcccc-cv/MC-VarNet
import os
import pickle
import random
import numpy as np
import cv2
import math
import torch

####################
# Files & IO
####################

###################### get image path list ######################
IMG_EXTENSIONS = ['.jpg', '.JPG', '.jpeg', '.JPEG', '.png', '.PNG', '.ppm', '.PPM', '.bmp', '.BMP']


def is_image_file(filename):
    return any(filename.endswith(extension) for extension in IMG_EXTENSIONS)


def _get_paths_from_images(path):
    '''get image path list from image folder'''
    assert os.path.isdir(path), '{:s} is not a valid directory'.format(path)
    images = []
    for dirpath, _, fnames in sorted(os.walk(path)):
        for fname in sorted(fnames):
            if is_image_file(fname):
                img_path = os.path.join(dirpath, fname)
                images.append(img_path)
    assert images, '{:s} has no valid image file'.format(path)
    return images


def _get_paths_from_lmdb(dataroot):
    '''get image path list from lmdb meta info'''
    meta_info = pickle.load(open(os.path.join(dataroot, 'meta_info.pkl'), 'rb'))
    paths = meta_info['keys']
    sizes = meta_info['resolution']
    if len(sizes) == 1:
        sizes = sizes * len(paths)
    return paths, sizes


def get_image_paths(data_type, dataroot):
    '''get image path list
    support lmdb or image files'''
    paths, sizes = None, None
    if dataroot is not None:
        if data_type == 'lmdb':
            paths, sizes = _get_paths_from_lmdb(dataroot)
        elif data_type == 'img':
            paths = sorted(_get_paths_from_images(dataroot))
        else:
            raise NotImplementedError('data_type [{:s}] is not recognized.'.format(data_type))
    return paths, sizes


###################### read images ######################
def _read_img_lmdb(env, key, size):
    '''read image from lmdb with key (w/ and w/o fixed size)
    size: (C, H, W) tuple'''
    with env.begin(write=False) as txn:
        buf = txn.get(key.encode('ascii'))
    img_flat = np.frombuffer(buf, dtype=np.uint8)
    C, H, W = size
    img = img_flat.reshape(H, W, C)
    return img


def read_img(env, path, size=None):
    '''read image by cv2 or from lmdb
    return: Numpy float32, HWC, BGR, [0,1]'''
    if env is None:  # img
        img = cv2.imread(path, cv2.IMREAD_UNCHANGED)
    else:
        img = _read_img_lmdb(env, path, size)
    img = img.astype(np.float32) / 255.
    if img.ndim == 2:
        img = np.expand_dims(img, axis=2)
    # some images have 4 channels
    if img.shape[2] > 3:
        img = img[:, :, :3]
    return img


####################
# image processing
# process on numpy image
####################


def augment(img_list, hflip=True, rot=True):
    # horizontal flip OR rotate
    hflip = hflip and random.random() < 0.5
    vflip = rot and random.random() < 0.5
    rot90 = rot and random.random() < 0.5

    def _augment(img):
        if hflip:
            img = img[:, ::-1, :]
        if vflip:
            img = img[::-1, :, :]
        if rot90:
            img = img.transpose(1, 0, 2)
        return img

    return [_augment(img) for img in img_list]

def channel_convert(in_c, tar_type, img_list):
    # conversion among BGR, gray and y
    if in_c == 3 and tar_type == 'gray':  # BGR to gray
        gray_list = [cv2.cvtColor(img, cv2.COLOR_BGR2GRAY) for img in img_list]
        return [np.expand_dims(img, axis=2) for img in gray_list]
    elif in_c == 3 and tar_type == 'y':  # BGR to y
        y_list = [bgr2ycbcr(img, only_y=True) for img in img_list]
        return [np.expand_dims(img, axis=2) for img in y_list]
    elif in_c == 1 and tar_type == 'RGB':  # gray/y to BGR
        return [cv2.cvtColor(img, cv2.COLOR_GRAY2BGR) for img in img_list]
    else:
        return img_list


def rgb2ycbcr(img, only_y=True):
    '''same as matlab rgb2ycbcr
    only_y: only return Y channel
    Input:
        uint8, [0, 255]
        float, [0, 1]
    '''
    in_img_type = img.dtype
    img.astype(np.float32)
    if in_img_type != np.uint8:
        img *= 255.
    # convert
    if only_y:
        rlt = np.dot(img, [65.481, 128.553, 24.966]) / 255.0 + 16.0
    else:
        rlt = np.matmul(img, [[65.481, -37.797, 112.0], [128.553, -74.203, -93.786],
                              [24.966, 112.0, -18.214]]) / 255.0 + [16, 128, 128]
    if in_img_type == np.uint8:
        rlt = rlt.round()
    else:
        rlt /= 255.
    return rlt.astype(in_img_type)


def bgr2ycbcr(img, only_y=True):
    '''bgr version of rgb2ycbcr
    only_y: only return Y channel
    Input:
        uint8, [0, 255]
        float, [0, 1]
    '''
    in_img_type = img.dtype
    img.astype(np.float32)
    if in_img_type != np.uint8:
        img *= 255.
    # convert
    if only_y:
        rlt = np.dot(img, [24.966, 128.553, 65.481]) / 255.0 + 16.0
    else:
        rlt = np.matmul(img, [[24.966, 112.0, -18.214], [128.553, -74.203, -93.786],
                              [65.481, -37.797, 112.0]]) / 255.0 + [16, 128, 128]
    if in_img_type == np.uint8:
        rlt = rlt.round()
    else:
        rlt /= 255.
    return rlt.astype(in_img_type)


def ycbcr2rgb(img):
    '''same as matlab ycbcr2rgb
    Input:
        uint8, [0, 255]
        float, [0, 1]
    '''
    in_img_type = img.dtype
    img.astype(np.float32)
    if in_img_type != np.uint8:
        img *= 255.
    # convert
    rlt = np.matmul(img, [[0.00456621, 0.00456621, 0.00456621], [0, -0.00153632, 0.00791071],
                          [0.00625893, -0.00318811, 0]]) * 255.0 + [-222.921, 135.576, -276.836]
    if in_img_type == np.uint8:
        rlt = rlt.round()
    else:
        rlt /= 255.
    return rlt.astype(in_img_type)


def modcrop(img_in, scale):
    # img_in: Numpy, HWC or HW
    img = np.copy(img_in)
    if img.ndim == 2:
        H, W = img.shape
        H_r, W_r = H % scale, W % scale
        img = img[:H - H_r, :W - W_r]
    elif img.ndim == 3:
        H, W, C = img.shape
        H_r, W_r = H % scale, W % scale
        img = img[:H - H_r, :W - W_r, :]
    else:
        raise ValueError('Wrong img ndim: [{:d}].'.format(img.ndim))
    return img

def cubic(x):
    absx = torch.abs(x)
    absx2 = absx**2
    absx3 = absx**3
    return (1.5 * absx3 - 2.5 * absx2 + 1) * (
        (absx <= 1).type_as(absx)) + (-0.5 * absx3 + 2.5 * absx2 - 4 * absx + 2) * ((
            (absx > 1) * (absx <= 2)).type_as(absx))

def calculate_weights_indices(in_length, out_length, scale, kernel, kernel_width, antialiasing):
    if (scale < 1) and (antialiasing):
        # Use a modified kernel to simultaneously interpolate and antialias- larger kernel width
        kernel_width = kernel_width / scale

    # Output-space coordinates
    x = torch.linspace(1, out_length, out_length)

    # Input-space coordinates. Calculate the inverse mapping such that 0.5
    # in output space maps to 0.5 in input space, and 0.5+scale in output
    # space maps to 1.5 in input space.
    u = x / scale + 0.5 * (1 - 1 / scale)

    # What is the left-most pixel that can be involved in the computation?
    left = torch.floor(u - kernel_width / 2)

    # What is the maximum number of pixels that can be involved in the
    # computation?  Note: it's OK to use an extra pixel here; if the
    # corresponding weights are all zero, it will be eliminated at the end
    # of this function.
    P = math.ceil(kernel_width) + 2

    # The indices of the input pixels involved in computing the k-th output
    # pixel are in row k of the indices matrix.
    indices = left.view(out_length, 1).expand(out_length, P) + torch.linspace(0, P - 1, P).view(
        1, P).expand(out_length, P)

    # The weights used to compute the k-th output pixel are in row k of the
    # weights matrix.
    distance_to_center = u.view(out_length, 1).expand(out_length, P) - indices
    # apply cubic kernel
    if (scale < 1) and (antialiasing):
        weights = scale * cubic(distance_to_center * scale)
    else:
        weights = cubic(distance_to_center)
    # Normalize the weights matrix so that each row sums to 1.
    weights_sum = torch.sum(weights, 1).view(out_length, 1)
    weights = weights / weights_sum.expand(out_length, P)

    # If a column in weights is all zero, get rid of it. only consider the first and last column.
    weights_zero_tmp = torch.sum((weights == 0), 0)
    if not math.isclose(weights_zero_tmp[0], 0, rel_tol=1e-6):
        indices = indices.narrow(1, 1, P - 2)
        weights = weights.narrow(1, 1, P - 2)
    if not math.isclose(weights_zero_tmp[-1], 0, rel_tol=1e-6):
        indices = indices.narrow(1, 0, P - 2)
        weights = weights.narrow(1, 0, P - 2)
    weights = weights.contiguous()
    indices = indices.contiguous()
    sym_len_s = -indices.min() + 1
    sym_len_e = indices.max() - in_length
    indices = indices + sym_len_s - 1
    return weights, indices, int(sym_len_s), int(sym_len_e)

def imresize_np(img, scale, antialiasing=True):
    # Now the scale should be the same for H and W
    # input: img: Numpy, HWC BGR [0,1]
    # output: HWC BGR [0,1] w/o round
    img = torch.from_numpy(img)

    in_H, in_W, in_C = img.size()
    _, out_H, out_W = in_C, math.ceil(in_H * scale), math.ceil(in_W * scale)
    kernel_width = 4
    kernel = 'cubic'

    # Return the desired dimension order for performing the resize.  The
    # strategy is to perform the resize first along the dimension with the
    # smallest scale factor.
    # Now we do not support this.

    # get weights and indices
    weights_H, indices_H, sym_len_Hs, sym_len_He = calculate_weights_indices(
        in_H, out_H, scale, kernel, kernel_width, antialiasing)
    weights_W, indices_W, sym_len_Ws, sym_len_We = calculate_weights_indices(
        in_W, out_W, scale, kernel, kernel_width, antialiasing)
    # process H dimension
    # symmetric copying
    img_aug = torch.FloatTensor(in_H + sym_len_Hs + sym_len_He, in_W, in_C)
    img_aug.narrow(0, sym_len_Hs, in_H).copy_(img)

    sym_patch = img[:sym_len_Hs, :, :]
    inv_idx = torch.arange(sym_patch.size(0) - 1, -1, -1).long()
    sym_patch_inv = sym_patch.index_select(0, inv_idx)
    img_aug.narrow(0, 0, sym_len_Hs).copy_(sym_patch_inv)

    sym_patch = img[-sym_len_He:, :, :]
    inv_idx = torch.arange(sym_patch.size(0) - 1, -1, -1).long()
    sym_patch_inv = sym_patch.index_select(0, inv_idx)
    img_aug.narrow(0, sym_len_Hs + in_H, sym_len_He).copy_(sym_patch_inv)

    out_1 = torch.FloatTensor(out_H, in_W, in_C)
    kernel_width = weights_H.size(1)
    for i in range(out_H):
        idx = int(indices_H[i][0])
        out_1[i, :, 0] = img_aug[idx:idx + kernel_width, :, 0].transpose(0, 1).mv(weights_H[i])
        out_1[i, :, 1] = img_aug[idx:idx + kernel_width, :, 1].transpose(0, 1).mv(weights_H[i])
        out_1[i, :, 2] = img_aug[idx:idx + kernel_width, :, 2].transpose(0, 1).mv(weights_H[i])

    # process W dimension
    # symmetric copying
    out_1_aug = torch.FloatTensor(out_H, in_W + sym_len_Ws + sym_len_We, in_C)
    out_1_aug.narrow(1, sym_len_Ws, in_W).copy_(out_1)

    sym_patch = out_1[:, :sym_len_Ws, :]
    inv_idx = torch.arange(sym_patch.size(1) - 1, -1, -1).long()
    sym_patch_inv = sym_patch.index_select(1, inv_idx)
    out_1_aug.narrow(1, 0, sym_len_Ws).copy_(sym_patch_inv)

    sym_patch = out_1[:, -sym_len_We:, :]
    inv_idx = torch.arange(sym_patch.size(1) - 1, -1, -1).long()
    sym_patch_inv = sym_patch.index_select(1, inv_idx)
    out_1_aug.narrow(1, sym_len_Ws + in_W, sym_len_We).copy_(sym_patch_inv)

    out_2 = torch.FloatTensor(out_H, out_W, in_C)
    kernel_width = weights_W.size(1)
    for i in range(out_W):
        idx = int(indices_W[i][0])
        out_2[:, i, 0] = out_1_aug[:, idx:idx + kernel_width, 0].mv(weights_W[i])
        out_2[:, i, 1] = out_1_aug[:, idx:idx + kernel_width, 1].mv(weights_W[i])
        out_2[:, i, 2] = out_1_aug[:, idx:idx + kernel_width, 2].mv(weights_W[i])

    return out_2.numpy()

data/data_util.py
lpcccc-cv/MC-VarNet
import torch
import numpy as np
import contextlib
import cv2

@contextlib.contextmanager
def temp_seed(rng, seed):
    state = rng.get_state()
    rng.seed(seed)
    try:
        yield
    finally:
        rng.set_state(state)
def mask_func_random_unique(shape, acc = 4, seed=42):
    """
    Args:
    shape:[320, 320, 2]

    Return:
    [1, 320, 1]非0即1的tensor
    """
    if len(shape) < 3:
        raise ValueError("Shape should have 3 or more dimensions")
    
    rng = np.random
    with temp_seed(rng, seed):
        num_cols = shape[-2]
        if acc == 4:
            center_fraction, acceleration = 0.08, 2#中心采样比例，加速比
        elif acc == 8:
            center_fraction, acceleration = 0.04, 8#中心采样比例，加速比
        else:
            assert('accelerate rate is not implmented')

        # create the mask
        num_low_freqs = int(round(num_cols * center_fraction))
        #
        #(采样条数-中心采样条数) / 所有未采样条数 计算每一行的采样概率
        #
        prob = (num_cols / acceleration - num_low_freqs) / (
            num_cols - num_low_freqs)
        mask = rng.uniform(size=num_cols) < prob
        pad = (num_cols - num_low_freqs + 1) // 2
        mask[pad: pad + num_low_freqs] = True

        # reshape the mask
        mask_shape = [1 for _ in shape]
        mask_shape[-2] = num_cols
        mask = torch.from_numpy(mask.reshape(*mask_shape).astype(np.float32)) # mask.shape=[col, 1]
        mask = mask.repeat(shape[0], 1, 1)   
    return mask[:, :, 0]

mask = mask_func_random_unique([240, 240, 2])
print(mask.shape)
cv2.imwrite('./mask_x4_brain.png', mask.numpy()*255)
data/generate_mask_random.py
lpcccc-cv/MC-VarNet
"""
Copyright (c) Facebook, Inc. and its affiliates.

This source code is licensed under the MIT license found in the
LICENSE file in the root directory of this source tree.
"""

import torch
import numpy as np

def complex_mul(x, y):
    """
    Complex multiplication.

    This multiplies two complex tensors assuming that they are both stored as
    real arrays with the last dimension being the complex dimension.

    Args:
        x (torch.Tensor): A PyTorch tensor with the last dimension of size 2.
        y (torch.Tensor): A PyTorch tensor with the last dimension of size 2.

    Returns:
        torch.Tensor: A PyTorch tensor with the last dimension of size 2.
    """
    assert x.shape[-1] == y.shape[-1] == 2
    re = x[..., 0] * y[..., 0] - x[..., 1] * y[..., 1]
    im = x[..., 0] * y[..., 1] + x[..., 1] * y[..., 0]

    return torch.stack((re, im), dim=-1)


def complex_conj(x):
    """
    Complex conjugate.

    This applies the complex conjugate assuming that the input array has the
    last dimension as the complex dimension.

    Args:
        x (torch.Tensor): A PyTorch tensor with the last dimension of size 2.
        y (torch.Tensor): A PyTorch tensor with the last dimension of size 2.

    Returns:
        torch.Tensor: A PyTorch tensor with the last dimension of size 2.
    """
    assert x.shape[-1] == 2

    return torch.stack((x[..., 0], -x[..., 1]), dim=-1)


def fft2c(data: torch.Tensor) -> torch.Tensor:
    """
    Apply centered 2 dimensional Fast Fourier Transform.

    Args:
        data: Complex valued input data containing at least 3 dimensions:
        dimensions -3 & -2 are spatial dimensions and dimension -1 has size
        2. All other dimensions are assumed to be batch dimensions.

    Returns:
        The FFT of the input.
    """
    if not data.shape[-1] == 2:
        raise ValueError("Tensor does not have separate complex dim.")

    data = ifftshift(data, dim=[-3, -2])
    data = torch.view_as_real(torch.fft.fftn(torch.view_as_complex(data), dim=(-2, -1), norm="ortho"))
    data = fftshift(data, dim=[-3, -2])

    return data

def ifft2c(data: torch.Tensor) -> torch.Tensor:
    """
    Apply centered 2-dimensional Inverse Fast Fourier Transform.

    Args:
    data: Complex valued input data containing at least 3 dimensions:
        dimensions -3 & -2 are spatial dimensions and dimension -1 has size
        2. All other dimensions are assumed to be batch dimensions.

    Returns:
        The IFFT of the input.
    """
    if not data.shape[-1] == 2:
        raise ValueError("Tensor does not have separate complex dim.")

    data = ifftshift(data, dim=[-3, -2])
    data = torch.view_as_real(torch.fft.ifftn(torch.view_as_complex(data), dim=(-2, -1), norm="ortho"))
    data = fftshift(data, dim=[-3, -2])

    return data


def complex_abs(data):
    """
    Compute the absolute value of a complex valued input tensor.

    Args:
        data (torch.Tensor): A complex valued tensor, where the size of the
            final dimension should be 2.

    Returns:
        torch.Tensor: Absolute value of data.
    """
    assert data.size(-1) == 2

    return (data ** 2).sum(dim=-1).sqrt()



def complex_abs_numpy(data):
    assert data.shape[-1] == 2

    return np.sqrt(np.sum(data ** 2, axis=-1))


def complex_abs_sq(data):#multi coil
    """
    Compute the squared absolute value of a complex tensor.

    Args:
        data (torch.Tensor): A complex valued tensor, where the size of the
            final dimension should be 2.

    Returns:
        torch.Tensor: Squared absolute value of data.
    """
    assert data.size(-1) == 2
    return (data ** 2).sum(dim=-1)


# Helper functions


def roll(x, shift, dim):
    """
    Similar to np.roll but applies to PyTorch Tensors.

    Args:
        x (torch.Tensor): A PyTorch tensor.
        shift (int): Amount to roll.
        dim (int): Which dimension to roll.

    Returns:
        torch.Tensor: Rolled version of x.
    """
    if isinstance(shift, (tuple, list)):
        assert len(shift) == len(dim)
        for s, d in zip(shift, dim):
            x = roll(x, s, d)
        return x
    shift = shift % x.size(dim)
    if shift == 0:
        return x
    left = x.narrow(dim, 0, x.size(dim) - shift)
    right = x.narrow(dim, x.size(dim) - shift, shift)
    return torch.cat((right, left), dim=dim)


def fftshift(x, dim=None):
    """
    Similar to np.fft.fftshift but applies to PyTorch Tensors

    Args:
        x (torch.Tensor): A PyTorch tensor.
        dim (int): Which dimension to fftshift.

    Returns:
        torch.Tensor: fftshifted version of x.
    """
    if dim is None:
        dim = tuple(range(x.dim()))
        shift = [dim // 2 for dim in x.shape]
    elif isinstance(dim, int):
        shift = x.shape[dim] // 2
    else:
        shift = [x.shape[i] // 2 for i in dim]

    return roll(x, shift, dim)


def ifftshift(x, dim=None):
    """
    Similar to np.fft.ifftshift but applies to PyTorch Tensors

    Args:
        x (torch.Tensor): A PyTorch tensor.
        dim (int): Which dimension to ifftshift.

    Returns:
        torch.Tensor: ifftshifted version of x.
    """
    if dim is None:
        dim = tuple(range(x.dim()))
        shift = [(dim + 1) // 2 for dim in x.shape]
    elif isinstance(dim, int):
        shift = (x.shape[dim] + 1) // 2
    else:
        shift = [(x.shape[i] + 1) // 2 for i in dim]

    return roll(x, shift, dim)


def tensor_to_complex_np(data):
    """
    Converts a complex torch tensor to numpy array.
    Args:
        data (torch.Tensor): Input data to be converted to numpy.

    Returns:
        np.array: Complex numpy version of data
    """
    data = data.numpy()
    return data[..., 0] + 1j * data[..., 1]

data/math_lib.py
lpcccc-cv/MC-VarNet
import os
import numpy as np
import cv2
import SimpleITK as sitk

# normalize to 0-255
def norm(data):
    data = data.astype(np.float32)
    # data = np.clip(data, a_min=-200, a_max=400)
    max = np.max(data)
    min = np.min(data)
    data = (data-min)/(max-min)
    return data*255.

# nii folder
file_path = '/home/lpc/dataset/IXI/IXI-PD'
# save image path
save_file_path = '/home/lpc/dataset/IXI/MCSR/test/PD'

step = 1   # sampling interval
skip = 10  # skip low-quality images at the begining and the end.
volume_num = 500  # set the first 500 volumes as the training data
totle_number=0  # total image number

file_list = sorted(os.listdir(file_path))

for number, name in enumerate(file_list):
    if 'IXI014' in name:
        continue  # skip it
    if number<volume_num:
        continue
    print(name)
    filename_1 = file_path+'/'+name
    img_1 = sitk.ReadImage(filename_1, sitk.sitkInt16)
    space = img_1.GetSpacing()
    img_1 = sitk.GetArrayFromImage(img_1)
    width, height, queue = img_1.shape
    # print(name, img.shape)
    # print(name, space)

    # normalize to 0-255
    data_1 = norm(img_1)
    for i in range(skip, width-skip, step):
        totle_number = totle_number+1
        img_arr1 = data_1[i, :, :]
        img_arr1 = np.expand_dims(img_arr1, axis=2)

        cv2.imwrite(save_file_path+'/{:08d}.png'.format(totle_number), img_arr1)
        print('Done!'+save_file_path+'/{:08d}.png'.format(totle_number))


data/read_nii_to_img.py
lpcccc-cv/MC-VarNet
import os
import pickle
import random
import numpy as np
import cv2
import math
import torch

####################
# Files & IO
####################

###################### get image path list ######################
IMG_EXTENSIONS = ['.jpg', '.JPG', '.jpeg', '.JPEG', '.png', '.PNG', '.ppm', '.PPM', '.bmp', '.BMP']


def is_image_file(filename):
    return any(filename.endswith(extension) for extension in IMG_EXTENSIONS)


def _get_paths_from_images(path):
    '''get image path list from image folder'''
    assert os.path.isdir(path), '{:s} is not a valid directory'.format(path)
    images = []
    for dirpath, _, fnames in sorted(os.walk(path)):
        for fname in sorted(fnames):
            if is_image_file(fname):
                img_path = os.path.join(dirpath, fname)
                images.append(img_path)
    assert images, '{:s} has no valid image file'.format(path)
    return images


def _get_paths_from_lmdb(dataroot):
    '''get image path list from lmdb meta info'''
    meta_info = pickle.load(open(os.path.join(dataroot, 'meta_info.pkl'), 'rb'))
    paths = meta_info['keys']
    sizes = meta_info['resolution']
    if len(sizes) == 1:
        sizes = sizes * len(paths)
    return paths, sizes


def get_image_paths(data_type, dataroot):
    '''get image path list
    support lmdb or image files'''
    paths, sizes = None, None
    if dataroot is not None:
        if data_type == 'lmdb':
            paths, sizes = _get_paths_from_lmdb(dataroot)
        elif data_type == 'img':
            paths = sorted(_get_paths_from_images(dataroot))
        else:
            raise NotImplementedError('data_type [{:s}] is not recognized.'.format(data_type))
    return paths, sizes


###################### read images ######################
def _read_img_lmdb(env, key, size):
    '''read image from lmdb with key (w/ and w/o fixed size)
    size: (C, H, W) tuple'''
    with env.begin(write=False) as txn:
        buf = txn.get(key.encode('ascii'))
    img_flat = np.frombuffer(buf, dtype=np.uint8)
    C, H, W = size
    img = img_flat.reshape(H, W, C)
    return img


def read_img(env, path, size=None):
    '''read image by cv2 or from lmdb
    return: Numpy float32, HWC, BGR, [0,1]'''
    if env is None:  # img
        img = cv2.imread(path, cv2.IMREAD_UNCHANGED)
    else:
        img = _read_img_lmdb(env, path, size)
    img = img.astype(np.float32) / 255.
    if img.ndim == 2:
        img = np.expand_dims(img, axis=2)
    return img


####################
# image processing
# process on numpy image
####################


def augment(img_list, hflip=True, rot=True):
    # horizontal flip OR rotate
    hflip = hflip and random.random() < 0.5
    vflip = rot and random.random() < 0.5
    rot90 = rot and random.random() < 0.5

    def _augment(img):
        if hflip:
            img = img[:, ::-1, :]
        if vflip:
            img = img[::-1, :, :]
        if rot90:
            img = img.transpose(1, 0, 2)
        return img

    return [_augment(img) for img in img_list]

def channel_convert(in_c, tar_type, img_list):
    # conversion among BGR, gray and y
    if in_c == 3 and tar_type == 'gray':  # BGR to gray
        gray_list = [cv2.cvtColor(img, cv2.COLOR_BGR2GRAY) for img in img_list]
        return [np.expand_dims(img, axis=2) for img in gray_list]
    elif in_c == 3 and tar_type == 'y':  # BGR to y
        y_list = [bgr2ycbcr(img, only_y=True) for img in img_list]
        return [np.expand_dims(img, axis=2) for img in y_list]
    elif in_c == 1 and tar_type == 'RGB':  # gray/y to BGR
        return [cv2.cvtColor(img, cv2.COLOR_GRAY2BGR) for img in img_list]
    else:
        return img_list


def rgb2ycbcr(img, only_y=True):
    '''same as matlab rgb2ycbcr
    only_y: only return Y channel
    Input:
        uint8, [0, 255]
        float, [0, 1]
    '''
    in_img_type = img.dtype
    img.astype(np.float32)
    if in_img_type != np.uint8:
        img *= 255.
    # convert
    if only_y:
        rlt = np.dot(img, [65.481, 128.553, 24.966]) / 255.0 + 16.0
    else:
        rlt = np.matmul(img, [[65.481, -37.797, 112.0], [128.553, -74.203, -93.786],
                              [24.966, 112.0, -18.214]]) / 255.0 + [16, 128, 128]
    if in_img_type == np.uint8:
        rlt = rlt.round()
    else:
        rlt /= 255.
    return rlt.astype(in_img_type)


def bgr2ycbcr(img, only_y=True):
    '''bgr version of rgb2ycbcr
    only_y: only return Y channel
    Input:
        uint8, [0, 255]
        float, [0, 1]
    '''
    in_img_type = img.dtype
    img.astype(np.float32)
    if in_img_type != np.uint8:
        img *= 255.
    # convert
    if only_y:
        rlt = np.dot(img, [24.966, 128.553, 65.481]) / 255.0 + 16.0
    else:
        rlt = np.matmul(img, [[24.966, 112.0, -18.214], [128.553, -74.203, -93.786],
                              [65.481, -37.797, 112.0]]) / 255.0 + [16, 128, 128]
    if in_img_type == np.uint8:
        rlt = rlt.round()
    else:
        rlt /= 255.
    return rlt.astype(in_img_type)


def ycbcr2rgb(img):
    '''same as matlab ycbcr2rgb
    Input:
        uint8, [0, 255]
        float, [0, 1]
    '''
    in_img_type = img.dtype
    img.astype(np.float32)
    if in_img_type != np.uint8:
        img *= 255.
    # convert
    rlt = np.matmul(img, [[0.00456621, 0.00456621, 0.00456621], [0, -0.00153632, 0.00791071],
                          [0.00625893, -0.00318811, 0]]) * 255.0 + [-222.921, 135.576, -276.836]
    if in_img_type == np.uint8:
        rlt = rlt.round()
    else:
        rlt /= 255.
    return rlt.astype(in_img_type)


def modcrop(img_in, scale):
    # img_in: Numpy, HWC or HW
    img = np.copy(img_in)
    if img.ndim == 2:
        H, W = img.shape
        H_r, W_r = H % scale, W % scale
        img = img[:H - H_r, :W - W_r]
    elif img.ndim == 3:
        H, W, C = img.shape
        H_r, W_r = H % scale, W % scale
        img = img[:H - H_r, :W - W_r, :]
    else:
        raise ValueError('Wrong img ndim: [{:d}].'.format(img.ndim))
    return img

def cubic(x):
    absx = torch.abs(x)
    absx2 = absx**2
    absx3 = absx**3
    return (1.5 * absx3 - 2.5 * absx2 + 1) * (
        (absx <= 1).type_as(absx)) + (-0.5 * absx3 + 2.5 * absx2 - 4 * absx + 2) * ((
            (absx > 1) * (absx <= 2)).type_as(absx))

def calculate_weights_indices(in_length, out_length, scale, kernel, kernel_width, antialiasing):
    if (scale < 1) and (antialiasing):
        # Use a modified kernel to simultaneously interpolate and antialias- larger kernel width
        kernel_width = kernel_width / scale

    # Output-space coordinates
    x = torch.linspace(1, out_length, out_length)

    # Input-space coordinates. Calculate the inverse mapping such that 0.5
    # in output space maps to 0.5 in input space, and 0.5+scale in output
    # space maps to 1.5 in input space.
    u = x / scale + 0.5 * (1 - 1 / scale)

    # What is the left-most pixel that can be involved in the computation?
    left = torch.floor(u - kernel_width / 2)

    # What is the maximum number of pixels that can be involved in the
    # computation?  Note: it's OK to use an extra pixel here; if the
    # corresponding weights are all zero, it will be eliminated at the end
    # of this function.
    P = math.ceil(kernel_width) + 2

    # The indices of the input pixels involved in computing the k-th output
    # pixel are in row k of the indices matrix.
    indices = left.view(out_length, 1).expand(out_length, P) + torch.linspace(0, P - 1, P).view(
        1, P).expand(out_length, P)

    # The weights used to compute the k-th output pixel are in row k of the
    # weights matrix.
    distance_to_center = u.view(out_length, 1).expand(out_length, P) - indices
    # apply cubic kernel
    if (scale < 1) and (antialiasing):
        weights = scale * cubic(distance_to_center * scale)
    else:
        weights = cubic(distance_to_center)
    # Normalize the weights matrix so that each row sums to 1.
    weights_sum = torch.sum(weights, 1).view(out_length, 1)
    weights = weights / weights_sum.expand(out_length, P)

    # If a column in weights is all zero, get rid of it. only consider the first and last column.
    weights_zero_tmp = torch.sum((weights == 0), 0)
    if not math.isclose(weights_zero_tmp[0], 0, rel_tol=1e-6):
        indices = indices.narrow(1, 1, P - 2)
        weights = weights.narrow(1, 1, P - 2)
    if not math.isclose(weights_zero_tmp[-1], 0, rel_tol=1e-6):
        indices = indices.narrow(1, 0, P - 2)
        weights = weights.narrow(1, 0, P - 2)
    weights = weights.contiguous()
    indices = indices.contiguous()
    sym_len_s = -indices.min() + 1
    sym_len_e = indices.max() - in_length
    indices = indices + sym_len_s - 1
    return weights, indices, int(sym_len_s), int(sym_len_e)

def imresize_np(img, scale, antialiasing=True):
    # Now the scale should be the same for H and W
    # input: img: Numpy, HWC BGR [0,1]
    # output: HWC BGR [0,1] w/o round
    img = torch.from_numpy(img)

    in_H, in_W, in_C = img.size()
    _, out_H, out_W = in_C, math.ceil(in_H * scale), math.ceil(in_W * scale)
    kernel_width = 4
    kernel = 'cubic'

    # Return the desired dimension order for performing the resize.  The
    # strategy is to perform the resize first along the dimension with the
    # smallest scale factor.
    # Now we do not support this.

    # get weights and indices
    weights_H, indices_H, sym_len_Hs, sym_len_He = calculate_weights_indices(
        in_H, out_H, scale, kernel, kernel_width, antialiasing)
    weights_W, indices_W, sym_len_Ws, sym_len_We = calculate_weights_indices(
        in_W, out_W, scale, kernel, kernel_width, antialiasing)
    # process H dimension
    # symmetric copying
    img_aug = torch.FloatTensor(in_H + sym_len_Hs + sym_len_He, in_W, in_C)
    img_aug.narrow(0, sym_len_Hs, in_H).copy_(img)

    sym_patch = img[:sym_len_Hs, :, :]
    inv_idx = torch.arange(sym_patch.size(0) - 1, -1, -1).long()
    sym_patch_inv = sym_patch.index_select(0, inv_idx)
    img_aug.narrow(0, 0, sym_len_Hs).copy_(sym_patch_inv)

    sym_patch = img[-sym_len_He:, :, :]
    inv_idx = torch.arange(sym_patch.size(0) - 1, -1, -1).long()
    sym_patch_inv = sym_patch.index_select(0, inv_idx)
    img_aug.narrow(0, sym_len_Hs + in_H, sym_len_He).copy_(sym_patch_inv)

    out_1 = torch.FloatTensor(out_H, in_W, in_C)
    kernel_width = weights_H.size(1)
    for i in range(out_H):
        idx = int(indices_H[i][0])
        out_1[i, :, 0] = img_aug[idx:idx + kernel_width, :, 0].transpose(0, 1).mv(weights_H[i])
        out_1[i, :, 1] = img_aug[idx:idx + kernel_width, :, 1].transpose(0, 1).mv(weights_H[i])
        out_1[i, :, 2] = img_aug[idx:idx + kernel_width, :, 2].transpose(0, 1).mv(weights_H[i])

    # process W dimension
    # symmetric copying
    out_1_aug = torch.FloatTensor(out_H, in_W + sym_len_Ws + sym_len_We, in_C)
    out_1_aug.narrow(1, sym_len_Ws, in_W).copy_(out_1)

    sym_patch = out_1[:, :sym_len_Ws, :]
    inv_idx = torch.arange(sym_patch.size(1) - 1, -1, -1).long()
    sym_patch_inv = sym_patch.index_select(1, inv_idx)
    out_1_aug.narrow(1, 0, sym_len_Ws).copy_(sym_patch_inv)

    sym_patch = out_1[:, -sym_len_We:, :]
    inv_idx = torch.arange(sym_patch.size(1) - 1, -1, -1).long()
    sym_patch_inv = sym_patch.index_select(1, inv_idx)
    out_1_aug.narrow(1, sym_len_Ws + in_W, sym_len_We).copy_(sym_patch_inv)

    out_2 = torch.FloatTensor(out_H, out_W, in_C)
    kernel_width = weights_W.size(1)
    for i in range(out_W):
        idx = int(indices_W[i][0])
        out_2[:, i, 0] = out_1_aug[:, idx:idx + kernel_width, 0].mv(weights_W[i])
        out_2[:, i, 1] = out_1_aug[:, idx:idx + kernel_width, 1].mv(weights_W[i])
        out_2[:, i, 2] = out_1_aug[:, idx:idx + kernel_width, 2].mv(weights_W[i])

    return out_2.numpy()



data/util.py
lpcccc-cv/MC-VarNet
import logging
logger = logging.getLogger('base')


def create_model(opt):
    model = opt['model']
    if model == 'joint-rec':
        from .joint_rec_base_model import BaseModel as M
    elif model == 'ref-rec':
        from .ref_rec_base_model import BaseModel as M
    else:
        raise NotImplementedError('Model [{:s}] not recognized.'.format(model))
    m = M(opt)
    logger.info('Model [{:s}] is created.'.format(m.__class__.__name__))
    return m

models/__init__.py
lpcccc-cv/MC-VarNet
import logging
from collections import OrderedDict
import torch.nn.functional as F
import torch
import torch.nn as nn
from torch.nn.parallel import DataParallel, DistributedDataParallel
import models.networks as networks
import models.lr_scheduler as lr_scheduler
from .base_model import BaseModel
from models.modules.loss import CharbonnierLoss, LapLoss, ncc_loss
from torchvision import models
import cv2
from matplotlib import pyplot as plt
import os

logger = logging.getLogger('base')


def vis_img(img, fname, ftype ,output_dir):
    os.makedirs(output_dir, exist_ok=True)
    plt.figure()
    plt.imshow(img, cmap='gray')
    figname = fname + '_' + ftype + '.png'
    figpath = os.path.join(output_dir, figname)
    plt.savefig(figpath)

class BaseModel(BaseModel):
    def __init__(self, opt):
        super(BaseModel, self).__init__(opt)

        if opt['dist']:
            self.rank = torch.distributed.get_rank()
        else:
            self.rank = -1  # non dist training
        train_opt = opt['train']
        self.dwt = opt["dwt"]
        self.which_dataset = opt["mode"]

        # define network and load pretrained models
        self.netG = networks.define_G(opt).to(self.device)
        
        if opt['dist']:
            self.netG = DistributedDataParallel(self.netG, device_ids=[torch.cuda.current_device()], find_unused_parameters=True)
        else:
            self.netG = DataParallel(self.netG)
        # print network
        self.print_network()
        self.load()

        if self.is_train:
            self.netG.train()

            #### loss
            loss_type = train_opt['pixel_criterion']
            if loss_type == 'l1':
                self.cri_pix = nn.L1Loss(reduction='mean').to(self.device)
            elif loss_type == 'l2':
                self.cri_pix = nn.MSELoss(reduction='sum').to(self.device)
            elif loss_type == 'cb':
                self.cri_pix = CharbonnierLoss().to(self.device)
            elif loss_type == 'lp':
                self.cri_pix = LapLoss(max_levels=5).to(self.device)

                raise NotImplementedError('Loss type [{:s}] is not recognized.'.format(loss_type))
            self.l_pix_w = train_opt['pixel_weight']
            self.l1_loss = nn.L1Loss(reduction='mean').to(self.device)
            self.ncc_loss = ncc_loss().to(self.device)


            #### optimizers
            wd_G = train_opt['weight_decay_G'] if train_opt['weight_decay_G'] else 0
            optim_params = []
            for k, v in self.netG.named_parameters():
                if opt['datasets']['train']['use_time'] == True and 'TMB' not in k:
                    v.requires_grad = False
                if v.requires_grad:
                    optim_params.append(v)
                else:
                    if self.rank <= 0:
                        logger.warning('Params [{:s}] will not optimize.'.format(k))

            self.optimizer_G = torch.optim.Adam(optim_params, lr=train_opt['lr_G'],
                                                weight_decay=wd_G,
                                                betas=(train_opt['beta1'], train_opt['beta2']))
            # self.optimizer_G = torch.optim.SGD(optim_params, lr=train_opt['lr_G'], momentum=0.9, weight_decay=wd_G)
            self.optimizers.append(self.optimizer_G)
            #### schedulers
            if train_opt['lr_scheme'] == 'MultiStepLR':
                for optimizer in self.optimizers:
                    self.schedulers.append(
                        lr_scheduler.MultiStepLR_Restart(optimizer, train_opt['lr_steps'],
                                                         restarts=train_opt['restarts'],
                                                         weights=train_opt['restart_weights'],
                                                         gamma=train_opt['lr_gamma'],
                                                         clear_state=train_opt['clear_state']))
            elif train_opt['lr_scheme'] == 'CosineAnnealingLR_Restart':
                for optimizer in self.optimizers:
                    self.schedulers.append(
                        lr_scheduler.CosineAnnealingLR_Restart(
                            optimizer, train_opt['T_period'], eta_min=train_opt['eta_min'],
                            restarts=train_opt['restarts'], weights=train_opt['restart_weights']))
            else:
                raise NotImplementedError()

            self.log_dict = OrderedDict()


    def feed_data(self, data):

        self.im1_lr = data['im1_LQ'].to(self.device)
        self.im1_gt = data['im1_GT'].to(self.device)
        self.im2_lr = data['im2_LQ'].to(self.device)
        self.im2_gt = data['im2_GT'].to(self.device)
        self.mask = data['mask'].to(self.device)


    def set_params_lr_zero(self):
        # fix normal module
        self.optimizers[0].param_groups[0]['lr'] = 0

    def optimize_parameters(self):
        self.optimizer_G.zero_grad()
        self.fake_H_1, self.fake_H_2 = self.netG(self.im1_lr, self.im2_lr, self.mask)
        
        # L1 loss
        l_pix_1 = self.cri_pix(self.fake_H_1, self.im1_gt)
        l_pix_2 = self.cri_pix(self.fake_H_2, self.im2_gt)

        total_loss = l_pix_1 + l_pix_2
        total_loss.backward()
        self.optimizer_G.step()

        # set log
        self.log_dict['l_pix'] = l_pix_1.item()
        self.log_dict['l_ref'] = l_pix_2.item()

    def test(self):
        self.netG.eval()
        with torch.no_grad():
            self.fake_H_1, self.fake_H_2 = self.netG(self.im1_lr, self.im2_lr, self.mask)
        self.netG.train()

    # def optimize_parameters(self):
    #     self.optimizer_G.zero_grad()
    #     self.fake_H = self.netG(self.lr)  
    #     # L1 loss
    #     l_pix = self.cri_pix(self.fake_H, self.gt)
    #     total_loss = l_pix 
    #     total_loss.backward()
    #     self.optimizer_G.step()
    #     # set log
    #     self.log_dict['l_pix'] = l_pix.item()

    # def test(self):
    #     self.netG.eval()
    #     with torch.no_grad():
    #         self.fake_H = self.netG(self.lr)
    #     self.netG.train()

    def get_current_log(self):
        return self.log_dict

    
    def get_current_visuals(self, need_GT=True):

        out_dict = OrderedDict()
        out_dict['im1_LQ'] = self.im1_lr.detach().float().cpu()
        out_dict['im2_LQ'] = self.im2_lr.detach().float().cpu()
        # if self.which_dataset == 'IXI':
        out_dict['im1_restore'] = self.fake_H_1.detach().float().cpu()
        out_dict['im1_GT'] = self.im1_gt.detach().float().cpu() 
        out_dict['im2_restore'] = self.fake_H_2.detach().float().cpu()
        out_dict['im2_GT'] = self.im2_gt.detach().float().cpu()   
  
        return out_dict

    def print_network(self):
        s, n = self.get_network_description(self.netG)
        if isinstance(self.netG, nn.DataParallel):
            net_struc_str = '{} - {}'.format(self.netG.__class__.__name__,
                                             self.netG.module.__class__.__name__)
        else:
            net_struc_str = '{}'.format(self.netG.__class__.__name__)
        if self.rank <= 0:
            logger.info('Network G structure: {}, with parameters: {:,d}'.format(net_struc_str, n))
            logger.info(s)

    def load(self):
        load_path_G = self.opt['path']['pretrain_model_G']
        if load_path_G is not None:
            logger.info('Loading model for G [{:s}] ...'.format(load_path_G))
            self.load_network(load_path_G, self.netG, self.opt['path']['strict_load'])

    def save(self, iter_label, epoch=None):
        if epoch != None:
            self.save_network(self.netG, 'G', iter_label, str(epoch))
        else:
            self.save_network(self.netG, 'G', iter_label)

models/joint_rec_base_model.py
lpcccc-cv/MC-VarNet
import torch
import torch.nn as nn
import torch.nn.functional as F
import functools
from .module_util import *
from data.IXI_dataset import complex_to_real, real_to_complex
import functools
from .module_util import *

class Encoding_Block(torch.nn.Module):
    def __init__(self, c_in, n_feat):
        super(Encoding_Block, self).__init__()

        self.conv1 = torch.nn.Conv2d(in_channels=c_in, out_channels=n_feat, kernel_size=3, padding=3 // 2)
        self.conv2 = torch.nn.Conv2d(in_channels=n_feat, out_channels=n_feat, kernel_size=3, padding=3 // 2)
        self.conv3 = torch.nn.Conv2d(in_channels=n_feat, out_channels=n_feat, kernel_size=3, padding=3 // 2)
        self.conv4 = torch.nn.Conv2d(in_channels=n_feat, out_channels=n_feat, kernel_size=3, stride = 2, padding=3 // 2)
        self.act =  torch.nn.PReLU()

        self.reset_parameters()

    def reset_parameters(self):
        for m in self.modules():
            if isinstance(m, torch.nn.Conv2d):
                torch.nn.init.xavier_normal_(m.weight)

    def forward(self, input):

        f_e = self.conv3(self.act(self.conv2(self.act(self.conv1(input)))))
        down = self.conv4(f_e)
        return f_e, down



class Decoding_Block(torch.nn.Module):
    def __init__(self,c_in, n_feat):
        super(Decoding_Block, self).__init__()
        self.conv1 = torch.nn.Conv2d(in_channels=n_feat*2, out_channels=n_feat, kernel_size=3, padding=3 // 2)
        self.conv2 = torch.nn.Conv2d(in_channels=n_feat, out_channels=n_feat, kernel_size=3, padding=3 // 2)
        self.conv3 = torch.nn.Conv2d(in_channels=n_feat, out_channels=n_feat, kernel_size=3, padding=3 // 2)
        self.conv4 = torch.nn.Conv2d(in_channels=n_feat, out_channels=n_feat, kernel_size=3, padding=3 // 2)
        self.up = torch.nn.ConvTranspose2d(c_in, n_feat, kernel_size=3, stride=2,padding=3 // 2)

        self.act =  torch.nn.PReLU()

        self.reset_parameters()

    def reset_parameters(self):
        for m in self.modules():
            if isinstance(m, torch.nn.Conv2d):
                torch.nn.init.xavier_normal_(m.weight)

    def forward(self, input, map):

        up = self.up(input, output_size = [input.shape[0],input.shape[1],map.shape[2],map.shape[3]])
        out1 = self.act(self.conv1(torch.cat((up, map), 1)))
        out = self.conv4(self.act(self.conv3(self.act(self.conv2(out1)))))

        return out

class Feature_Decoding_End(torch.nn.Module):
    def __init__(self, n_feat, c_out):
        super(Feature_Decoding_End, self).__init__()
        self.conv1 = torch.nn.Conv2d(in_channels=n_feat, out_channels=n_feat, kernel_size=3, padding=3 // 2)
        self.conv2 = torch.nn.Conv2d(in_channels=n_feat, out_channels=n_feat, kernel_size=3, padding=3 // 2)
        self.conv3 = torch.nn.Conv2d(in_channels=n_feat, out_channels=n_feat, kernel_size=3, padding=3 // 2)
        self.conv_out = torch.nn.Conv2d(in_channels=n_feat, out_channels=c_out, kernel_size=3, padding=3 // 2)
        
        self.act =  torch.nn.PReLU()
        
        self.reset_parameters()

    def reset_parameters(self):
        for m in self.modules():
            if isinstance(m, torch.nn.Conv2d):
                torch.nn.init.xavier_normal_(m.weight)

    def forward(self, input):
        out = self.conv3(self.act(self.conv2(self.act(self.conv1(input)))))
        out = self.conv_out(out)
        return out

class Unet_denoise(torch.nn.Module):
    def __init__(self, cin, n_feat):
        super(Unet_denoise, self).__init__()

        self.Encoding_block1 = Encoding_Block(cin, n_feat)
        self.Encoding_block2 = Encoding_Block(n_feat, n_feat)
        self.Encoding_block3 = Encoding_Block(n_feat, n_feat)

        self.Decoding_block1 = Decoding_Block(n_feat, n_feat)
        self.Decoding_block2 = Decoding_Block(n_feat, n_feat)
        self.Decoding_block_End = Feature_Decoding_End(n_feat, cin)

        self.acti = torch.nn.PReLU()
        self.reset_parameters()
    def reset_parameters(self):
        for m in self.modules():
            if isinstance(m, torch.nn.Conv2d):
                torch.nn.init.xavier_normal_(m.weight)

    def forward(self, x):

        encode0, down0 = self.Encoding_block1(x)
        encode1, down1 = self.Encoding_block2(down0)
        encode2, _ = self.Encoding_block3(down1)

        decode1 = self.Decoding_block1(encode2, encode1)
        decode0 = self.Decoding_block2(decode1, encode0)
        decode0 = self.Decoding_block_End(decode0)
        return decode0+x



class VANet(nn.Module):
    def __init__(self):
        super(VANet, self).__init__()

        self.in_channel = 32
        self.channel_fea = 32
        self.iter_num = 4

        ## map to image domain
        self.map_rec_x = nn.Conv2d(in_channels=self.in_channel, out_channels=1, kernel_size=3, padding=3 // 2)
        self.map_rec_y = nn.Conv2d(in_channels=self.in_channel, out_channels=1, kernel_size=3, padding=3 // 2)

        self.init_x = Unet_denoise(cin=self.in_channel, n_feat=self.channel_fea)
        self.init_c = Unet_denoise(cin=self.in_channel, n_feat=self.channel_fea)

        ## denoise networks
        self.denoise_x = Unet_denoise(cin=self.in_channel, n_feat=self.channel_fea)
        self.denoise_c = Unet_denoise(cin=self.in_channel, n_feat=self.channel_fea)
        self.denoise_u = Unet_denoise(cin=self.in_channel, n_feat=self.channel_fea)

        ## CDic layer
        self.E_x = nn.Sequential(*[nn.Conv2d(in_channels=self.in_channel, out_channels=self.channel_fea, kernel_size=3, padding=3 // 2), \
                                    nn.ReLU(), \
                                   nn.Conv2d(in_channels=self.channel_fea, out_channels=self.channel_fea, kernel_size=3, padding=3 // 2)])
        self.D_x = nn.Sequential(*[nn.Conv2d(in_channels=self.channel_fea, out_channels=self.channel_fea, kernel_size=3, padding=3 // 2), \
                                    nn.ReLU(), \
                                   nn.Conv2d(in_channels=self.channel_fea, out_channels=self.channel_fea, kernel_size=3, padding=3 // 2)])
        self.E_c = nn.Sequential(*[nn.Conv2d(in_channels=self.in_channel, out_channels=self.channel_fea, kernel_size=3, padding=3 // 2), \
                                    nn.ReLU(), \
                                   nn.Conv2d(in_channels=self.channel_fea, out_channels=self.channel_fea, kernel_size=3, padding=3 // 2)])
        self.D_c = nn.Sequential(*[nn.Conv2d(in_channels=self.channel_fea, out_channels=self.channel_fea, kernel_size=3, padding=3 // 2), \
                                    nn.ReLU(), \
                                   nn.Conv2d(in_channels=self.channel_fea, out_channels=self.channel_fea, kernel_size=3, padding=3 // 2)])
        self.T_c = nn.Sequential(*[nn.Conv2d(in_channels=self.in_channel, out_channels=self.in_channel, kernel_size=3, padding=3 // 2), \
                                    nn.ReLU(), \
                                   nn.Conv2d(in_channels=self.in_channel, out_channels=self.in_channel, kernel_size=3, padding=3 // 2)])
        
        self.mu_x = [nn.Parameter(torch.tensor(0.1)) for _ in range(self.iter_num)]
        self.mu_c = [nn.Parameter(torch.tensor(0.1)) for _ in range(self.iter_num)]

        self.lamda_1 = [nn.Parameter(torch.tensor(1.0)) for _ in range(self.iter_num)]
        self.gamma_x = [nn.Parameter(torch.tensor(1.0)) for _ in range(self.iter_num)]
        self.gamma_c = [nn.Parameter(torch.tensor(1.0)) for _ in range(self.iter_num)]
        self.alpha = [nn.Parameter(torch.tensor(1.0)) for _ in range(self.iter_num)]
        self.beta = [nn.Parameter(torch.tensor(1.0)) for _ in range(self.iter_num)]

    
    def data_consistency_layer(self, generated, X_k, mask):

        gene_complex = real_to_complex(generated)
        output_complex = X_k + gene_complex * (1.0 - mask)
        output_img = complex_to_real(output_complex)
        
        return output_img
    
    def forward(self, x, y, mask):
        # x: input low-quality target image [b,1,h,w]
        # y: input high-quality ref image [b,1,h,w]
        # mask: undersampling mask [b,1,h,w]

        # repeat input image to high-dimension
        x = x.repeat(1,self.in_channel,1,1)
        y = y.repeat(1,self.in_channel,1,1)

        # x_k
        mask = mask.repeat(1,self.in_channel,1,1)
        x_k = real_to_complex(x)

        # initialize variables
        x_t0 = self.init_x(x)
        c_t0 = self.init_c(y)
        u_t0 = y-c_t0
        P_t0 = c_t0
        Q_t0 = u_t0

        for i in range(self.iter_num):
            # update x
            DC_x = self.data_consistency_layer(x_t0, x_k, mask)  ## DC
            refine = self.gamma_x[i]*self.D_x(self.E_x(x_t0)-self.E_c(c_t0)) + self.lamda_1[i]*self.denoise_x(x_t0)  # refine
            x_t1 = DC_x - self.mu_x[i]*refine # update x
            if i != self.iter_num-1:
                # update C and U
                c_t1 = self.T_c(y - u_t0 + self.gamma_c[i]*self.D_c(self.E_x(x_t0)) + self.alpha[i]*P_t0)
                u_t1 = (y + self.beta[i]*Q_t0 - c_t0)/(1+self.beta[i])
                # update P and Q
                P_t1 = self.denoise_c(c_t1)
                Q_t1 = self.denoise_u(u_t1)

                c_t0 = c_t1
                u_t0 = u_t1
                P_t0 = P_t1
                Q_t0 = Q_t1
            x_t0 = x_t1
        
        # map feature to image
        x_out = self.map_rec_x(x_t0)
        # y_out = self.map_rec_y(c_t0+u_t0)

        return x_out
        # return x_out, y_out

models/modules/VANet.py
lpcccc-cv/MC-VarNet
import math
import torch
from torch import nn as nn
from torch.nn import functional as F
from torch.nn import init as init
from torch.nn.modules.batchnorm import _BatchNorm

# from basicsr.utils import get_root_logger

try:
    from basicsr.models.ops.dcn import (ModulatedDeformConvPack,
                                        modulated_deform_conv)
except ImportError:
    print('Cannot import dcn. Ignore this warning if dcn is not used. '
          'Otherwise install BasicSR with compiling dcn.')
    ModulatedDeformConvPack = object
    modulated_deform_conv = None


@torch.no_grad()
def default_init_weights(module_list, scale=1, bias_fill=0, **kwargs):
    """Initialize network weights.

    Args:
        module_list (list[nn.Module] | nn.Module): Modules to be initialized.
        scale (float): Scale initialized weights, especially for residual
            blocks. Default: 1.
        bias_fill (float): The value to fill bias. Default: 0
        kwargs (dict): Other arguments for initialization function.
    """
    if not isinstance(module_list, list):
        module_list = [module_list]
    for module in module_list:
        for m in module.modules():
            if isinstance(m, nn.Conv2d):
                init.kaiming_normal_(m.weight, **kwargs)
                m.weight.data *= scale
                if m.bias is not None:
                    m.bias.data.fill_(bias_fill)
            elif isinstance(m, nn.Linear):
                init.kaiming_normal_(m.weight, **kwargs)
                m.weight.data *= scale
                if m.bias is not None:
                    m.bias.data.fill_(bias_fill)
            elif isinstance(m, _BatchNorm):
                init.constant_(m.weight, 1)
                if m.bias is not None:
                    m.bias.data.fill_(bias_fill)


def make_layer(basic_block, num_basic_block, **kwarg):
    """Make layers by stacking the same blocks.

    Args:
        basic_block (nn.module): nn.module class for basic block.
        num_basic_block (int): number of blocks.

    Returns:
        nn.Sequential: Stacked blocks in nn.Sequential.
    """
    layers = []
    for _ in range(num_basic_block):
        layers.append(basic_block(**kwarg))
    return nn.Sequential(*layers)


class ResidualBlockNoBN(nn.Module):
    """Residual block without BN.

    It has a style of:
        ---Conv-ReLU-Conv-+-
         |________________|

    Args:
        num_feat (int): Channel number of intermediate features.
            Default: 64.
        res_scale (float): Residual scale. Default: 1.
        pytorch_init (bool): If set to True, use pytorch default init,
            otherwise, use default_init_weights. Default: False.
    """

    def __init__(self, num_feat=64, res_scale=1, pytorch_init=False):
        super(ResidualBlockNoBN, self).__init__()
        self.res_scale = res_scale
        self.conv1 = nn.Conv2d(num_feat, num_feat, 3, 1, 1, bias=True)
        self.conv2 = nn.Conv2d(num_feat, num_feat, 3, 1, 1, bias=True)
        self.relu = nn.ReLU(inplace=True)

        if not pytorch_init:
            default_init_weights([self.conv1, self.conv2], 0.1)

    def forward(self, x):
        identity = x
        out = self.conv2(self.relu(self.conv1(x)))
        return identity + out * self.res_scale


class Upsample(nn.Sequential):
    """Upsample module.

    Args:
        scale (int): Scale factor. Supported scales: 2^n and 3.
        num_feat (int): Channel number of intermediate features.
    """

    def __init__(self, scale, num_feat):
        m = []
        if (scale & (scale - 1)) == 0:  # scale = 2^n
            for _ in range(int(math.log(scale, 2))):
                m.append(nn.Conv2d(num_feat, 4 * num_feat, 3, 1, 1))
                m.append(nn.PixelShuffle(2))
        elif scale == 3:
            m.append(nn.Conv2d(num_feat, 9 * num_feat, 3, 1, 1))
            m.append(nn.PixelShuffle(3))
        else:
            raise ValueError(f'scale {scale} is not supported. '
                             'Supported scales: 2^n and 3.')
        super(Upsample, self).__init__(*m)


def flow_warp(x,
              flow,
              interp_mode='bilinear',
              padding_mode='zeros',
              align_corners=True):
    """Warp an image or feature map with optical flow.

    Args:
        x (Tensor): Tensor with size (n, c, h, w).
        flow (Tensor): Tensor with size (n, h, w, 2), normal value.
        interp_mode (str): 'nearest' or 'bilinear'. Default: 'bilinear'.
        padding_mode (str): 'zeros' or 'border' or 'reflection'.
            Default: 'zeros'.
        align_corners (bool): Before pytorch 1.3, the default value is
            align_corners=True. After pytorch 1.3, the default value is
            align_corners=False. Here, we use the True as default.

    Returns:
        Tensor: Warped image or feature map.
    """
    assert x.size()[-2:] == flow.size()[1:3]
    _, _, h, w = x.size()
    # create mesh grid
    grid_y, grid_x = torch.meshgrid(
        torch.arange(0, h).type_as(x),
        torch.arange(0, w).type_as(x))
    grid = torch.stack((grid_x, grid_y), 2).float()  # W(x), H(y), 2
    grid.requires_grad = False

    vgrid = grid + flow
    # scale grid to [-1,1]
    vgrid_x = 2.0 * vgrid[:, :, :, 0] / max(w - 1, 1) - 1.0
    vgrid_y = 2.0 * vgrid[:, :, :, 1] / max(h - 1, 1) - 1.0
    vgrid_scaled = torch.stack((vgrid_x, vgrid_y), dim=3)
    output = F.grid_sample(
        x,
        vgrid_scaled,
        mode=interp_mode,
        padding_mode=padding_mode,
        align_corners=align_corners)

    # TODO, what if align_corners=False
    return output


def resize_flow(flow,
                size_type,
                sizes,
                interp_mode='bilinear',
                align_corners=False):
    """Resize a flow according to ratio or shape.

    Args:
        flow (Tensor): Precomputed flow. shape [N, 2, H, W].
        size_type (str): 'ratio' or 'shape'.
        sizes (list[int | float]): the ratio for resizing or the final output
            shape.
            1) The order of ratio should be [ratio_h, ratio_w]. For
            downsampling, the ratio should be smaller than 1.0 (i.e., ratio
            < 1.0). For upsampling, the ratio should be larger than 1.0 (i.e.,
            ratio > 1.0).
            2) The order of output_size should be [out_h, out_w].
        interp_mode (str): The mode of interpolation for resizing.
            Default: 'bilinear'.
        align_corners (bool): Whether align corners. Default: False.

    Returns:
        Tensor: Resized flow.
    """
    _, _, flow_h, flow_w = flow.size()
    if size_type == 'ratio':
        output_h, output_w = int(flow_h * sizes[0]), int(flow_w * sizes[1])
    elif size_type == 'shape':
        output_h, output_w = sizes[0], sizes[1]
    else:
        raise ValueError(
            f'Size type should be ratio or shape, but got type {size_type}.')

    input_flow = flow.clone()
    ratio_h = output_h / flow_h
    ratio_w = output_w / flow_w
    input_flow[:, 0, :, :] *= ratio_w
    input_flow[:, 1, :, :] *= ratio_h
    resized_flow = F.interpolate(
        input=input_flow,
        size=(output_h, output_w),
        mode=interp_mode,
        align_corners=align_corners)
    return resized_flow


# TODO: may write a cpp file
def pixel_unshuffle(x, scale):
    """ Pixel unshuffle.

    Args:
        x (Tensor): Input feature with shape (b, c, hh, hw).
        scale (int): Downsample ratio.

    Returns:
        Tensor: the pixel unshuffled feature.
    """
    b, c, hh, hw = x.size()
    out_channel = c * (scale**2)
    assert hh % scale == 0 and hw % scale == 0
    h = hh // scale
    w = hw // scale
    x_view = x.view(b, c, h, scale, w, scale)
    return x_view.permute(0, 1, 3, 5, 2, 4).reshape(b, out_channel, h, w)


class DCNv2Pack(ModulatedDeformConvPack):
    """Modulated deformable conv for deformable alignment.

    Different from the official DCNv2Pack, which generates offsets and masks
    from the preceding features, this DCNv2Pack takes another different
    features to generate offsets and masks.

    Ref:
        Delving Deep into Deformable Alignment in Video Super-Resolution.
    """

    def forward(self, x, feat):
        out = self.conv_offset(feat)
        o1, o2, mask = torch.chunk(out, 3, dim=1)
        offset = torch.cat((o1, o2), dim=1)
        mask = torch.sigmoid(mask)

        offset_absmean = torch.mean(torch.abs(offset))
        if offset_absmean > 50:
            print('Offset abs mean is larger than 50.')
            # logger = get_root_logger()
            # logger.warning(
            #     f'Offset abs mean is {offset_absmean}, larger than 50.')

        return modulated_deform_conv(x, offset, mask, self.weight, self.bias,
                                     self.stride, self.padding, self.dilation,
                                     self.groups, self.deformable_groups)

models/modules/arch_util.py
lpcccc-cv/MC-VarNet
import math

import torch
import torch.nn as nn
import torch.nn.functional as F

from torch.autograd import Variable

def default_conv(in_channels, out_channels, kernel_size, bias=True):
    return nn.Conv2d(
        in_channels, out_channels, kernel_size,
        padding=(kernel_size//2), bias=bias)

class MeanShift(nn.Conv2d):
    def __init__(self, rgb_range, rgb_mean, rgb_std, sign=-1):
        super(MeanShift, self).__init__(3, 3, kernel_size=1)
        std = torch.Tensor(rgb_std)
        self.weight.data = torch.eye(3).view(3, 3, 1, 1)
        self.weight.data.div_(std.view(3, 1, 1, 1))
        self.bias.data = sign * rgb_range * torch.Tensor(rgb_mean)
        self.bias.data.div_(std)
        self.requires_grad = False

class BasicBlock(nn.Sequential):
    def __init__(
        self, in_channels, out_channels, kernel_size, stride=1, bias=False,
        bn=True, act=nn.ReLU(True)):

        m = [nn.Conv2d(
            in_channels, out_channels, kernel_size,
            padding=(kernel_size//2), stride=stride, bias=bias)
        ]
        if bn: m.append(nn.BatchNorm2d(out_channels))
        if act is not None: m.append(act)
        super(BasicBlock, self).__init__(*m)

class ResBlock(nn.Module):
    def __init__(
        self, conv, n_feat, kernel_size,
        bias=True, bn=False, act=nn.ReLU(True), res_scale=1):

        super(ResBlock, self).__init__()
        m = []
        for i in range(2):
            m.append(conv(n_feat, n_feat, kernel_size, bias=bias))
            if bn: m.append(nn.BatchNorm2d(n_feat))
            if i == 0: m.append(act)

        self.body = nn.Sequential(*m)
        self.res_scale = res_scale

    def forward(self, x):
        res = self.body(x).mul(self.res_scale)
        res += x

        return res

class DenseBlock(nn.Module):
    def __init__(self, channel_in, channel_out, init='xavier', gc=32, bias=True):
        super(DenseBlock, self).__init__()
        self.conv1 = nn.Conv2d(channel_in, gc, 3, 1, 1, bias=bias)
        self.conv2 = nn.Conv2d(channel_in + gc, gc, 3, 1, 1, bias=bias)
        self.conv3 = nn.Conv2d(channel_in + 2 * gc, gc, 3, 1, 1, bias=bias)
        self.conv4 = nn.Conv2d(channel_in + 3 * gc, gc, 3, 1, 1, bias=bias)
        self.conv5 = nn.Conv2d(channel_in + 4 * gc, channel_out, 3, 1, 1, bias=bias)
        self.lrelu = nn.LeakyReLU(negative_slope=0.2, inplace=True)

        # if init == 'xavier':
        #     mutil.initialize_weights_xavier([self.conv1, self.conv2, self.conv3, self.conv4], 0.1)
        # else:
        #     mutil.initialize_weights([self.conv1, self.conv2, self.conv3, self.conv4], 0.1)
        # mutil.initialize_weights(self.conv5, 0)

    def forward(self, x):
        x1 = self.lrelu(self.conv1(x))
        x2 = self.lrelu(self.conv2(torch.cat((x, x1), 1)))
        x3 = self.lrelu(self.conv3(torch.cat((x, x1, x2), 1)))
        x4 = self.lrelu(self.conv4(torch.cat((x, x1, x2, x3), 1)))
        x5 = self.conv5(torch.cat((x, x1, x2, x3, x4), 1))
        return x5

class Upsampler(nn.Sequential):
    def __init__(self, conv, scale, n_feat, bn=False, act=False, bias=True):

        m = []
        if (scale & (scale - 1)) == 0:    # Is scale = 2^n?
            for _ in range(int(math.log(scale, 2))):
                m.append(conv(n_feat, 4 * n_feat, 3, bias))
                m.append(nn.PixelShuffle(2))
                if bn: m.append(nn.BatchNorm2d(n_feat))
                if act: m.append(act())
        elif scale == 3:
            m.append(conv(n_feat, 9 * n_feat, 3, bias))
            m.append(nn.PixelShuffle(3))
            if bn: m.append(nn.BatchNorm2d(n_feat))
            if act: m.append(act())
        else:
            raise NotImplementedError

        super(Upsampler, self).__init__(*m)


class Conv_up(nn.Module):
    def __init__(self, c_in, mid_c, up_factor):
        super(Conv_up, self).__init__()

        body = [nn.Conv2d(in_channels=c_in, out_channels=mid_c, kernel_size=3, padding=3 // 2), nn.ReLU(),
                # nn.Conv2d(in_channels=64, out_channels=64, kernel_size=3, padding=3 // 2), nn.ReLU(),
                # nn.Conv2d(in_channels=64, out_channels=64, kernel_size=3, padding=3 // 2), nn.ReLU(),
                ]
        self.body = nn.Sequential(*body)
        conv = default_conv
        ## x3 00
        ## x2 11
        if up_factor == 2:
            modules_tail = [
                # nn.ConvTranspose2d(64, 64, kernel_size=3, stride=up_factor, padding=1, output_padding=1),
                nn.Upsample(scale_factor=2),
                conv(mid_c, c_in,3),
                conv(c_in, c_in, 3)]

        elif up_factor == 3:
            modules_tail = [
                # nn.ConvTranspose2d(64, 64, kernel_size=3, stride=up_factor, padding=0, output_padding=0),
                nn.Upsample(scale_factor=3),
                conv(mid_c, c_in,3),
                conv(c_in, c_in, 3)]

        elif up_factor == 4:
            modules_tail = [
                # nn.ConvTranspose2d(64, 64, kernel_size=3, stride=2, padding=1, output_padding=1),
                # nn.ConvTranspose2d(64, 64, kernel_size=3, stride=2, padding=1, output_padding=1),
                nn.Upsample(scale_factor=4),
                conv(mid_c, c_in,3),
                conv(c_in, c_in, 3)]
        self.tail = nn.Sequential(*modules_tail)

    def forward(self, input):

        out = self.body(input)
        out = self.tail(out)
        return out


class Conv_down(nn.Module):
    def __init__(self, c_in,mid_c, up_factor):
        super(Conv_down, self).__init__()

        body = [nn.Conv2d(in_channels=c_in, out_channels=mid_c, kernel_size=3, padding=3 // 2), nn.ReLU(),
                # nn.Conv2d(in_channels=64, out_channels=64, kernel_size=3, padding=3 // 2), nn.ReLU(),
                # nn.Conv2d(in_channels=64, out_channels=64, kernel_size=3, padding=3 // 2), nn.ReLU(),
                ]
        self.body = nn.Sequential(*body)
        conv = default_conv
        if up_factor == 4:
            modules_tail = [
                # nn.Conv2d(in_channels=64, out_channels=64, kernel_size=3, padding=1, stride=2),
                # nn.Conv2d(in_channels=64, out_channels=64, kernel_size=3, padding=1, stride=2),
                nn.MaxPool2d(4),
                conv(mid_c, c_in,3),
                conv(c_in, c_in, 3)]

        elif up_factor == 3:
            modules_tail = [
                # nn.Conv2d(in_channels=64, out_channels=64, kernel_size=3, padding=1, stride=up_factor),
                nn.MaxPool2d(3),
                conv(mid_c, c_in,3),
                conv(c_in, c_in, 3)]

        elif up_factor == 2:
            modules_tail = [
                # nn.Conv2d(in_channels=64, out_channels=64, kernel_size=3, padding=1, stride=up_factor),
                nn.MaxPool2d(2),
                conv(mid_c, c_in,3),
                conv(c_in, c_in, 3)]
                
        self.tail = nn.Sequential(*modules_tail)

    def forward(self, input):

        out = self.body(input)
        out = self.tail(out)
        return out

class HinResBlock(nn.Module):
    def __init__(self, channel_in, channel_out):
        super(HinResBlock, self).__init__()
        feature = 64
        self.conv1 = nn.Conv2d(channel_in, feature, kernel_size=3, padding=1)
        self.relu1 = nn.LeakyReLU(negative_slope=0.2, inplace=True)
        self.conv2 = nn.Conv2d(feature, feature, kernel_size=3, padding=1)
        self.conv3 = nn.Conv2d((feature+channel_in), channel_out, kernel_size=3, padding=1)
        self.norm = nn.InstanceNorm2d(feature // 2, affine=True)

    def forward(self, x):
        residual = self.relu1(self.conv1(x))

        out_1, out_2 = torch.chunk(residual, 2, dim=1)
        residual = torch.cat([self.norm(out_1), out_2], dim=1)

        residual = self.relu1(self.conv2(residual))
        input = torch.cat((x, residual), dim=1)
        out = self.conv3(input)
        return out


def subnet(net_structure, init='xavier'):
    def constructor(channel_in, channel_out):
        if net_structure == 'DBNet':
            if init == 'xavier':
                return DenseBlock(channel_in, channel_out, init)
            else:
                return DenseBlock(channel_in, channel_out)
        elif net_structure == 'Resnet':
            return ResBlock(channel_in, channel_out)
        elif net_structure == 'HinResnet':
            return HinResBlock(channel_in, channel_out)
        else:
            return None
    return constructor

class InvBlock(nn.Module):
    def __init__(self, subnet_constructor, channel_num1, channel_num2, clamp=0.8):
        super(InvBlock, self).__init__()
        # channel_num: 3
        # channel_split_num: 1

        self.split_len1 = channel_num1  # 1
        self.split_len2 = channel_num2  # 2

        self.clamp = clamp

        self.F = subnet_constructor(self.split_len2, self.split_len1)
        self.G = subnet_constructor(self.split_len1, self.split_len2)
        self.H = subnet_constructor(self.split_len1, self.split_len2)

        #in_channels = 3
        # self.invconv = InvertibleConv1x1(channel_num, LU_decomposed=True)
        # self.flow_permutation = lambda z, logdet, rev: self.invconv(z, logdet, rev)

    # def forward(self, x, rev=False):
    def forward(self, x, rev=False):
        x1, x2 = (x.narrow(1, 0, self.split_len1), x.narrow(1, self.split_len1, self.split_len2))
        if not rev:
            # invert1x1conv
            # x, logdet = self.flow_permutation(x, logdet=0, rev=False)

            # split to 1 channel and 2 channel.
            # x1, x2 = (x.narrow(1, 0, self.split_len1), x.narrow(1, self.split_len1, self.split_len2))

            y1 = x1 + self.F(x2)  # 1 channel
            self.s = self.clamp * (torch.sigmoid(self.H(y1)) * 2 - 1)
            y2 = x2.mul(torch.exp(self.s)) + self.G(y1)  # 2 channel
            
        else:
            # split.
            # x1, x2 = (x.narrow(1, 0, self.split_len1), x.narrow(1, self.split_len1, self.split_len2))
            self.s = self.clamp * (torch.sigmoid(self.H(x1)) * 2 - 1)
            y2 = (x2 - self.G(x1)).div(torch.exp(self.s))
            y1 = x1 - self.F(y2)

            # x = torch.cat((y1, y2), 1)
            # print("rev_inn")
            # inv permutation
            # out = x
        out = torch.cat((y1, y2), 1)
        return out
models/modules/common.py
lpcccc-cv/MC-VarNet
from turtle import forward
import torch
import torch.nn as nn
import numpy as np
import torch.nn.functional as fnn
from torch.autograd import Variable
import math

class CharbonnierLoss(nn.Module):
    """Charbonnier Loss (L1)"""

    def __init__(self, eps=1e-6):
        super(CharbonnierLoss, self).__init__()
        self.eps = eps

    def forward(self, x, y):
        diff = x - y
        loss = torch.sum(torch.sqrt(diff * diff + self.eps))
        return loss

class fftLoss(nn.Module):
    def __init__(self):
        super(fftLoss, self).__init__()

    def forward(self, x, y, mask):
        diff = torch.fft.fft2(x) - torch.fft.fft2(y)
        loss = torch.mean(abs(diff*mask))
        return loss

class ncc_loss(nn.Module):
    def __init__(self):
        super(ncc_loss, self).__init__()

    def compute_local_sums(self, I, J, filt, stride, padding, win):
        I2 = I * I
        J2 = J * J
        IJ = I * J

        I_sum = fnn.conv2d(I, filt, stride=stride, padding=padding)
        J_sum = fnn.conv2d(J, filt, stride=stride, padding=padding)
        I2_sum = fnn.conv2d(I2, filt, stride=stride, padding=padding)
        J2_sum = fnn.conv2d(J2, filt, stride=stride, padding=padding)
        IJ_sum = fnn.conv2d(IJ, filt, stride=stride, padding=padding)

        win_size = np.prod(win)
        u_I = I_sum / win_size
        u_J = J_sum / win_size

        cross = IJ_sum - u_J * I_sum - u_I * J_sum + u_I * u_J * win_size
        I_var = I2_sum - 2 * u_I * I_sum + u_I * u_I * win_size
        J_var = J2_sum - 2 * u_J * J_sum + u_J * u_J * win_size

        return I_var, J_var, cross

    def forward(self, I, J, win=None):
        ndims = len(list(I.size())) - 2
        assert ndims in [1, 2, 3], "volumes should be 1 to 3 dimensions. found: %d" % ndims

        if win is None:
            win = [9] * ndims
        else:
            win = win * ndims
        conv_fn = getattr(fnn, 'conv%dd' % ndims)
        I2 = I * I
        J2 = J * J
        IJ = I * J
        sum_filt = torch.ones([1, 1, *win]).to("cuda")
        pad_no = math.floor(win[0] / 2)
        if ndims == 1:
            stride = (1)
            padding = (pad_no)
        elif ndims == 2:
            stride = (1, 1)
            padding = (pad_no, pad_no)
        else:
            stride = (1, 1, 1)
            padding = (pad_no, pad_no, pad_no)
        I_var, J_var, cross = self.compute_local_sums(I, J, sum_filt, stride, padding, win)
        cc = cross * cross / (I_var * J_var + 1e-5)

        return torch.mean(cc)


def build_gauss_kernel(size=5, sigma=1.0, n_channels=1, cuda=False):
    if size % 2 != 1:
        raise ValueError("kernel size must be uneven")
    grid = np.float32(np.mgrid[0:size,0:size].T)
    gaussian = lambda x: np.exp((x - size//2)**2/(-2*sigma**2))**2
    kernel = np.sum(gaussian(grid), axis=2)
    kernel /= np.sum(kernel)
    # repeat same kernel across depth dimension
    kernel = np.tile(kernel, (n_channels, 1, 1))
    # conv weight should be (out_channels, groups/in_channels, h, w), 
    # and since we have depth-separable convolution we want the groups dimension to be 1
    kernel = torch.FloatTensor(kernel[:, None, :, :])
    if cuda:
        kernel = kernel.cuda()
    return Variable(kernel, requires_grad=False)


def conv_gauss(img, kernel):
    """ convolve img with a gaussian kernel that has been built with build_gauss_kernel """
    n_channels, _, kw, kh = kernel.shape
    img = fnn.pad(img, (kw//2, kh//2, kw//2, kh//2), mode='replicate')
    return fnn.conv2d(img, kernel, groups=n_channels)


def laplacian_pyramid(img, kernel, max_levels=5):
    current = img
    pyr = []

    for level in range(max_levels):
        filtered = conv_gauss(current, kernel)
        diff = current - filtered
        pyr.append(diff)
        current = fnn.avg_pool2d(filtered, 2)

    pyr.append(current)
    return pyr

class LapLoss(nn.Module):
    def __init__(self, max_levels=5, k_size=5, sigma=2.0):
        super(LapLoss, self).__init__()
        self.max_levels = max_levels
        self.k_size = k_size
        self.sigma = sigma
        self._gauss_kernel = None
        
    def forward(self, input, target):
        # input shape :[B, N, C, H, W]
        if len(input.shape) == 5:
            B,N,C,H,W = input.size()
            input = input.view(-1, C, H , W)
            target = target.view(-1, C, H, W)
        if self._gauss_kernel is None or self._gauss_kernel.shape[1] != input.shape[1]:
            self._gauss_kernel = build_gauss_kernel(
                size=self.k_size, sigma=self.sigma, 
                n_channels=input.shape[1], cuda=input.is_cuda
            )
        pyr_input  = laplacian_pyramid(input, self._gauss_kernel, self.max_levels)
        pyr_target = laplacian_pyramid(target, self._gauss_kernel, self.max_levels)
        return sum(fnn.l1_loss(a, b) for a, b in zip(pyr_input, pyr_target))



# if __name__ == "__main__":  

models/modules/loss.py
lpcccc-cv/MC-VarNet
import torch
import torch.nn as nn
import torch.nn.init as init
import torch.nn.functional as F


def initialize_weights(net_l, scale=1):
    if not isinstance(net_l, list):
        net_l = [net_l]
    for net in net_l:
        for m in net.modules():
            if isinstance(m, nn.Conv2d):
                init.kaiming_normal_(m.weight, a=0, mode='fan_in')
                m.weight.data *= scale  # for residual block
                if m.bias is not None:
                    m.bias.data.zero_()
            elif isinstance(m, nn.Linear):
                init.kaiming_normal_(m.weight, a=0, mode='fan_in')
                m.weight.data *= scale
                if m.bias is not None:
                    m.bias.data.zero_()
            elif isinstance(m, nn.BatchNorm2d):
                init.constant_(m.weight, 1)
                init.constant_(m.bias.data, 0.0)


def make_layer(block, n_layers):
    layers = []
    for _ in range(n_layers):
        layers.append(block())
    return nn.Sequential(*layers)


class ResidualBlock_noBN(nn.Module):
    '''Residual block w/o BN
    ---Conv-ReLU-Conv-+-
     |________________|
    '''

    def __init__(self, nf=64):
        super(ResidualBlock_noBN, self).__init__()
        self.conv1 = nn.Conv2d(nf, nf, 3, 1, 1, bias=True)
        self.conv2 = nn.Conv2d(nf, nf, 3, 1, 1, bias=True)

        # initialization
        initialize_weights([self.conv1, self.conv2], 0.1)

    def forward(self, x):
        identity = x
        out = F.relu(self.conv1(x), inplace=True)
        out = self.conv2(out)
        return identity + out

class ResidualBlock_noBN_noAct(nn.Module):
    '''Residual block w/o BN
    ---Conv-ReLU-Conv-+-
     |________________|
    '''

    def __init__(self, nf=64):
        super(ResidualBlock_noBN_noAct, self).__init__()
        self.conv1 = nn.Conv2d(nf, nf, 3, 1, 1, bias=True)
        self.conv2 = nn.Conv2d(nf, nf, 3, 1, 1, bias=True)

        # initialization
        initialize_weights([self.conv1, self.conv2], 0.1)

    def forward(self, x):
        identity = x
        out = F.relu(self.conv1(x), inplace=True)
        out = self.conv2(out)
        return identity + out

class BasicConv(nn.Module):
    def __init__(self, in_channel, out_channel, kernel_size, stride, bias=False, norm=False, relu=True, transpose=False,
                 channel_shuffle_g=0, norm_method=nn.BatchNorm2d, groups=1):
        super(BasicConv, self).__init__()
        self.channel_shuffle_g = channel_shuffle_g
        self.norm = norm
        if bias and norm:
            bias = False

        padding = kernel_size // 2
        layers = list()
        if transpose:
            padding = kernel_size // 2 - 1
            layers.append(
                nn.ConvTranspose2d(in_channel, out_channel, kernel_size, padding=padding, stride=stride, bias=bias, groups=groups))
        else:
            layers.append(
                nn.Conv2d(in_channel, out_channel, kernel_size, padding=padding, stride=stride, bias=bias, groups=groups))
        if norm:
            layers.append(norm_method(out_channel))
        elif relu:
            layers.append(nn.ReLU(inplace=True))

        self.main = nn.Sequential(*layers)

    def forward(self, x):
        return self.main(x)

class ResBlock_fft_bench(nn.Module):
    def __init__(self, n_feat=64, norm='backward'): # 'ortho'
        super(ResBlock_fft_bench, self).__init__()
        self.main = nn.Sequential(
            BasicConv(n_feat, n_feat, kernel_size=3, stride=1, relu=True),
            BasicConv(n_feat, n_feat, kernel_size=3, stride=1, relu=False)
        )
        self.main_fft = nn.Sequential(
            BasicConv(n_feat*2, n_feat*2, kernel_size=1, stride=1, relu=True),
            BasicConv(n_feat*2, n_feat*2, kernel_size=1, stride=1, relu=False)
        )
        self.dim = n_feat
        self.norm = norm
    def forward(self, x):
        _, _, H, W = x.shape
        dim = 1
        y = torch.fft.rfft2(x, norm=self.norm)
        y_imag = y.imag
        y_real = y.real
        y_f = torch.cat([y_real, y_imag], dim=dim)
        y = self.main_fft(y_f)
        y_real, y_imag = torch.chunk(y, 2, dim=dim)
        y = torch.complex(y_real, y_imag)
        y = torch.fft.irfft2(y, s=(H, W), norm=self.norm)
        return self.main(x) + x + y

def flow_warp(x, flow, interp_mode='bilinear', padding_mode='zeros'):
    """Warp an image or feature map with optical flow
    Args:
        x (Tensor): size (N, C, H, W)
        flow (Tensor): size (N, H, W, 2), normal value
        interp_mode (str): 'nearest' or 'bilinear'
        padding_mode (str): 'zeros' or 'border' or 'reflection'

    Returns:
        Tensor: warped image or feature map
    """
    print(x.size()[-2:])
    print(flow.size()[1:3])
    assert x.size()[-2:] == flow.size()[1:3]
    B, C, H, W = x.size()
    # mesh grid
    grid_y, grid_x = torch.meshgrid(torch.arange(0, H), torch.arange(0, W))
    grid = torch.stack((grid_x, grid_y), 2).float()  # W(x), H(y), 2
    grid.requires_grad = False
    grid = grid.type_as(x)
    vgrid = grid + flow
    # scale grid to [-1,1]
    vgrid_x = 2.0 * vgrid[:, :, :, 0] / max(W - 1, 1) - 1.0
    vgrid_y = 2.0 * vgrid[:, :, :, 1] / max(H - 1, 1) - 1.0
    vgrid_scaled = torch.stack((vgrid_x, vgrid_y), dim=3)
    output = F.grid_sample(x, vgrid_scaled, mode=interp_mode, padding_mode=padding_mode)
    return output


class GridNet(nn.Module):
    def __init__(self):
        super(GridNet, self).__init__()
        self.lateral = nn.ModuleList([
            nn.ModuleList([
                LateralBlock(24),
                LateralBlock(24),
                LateralBlock(24),
                LateralBlock(24),
                LateralBlock(24),
                LateralBlock(24),
            ]),
            nn.ModuleList([
                LateralBlock(48),
                LateralBlock(48),
                LateralBlock(48),
                LateralBlock(48),
                LateralBlock(48),
            ]),
            nn.ModuleList([
                LateralBlock(96),
                LateralBlock(96),
                LateralBlock(96),
                LateralBlock(96),
                LateralBlock(96),
            ])
        ])
        self.down = nn.ModuleList([
            nn.ModuleList([
                DownBlock(24, 48),
                DownBlock(24, 48),
                DownBlock(24, 48),
            ]),
            nn.ModuleList([
                DownBlock(48, 96),
                DownBlock(48, 96),
                DownBlock(48, 96),
            ])
        ])
        self.up = nn.ModuleList([
            nn.ModuleList([
                UpBlock(48, 24),
                UpBlock(48, 24),
                UpBlock(48, 24),
            ]),
            nn.ModuleList([
                UpBlock(96, 48),
                UpBlock(96, 48),
                UpBlock(96, 48),
            ])
        ])
        self.compress = nn.ModuleList([
            nn.Conv2d(24*3, 24, 3, 1, 1),
            nn.Conv2d(48*3, 48, 3, 1, 1),
            nn.Conv2d(96*3, 96, 3, 1, 1)
        ])
        self.to_RGB = nn.Sequential(
            nn.PReLU(),
            nn.Conv2d(24, 1, 3, 1, 1)
        )

    def forward(self, pyramid):
        x_0_0, x_1_0, x_2_0 = pyramid

        # compress dim 32 x 2 -> 32
        x_0_0 = self.compress[0](x_0_0)
        x_1_0 = self.compress[1](x_1_0)
        x_2_0 = self.compress[2](x_2_0)

        # first half: down & lateral
        x_0_1 = self.lateral[0][0](x_0_0)
        x_0_2 = self.lateral[0][1](x_0_1)
        x_0_3 = self.lateral[0][2](x_0_2)

        x_1_0 = x_1_0 + self.down[0][0](x_0_0)
        x_2_0 = x_2_0 + self.down[1][0](x_1_0)

        x_1_1 = self.lateral[1][0](x_1_0)
        x_2_1 = self.lateral[2][0](x_2_0)

        x_1_1 = x_1_1 + self.down[0][1](x_0_1)
        x_2_1 = x_2_1 + self.down[1][1](x_1_1)

        x_1_2 = self.lateral[1][1](x_1_1)
        x_2_2 = self.lateral[2][1](x_2_1)

        x_1_2 = x_1_2 + self.down[0][2](x_0_2)
        x_2_2 = x_2_2 + self.down[1][2](x_1_2)

        x_1_3 = self.lateral[1][2](x_1_2)
        x_2_3 = self.lateral[2][2](x_2_2)

        # second half: up & lateral
        x_2_4 = self.lateral[2][3](x_2_3)
        x_2_5 = self.lateral[2][4](x_2_4)

        x_1_3 = x_1_3 + self.up[1][0](x_2_3)
        x_0_3 = x_0_3 + self.up[0][0](x_1_3)

        x_1_4 = self.lateral[1][3](x_1_3)
        x_0_4 = self.lateral[0][3](x_0_3)

        x_1_4 = x_1_4 + self.up[1][1](x_2_4)
        x_0_4 = x_0_4 + self.up[0][1](x_1_4)

        x_1_5 = self.lateral[1][4](x_1_4)
        x_0_5 = self.lateral[0][4](x_0_4)

        x_1_5 = x_1_5 + self.up[1][2](x_2_5)
        x_0_5 = x_0_5 + self.up[0][2](x_1_5)

        # final synthesis
        output = self.lateral[0][5](x_0_5)
        output = self.to_RGB(output)

        return output



class LateralBlock(nn.Module):
    def __init__(self, dim):
        super(LateralBlock, self).__init__()
        self.layers = nn.Sequential(
            nn.PReLU(),
            nn.Conv2d(dim, dim, 3, 1, 1),
            nn.PReLU(),
            nn.Conv2d(dim, dim, 3, 1, 1)
        )

    def forward(self, x):
        res = x
        x = self.layers(x)
        return x + res


class DownBlock(nn.Module):
    def __init__(self, in_dim, out_dim):
        super(DownBlock, self).__init__()
        self.layers = nn.Sequential(
            nn.PReLU(),
            nn.Conv2d(in_dim, out_dim, 3, 2, 1),
            nn.PReLU(),
            nn.Conv2d(out_dim, out_dim, 3, 1, 1)
        )

    def forward(self, x):
        return self.layers(x)


class UpBlock(nn.Module):
    def __init__(self, in_dim, out_dim):
        super(UpBlock, self).__init__()
        self.layers = nn.Sequential(
            nn.PReLU(),
            nn.Conv2d(in_dim, out_dim, 3, 1, 1),
            nn.PReLU(),
            nn.Conv2d(out_dim, out_dim, 3, 1, 1)
        )

    def forward(self, x):
        x = F.interpolate(x, scale_factor=2.0, mode='bilinear', align_corners=False)
        return self.layers(x)


models/modules/module_util.py
lpcccc-cv/MC-VarNet
import models.modules.VANet as VANet
####################
# define network
####################
# Generator
def define_G(opt):
    opt_net = opt['network_G']
    which_model = opt_net['which_model_G']

    if which_model == 'VANet':
        netG = VANet.VANet()
        
    else:
        raise NotImplementedError('Generator model [{:s}] not recognized'.format(which_model))

    return netG

models/networks.py
lpcccc-cv/MC-VarNet
import logging
from collections import OrderedDict
import torch.nn.functional as F
import torch
import torch.nn as nn
from torch.nn.parallel import DataParallel, DistributedDataParallel
import models.networks as networks
import models.lr_scheduler as lr_scheduler
from .base_model import BaseModel
from models.modules.loss import CharbonnierLoss, LapLoss, ncc_loss, fftLoss
from torchvision import models
import cv2
from matplotlib import pyplot as plt
import os
import numpy as np

logger = logging.getLogger('base')


def vis_img(img, fname, ftype ,output_dir):
    os.makedirs(output_dir, exist_ok=True)
    plt.figure()
    plt.imshow(img, cmap='gray')
    figname = fname + '_' + ftype + '.png'
    figpath = os.path.join(output_dir, figname)
    plt.savefig(figpath)

class BaseModel(BaseModel):
    def __init__(self, opt):
        super(BaseModel, self).__init__(opt)

        if opt['dist']:
            self.rank = torch.distributed.get_rank()
        else:
            self.rank = -1  # non dist training
        train_opt = opt['train']
        self.dwt = opt["dwt"]
        self.which_dataset = opt["mode"]
        self.modal_1 = opt['data1']
        self.kspace_loss = opt['k_loss']

        # define network and load pretrained models
        self.netG = networks.define_G(opt).to(self.device)
        
        if opt['dist']:
            self.netG = DistributedDataParallel(self.netG, device_ids=[torch.cuda.current_device()], find_unused_parameters=True)
        else:
            self.netG = DataParallel(self.netG)
        # print network
        self.print_network()
        self.load()

        if self.is_train:
            self.netG.train()

            #### loss
            loss_type = train_opt['pixel_criterion']
            if loss_type == 'l1':
                self.cri_pix = nn.L1Loss(reduction='mean').to(self.device)
            elif loss_type == 'l2':
                self.cri_pix = nn.MSELoss(reduction='sum').to(self.device)
            elif loss_type == 'cb':
                self.cri_pix = CharbonnierLoss().to(self.device)
            elif loss_type == 'lp':
                self.cri_pix = LapLoss(max_levels=5).to(self.device)

                raise NotImplementedError('Loss type [{:s}] is not recognized.'.format(loss_type))
            self.l_pix_w = train_opt['pixel_weight']
            self.l1_loss = nn.L1Loss(reduction='mean').to(self.device)
            self.fft_loss = fftLoss().to(self.device)


            #### optimizers
            wd_G = train_opt['weight_decay_G'] if train_opt['weight_decay_G'] else 0
            optim_params = []
            for k, v in self.netG.named_parameters():
                if opt['datasets']['train']['use_time'] == True and 'TMB' not in k:
                    v.requires_grad = False
                if v.requires_grad:
                    optim_params.append(v)
                else:
                    if self.rank <= 0:
                        logger.warning('Params [{:s}] will not optimize.'.format(k))

            self.optimizer_G = torch.optim.Adam(optim_params, lr=train_opt['lr_G'],
                                                weight_decay=wd_G,
                                                betas=(train_opt['beta1'], train_opt['beta2']))
            # self.optimizer_G = torch.optim.SGD(optim_params, lr=train_opt['lr_G'], momentum=0.9, weight_decay=wd_G)
            self.optimizers.append(self.optimizer_G)
            #### schedulers
            if train_opt['lr_scheme'] == 'MultiStepLR':
                for optimizer in self.optimizers:
                    self.schedulers.append(
                        lr_scheduler.MultiStepLR_Restart(optimizer, train_opt['lr_steps'],
                                                         restarts=train_opt['restarts'],
                                                         weights=train_opt['restart_weights'],
                                                         gamma=train_opt['lr_gamma'],
                                                         clear_state=train_opt['clear_state']))
            elif train_opt['lr_scheme'] == 'CosineAnnealingLR_Restart':
                for optimizer in self.optimizers:
                    self.schedulers.append(
                        lr_scheduler.CosineAnnealingLR_Restart(
                            optimizer, train_opt['T_period'], eta_min=train_opt['eta_min'],
                            restarts=train_opt['restarts'], weights=train_opt['restart_weights']))
            else:
                raise NotImplementedError()

            self.log_dict = OrderedDict()

    
    def feed_data(self, data):   

        self.im1_lr = data['im1_LQ'].to(self.device)
        self.im1_gt = data['im1_GT'].to(self.device)
        self.im2_lr = data['im2_LQ'].to(self.device)
        self.im2_gt = data['im2_GT'].to(self.device)
        self.mask = data['mask'].to(self.device)


    def set_params_lr_zero(self):
        # fix normal module
        self.optimizers[0].param_groups[0]['lr'] = 0

    def optimize_parameters(self):
        self.optimizer_G.zero_grad()
        if self.modal_1:
            self.fake_1 = self.netG(self.im1_lr, self.im2_gt, self.mask)
            l_pix = self.cri_pix(self.fake_1, self.im1_gt)
            # self.fake_1, self.fake_2 = self.netG(self.im1_lr, self.im2_gt, self.mask)
            # l_pix = self.cri_pix(self.fake_1, self.im1_gt) + 0.1 * self.cri_pix(self.fake_2, self.im2_gt)
        else:
            self.fake_2 = self.netG(self.im2_lr, self.im1_gt, self.mask)
            l_pix = self.cri_pix(self.fake_2, self.im2_gt)
            # self.fake_2, self.fake_1 = self.netG(self.im2_lr, self.im1_gt, self.mask)
            # l_pix = self.cri_pix(self.fake_2, self.im2_gt) + 0.1 * self.cri_pix(self.fake_1, self.im1_gt)

        total_loss = l_pix
        total_loss.backward()
        self.optimizer_G.step()
        self.log_dict['l_pix'] = l_pix.item()


    def test(self):
        self.netG.eval()
        with torch.no_grad():
            if self.modal_1:
                self.fake_H = self.netG(self.im1_lr, self.im2_gt, self.mask)
                # self.fake_H,_ = self.netG(self.im1_lr, self.im2_gt, self.mask)
            else:
                self.fake_H = self.netG(self.im2_lr, self.im1_gt, self.mask)
                # self.fake_H ,_= self.netG(self.im2_lr, self.im1_gt, self.mask)
        self.netG.train()

    def get_current_log(self):
        return self.log_dict

    
    def get_current_visuals(self, need_GT=True):

        out_dict = OrderedDict()
        out_dict['im1_restore'] = self.fake_H.detach().float().cpu()
        if self.modal_1:
            out_dict['im1_GT'] = self.im1_gt.detach().float().cpu() 
            if self.which_dataset == 'fastmri':
                out_dict['im1_restore'] = out_dict['im1_restore']*self.im1_std.cpu() + self.im1_mean.cpu()
                out_dict['im1_GT'] = out_dict['im1_GT']*self.im1_std.cpu() + self.im1_mean.cpu()
        else:
            out_dict['im1_GT'] = self.im2_gt.detach().float().cpu() 
            if self.which_dataset == 'fastmri':
                # print(out_dict['im1_restore'].shape, self.im2_std.cpu().shape)
                out_dict['im1_restore'] = out_dict['im1_restore']*self.im2_std.cpu() + self.im2_mean.cpu()
                out_dict['im1_GT'] = out_dict['im1_GT']*self.im2_std.cpu() + self.im2_mean.cpu()
  
        return out_dict

    def print_network(self):
        s, n = self.get_network_description(self.netG)
        if isinstance(self.netG, nn.DataParallel):
            net_struc_str = '{} - {}'.format(self.netG.__class__.__name__,
                                             self.netG.module.__class__.__name__)
        else:
            net_struc_str = '{}'.format(self.netG.__class__.__name__)
        if self.rank <= 0:
            logger.info('Network G structure: {}, with parameters: {:,d}'.format(net_struc_str, n))
            logger.info(s)

    def load(self):
        load_path_G = self.opt['path']['pretrain_model_G']
        if load_path_G is not None:
            logger.info('Loading model for G [{:s}] ...'.format(load_path_G))
            self.load_network(load_path_G, self.netG, self.opt['path']['strict_load'])

    def save(self, iter_label, epoch=None):
        if epoch != None:
            self.save_network(self.netG, 'G', iter_label, str(epoch))
        else:
            self.save_network(self.netG, 'G', iter_label)

models/ref_rec_base_model.py
lpcccc-cv/MC-VarNet
import os
import os.path as osp
import logging
import yaml
from utils.util import OrderedYaml
Loader, Dumper = OrderedYaml()


def parse(opt_path, is_train=True):
    with open(opt_path, mode='r') as f:
        opt = yaml.load(f, Loader=Loader)
    # export CUDA_VISIBLE_DEVICES
    gpu_list = ','.join(str(x) for x in opt['gpu_ids'])
    os.environ['CUDA_VISIBLE_DEVICES'] = gpu_list
    print('export CUDA_VISIBLE_DEVICES=' + gpu_list)

    opt['is_train'] = is_train
    if opt['distortion'] == 'sr' or opt['distortion'] == 'isr':
        scale = opt['scale']

    # datasets
    for phase, dataset in opt['datasets'].items():
        phase = phase.split('_')[0]
        dataset['phase'] = phase
        if opt['distortion'] == 'sr' or opt['distortion'] == 'isr':
            dataset['scale'] = scale
        is_lmdb = False
        if dataset.get('dataroot_GT', None) is not None:
            dataset['dataroot_GT'] = os.path.expanduser(dataset['dataroot_GT'])
            if dataset['dataroot_GT'].endswith('lmdb'):
                is_lmdb = True
        if dataset.get('dataroot_LQ', None) is not None:
            dataset['dataroot_LQ'] = os.path.expanduser(dataset['dataroot_LQ'])
            if dataset['dataroot_LQ'].endswith('lmdb'):
                is_lmdb = True
        dataset['data_type'] = 'lmdb' if is_lmdb else 'img'
        if dataset['mode'].endswith('mc'):  # for memcached
            dataset['data_type'] = 'mc'
            dataset['mode'] = dataset['mode'].replace('_mc', '')

    # path
    for key, path in opt['path'].items():
        if path and key in opt['path'] and key != 'strict_load':
            opt['path'][key] = osp.expanduser(path)
    opt['path']['root'] = osp.abspath(osp.join(__file__, osp.pardir, osp.pardir)) # opt['path']['root'] = osp.abspath(osp.join(__file__, osp.pardir, osp.pardir, osp.pardir))
    if is_train:
        experiments_root = os.path.join(opt['path']['root'], 'experiments', opt['name'])
        opt['path']['experiments_root'] = experiments_root
        opt['path']['models'] = os.path.join(experiments_root, 'models')
        opt['path']['training_state'] = os.path.join(experiments_root, 'training_state')
        opt['path']['log'] = experiments_root
        opt['path']['val_images'] = os.path.join(experiments_root, 'val_images')

        # # change some options for debug mode
        # if 'debug' in opt['name']:
        #     opt['train']['val_freq'] = 8
        #     opt['logger']['print_freq'] = 1
        #     opt['logger']['save_checkpoint_freq'] = 8
    else:  # test
        results_root = os.path.join(opt['path']['root'], 'results', opt['name'])
        opt['path']['results_root'] = results_root
        opt['path']['log'] = results_root

    # network
    if opt['distortion'] == 'sr' or opt['distortion'] == 'isr':
        opt['network_G']['scale'] = scale

    return opt


def dict2str(opt, indent_l=1):
    '''dict to string for logger'''
    msg = ''
    for k, v in opt.items():
        if isinstance(v, dict):
            msg += ' ' * (indent_l * 2) + k + ':[\n'
            msg += dict2str(v, indent_l + 1)
            msg += ' ' * (indent_l * 2) + ']\n'
        else:
            msg += ' ' * (indent_l * 2) + k + ': ' + str(v) + '\n'
    return msg


# convert to NoneDict, which return None for missing key.
class NoneDict(dict):
    def __missing__(self, key):
        return None


def dict_to_nonedict(opt):
    if isinstance(opt, dict):
        new_opt = dict()
        for key, sub_opt in opt.items():
            new_opt[key] = dict_to_nonedict(sub_opt)
        return NoneDict(**new_opt)
    elif isinstance(opt, list):
        return [dict_to_nonedict(sub_opt) for sub_opt in opt]
    else:
        return opt


def check_resume(opt, resume_iter):
    '''Check resume states and pretrain_model paths'''
    logger = logging.getLogger('base')
    if opt['path']['resume_state']:
        if opt['path'].get('pretrain_model_G', None) is not None or opt['path'].get(
                'pretrain_model_D', None) is not None:
            logger.warning('pretrain_model path will be ignored when resuming training.')

        if opt['path']['pretrain_model_G'] is None:
            opt['path']['pretrain_model_G'] = osp.join(opt['path']['models'], '{}_G.pth'.format(resume_iter))
        else:
            pass

        logger.info('Set [pretrain_model_G] to ' + opt['path']['pretrain_model_G'])
        if 'gan' in opt['model']:
            opt['path']['pretrain_model_D'] = osp.join(opt['path']['models'],
                                                       '{}_D.pth'.format(resume_iter))
            logger.info('Set [pretrain_model_D] to ' + opt['path']['pretrain_model_D'])

options/options.py
lpcccc-cv/MC-VarNet
import torch
from utils import util
from tqdm import tqdm
import cv2, os
import models.modules.VANet as VANet

import time

def main():
    mode = 'IXI'  #'IXI' or 'brain'
    dataset_opt = {}
    dataset_opt['task'] = 'rec'   # 'rec' or 'sr'
    dataset_opt['scale'] = 4
    dataset_opt['hr_in'] = True   # default=True
    dataset_opt['crop_size'] = 0
    dataset_opt['test_size'] = 240
    joint_rec = False
    save_result = True
    #### create train and val dataloader
    if mode == 'IXI':
        data1 = True
        from data.IXI_dataset import IXI_train as D
        dataset_opt['dataroot_GT'] = '/home/lpc/dataset/IXI/MC_MRI/test/T2'
    elif mode == 'brain':
        data1 = False
        from data.brain_dataset import brain_train as D
        dataset_opt['dataroot_GT'] = '/home/lpc/dataset/BrainTS/MCSR/T1_test'
    
    val_set = D(dataset_opt, train=False)
    val_loader = torch.utils.data.DataLoader(val_set, batch_size=1, shuffle=False, num_workers=1,pin_memory=True)
    print('Number of val images: {:d}'.format(len(val_set)))

    # creat model
    model_path = '~~~~~~~~~~'
    model = VANet.VANet().cuda()

    model_params = util.get_model_total_params(model)
    print('Model_params: ', model_params)
    model.load_state_dict(torch.load(model_path), strict=True)
    model.eval()

    with torch.no_grad():
        #### validation
        if joint_rec:
            avg_psnr_im1 = 0.0
            avg_psnr_im2 = 0.0
            avg_ssim_im1 = 0.0
            avg_ssim_im2 = 0.0
            avg_rmse_im1 = 0.0
            avg_rmse_im2 = 0.0
            idx = 0
            for i, val_data in enumerate(tqdm(val_loader)): 
                im1_lr = val_data['im1_LQ'].cuda()
                im1_gt = val_data['im1_GT'].cuda()
                im2_lr = val_data['im2_LQ'].cuda()
                im2_gt = val_data['im2_GT'].cuda()
                mask = val_data['mask'].cuda()

                sr_img_1, sr_img_2= model(im1_lr, im2_lr, mask)
                
                sr_img_1 = sr_img_1[0,0].cpu().detach().numpy()*255.
                im1_gt = im1_gt[0,0].cpu().detach().numpy()*255.
                sr_img_2 = sr_img_2[0,0].cpu().detach().numpy()*255.
                im2_gt = im2_gt[0,0].cpu().detach().numpy()*255.
                # calculate PSNR
                cur_psnr_im1 = util.calculate_psnr(sr_img_1, im1_gt)
                avg_psnr_im1 += cur_psnr_im1
                cur_psnr_im2 = util.calculate_psnr(sr_img_2, im2_gt)
                avg_psnr_im2 += cur_psnr_im2
                # ssim
                cur_ssim_im1 = util.calculate_ssim(sr_img_1, im1_gt)
                avg_ssim_im1 += cur_ssim_im1
                cur_ssim_im2 = util.calculate_ssim(sr_img_2, im2_gt)
                avg_ssim_im2 += cur_ssim_im2
                # rmse
                cur_rmse_im1 = util.calculate_rmse(sr_img_1, im1_gt)
                avg_rmse_im1 += cur_rmse_im1
                cur_rmse_im2 = util.calculate_rmse(sr_img_2, im2_gt)
                avg_rmse_im2 += cur_rmse_im2
                
                idx += 1

            avg_psnr_im1 = avg_psnr_im1 / idx
            avg_psnr_im2 = avg_psnr_im2 / idx
            avg_ssim_im1 = avg_ssim_im1 / idx
            avg_ssim_im2 = avg_ssim_im2 / idx
            avg_rmse_im1 = avg_rmse_im1 / idx
            avg_rmse_im2 = avg_rmse_im2 / idx
            # log
            print("# image1 Validation # PSNR: {:.6f}".format(avg_psnr_im1))
            print("# image2 Validation # PSNR: {:.6f}".format(avg_psnr_im2))
            print("# image1 Validation # SSIM: {:.6f}".format(avg_ssim_im1))
            print("# image2 Validation # SSIM: {:.6f}".format(avg_ssim_im2))
            print("# image1 Validation # RMSE: {:.6f}".format(avg_rmse_im1))
            print("# image2 Validation # RMSE: {:.6f}".format(avg_rmse_im2))
        else:
            avg_psnr_im1 = 0.0
            avg_ssim_im1 = 0.0
            avg_rmse_im1 = 0.0
            idx = 0
            time_begin = time.time()
            for i,val_data in enumerate(tqdm(val_loader)):
                im1_lr = val_data['im1_LQ'].cuda()
                im1_gt = val_data['im1_GT'].cuda()
                im2_lr = val_data['im2_LQ'].cuda()
                im2_gt = val_data['im2_GT'].cuda()
                mask = val_data['mask'].cuda()

                if data1: 
                    sr_img_1 = model(im1_lr, im2_gt, mask)
                    sr_img_1 = sr_img_1[0,0].cpu().detach().numpy()*255.
                    im1_gt = im1_gt[0,0].cpu().detach().numpy()*255.
               
                    # calculate PSNR
                    cur_psnr_im1 = util.calculate_psnr(sr_img_1, im1_gt)
                    avg_psnr_im1 += cur_psnr_im1
                    cur_ssim_im1 = util.calculate_ssim(sr_img_1, im1_gt)
                    avg_ssim_im1 += cur_ssim_im1
                    cur_rmse_im1 = util.calculate_rmse(sr_img_1, im1_gt)
                    avg_rmse_im1 += cur_rmse_im1
                else:
                    sr_img_1= model(im2_lr, im1_gt, mask)
                    sr_img_1 = sr_img_1[0,0].cpu().detach().numpy()*255.
                    im1_gt = im2_gt[0,0].cpu().detach().numpy()*255.
                    # calculate PSNR
                    cur_psnr_im1 = util.calculate_psnr(sr_img_1, im1_gt)
                    avg_psnr_im1 += cur_psnr_im1
                    cur_ssim_im1 = util.calculate_ssim(sr_img_1, im1_gt)
                    avg_ssim_im1 += cur_ssim_im1
                    cur_rmse_im1 = util.calculate_rmse(sr_img_1, im1_gt)
                    avg_rmse_im1 += cur_rmse_im1

                # save images
                # save_path_1 = '/home/lpc/program/MC_MRI/test_result/IXI/rec/'+'IXI_guided_T2_MDUNet_no_align'
                # if not os.path.exists(save_path_1):
                #     os.makedirs(save_path_1)
                # if save_result:
                #     cv2.imwrite(os.path.join(save_path_1, '{:08d}.png'.format(i+1)), sr_img_1)


                idx += 1 
            time_end = time.time()   
            avg_psnr_im1 = avg_psnr_im1 / idx
            avg_ssim_im1 = avg_ssim_im1 / idx
            avg_rmse_im1 = avg_rmse_im1 / idx
            # log
            print("# image1 Validation # PSNR: {:.6f}".format(avg_psnr_im1))
            print("# image1 Validation # SSIM: {:.6f}".format(avg_ssim_im1))
            print("# image1 Validation # RMSE: {:.6f}".format(avg_rmse_im1))
            print('Total time:', time_end-time_begin)

### CUDA_VISIBLE_DEVICES=2 python test_PSNR.py
if __name__ == '__main__':
    main()

test_PSNR.py
lpcccc-cv/MC-VarNet
import os
import math
import argparse
import random
import logging
import time
import torch
import torch.distributed as dist
import torch.multiprocessing as mp
from data.data_sampler import DistIterSampler

import options.options as option
from utils import util
from data import create_dataloader, create_dataset
from models import create_model

import data.util as data_util
import cv2
import numpy as np
import skimage.metrics as sm
import csv
from tqdm import tqdm

def init_dist(backend='nccl', **kwargs):
    ''' initialization for distributed training'''
    # if mp.get_start_method(allow_none=True) is None:
    if mp.get_start_method(allow_none=True) != 'spawn':
        mp.set_start_method('spawn')
    rank = int(os.environ['RANK'])
    num_gpus = torch.cuda.device_count()
    torch.cuda.set_device(rank % num_gpus)
    dist.init_process_group(backend=backend, **kwargs)

def main():
    #### options
    parser = argparse.ArgumentParser()
    parser.add_argument('-opt', type=str, help='Path to option YAML file.')
    parser.add_argument('--launcher', choices=['none', 'pytorch'], default='none',
                        help='job launcher')
    parser.add_argument('--local_rank', type=int, default=0)
    args = parser.parse_args()
    opt = option.parse(args.opt, is_train=True)
    
    #### distributed training settings
    if args.launcher == 'none':  # disabled distributed training
        opt['dist'] = False
        rank = -1
        print('Disabled distributed training.')
    else:
        opt['dist'] = True
        init_dist()
        world_size = torch.distributed.get_world_size()
        rank = torch.distributed.get_rank()

    #### loading resume state if exists
    if opt['path'].get('resume_state', None):
        # distributed resuming: all load into default GPU
        device_id = torch.cuda.current_device()
        resume_state = torch.load(opt['path']['resume_state'], map_location=lambda storage, loc: storage.cuda(device_id))
        option.check_resume(opt, resume_state['iter'])  # check resume options
    else:
        resume_state = None

    #### mkdir and loggers
    if rank <= 0:  # normal training (rank -1) OR distributed training (rank 0)
        if resume_state is None:
            util.mkdir_and_rename(
                opt['path']['experiments_root'])  # rename experiment folder if exists
            util.mkdirs((path for key, path in opt['path'].items() if not key == 'experiments_root'
                         and 'pretrain_model' not in key and 'resume' not in key))

        # config loggers. Before it, the log will not work
        util.setup_logger('base', opt['path']['log'], 'train_' + opt['name'], level=logging.INFO,
                          screen=True, tofile=True)
        logger = logging.getLogger('base')
        logger.info(option.dict2str(opt))
        # tensorboard logger
        if opt['use_tb_logger'] and 'debug' not in opt['name']:
            version = float(torch.__version__[0:3])
            if version >= 1.1:  # PyTorch 1.1
                from torch.utils.tensorboard import SummaryWriter
            else:
                logger.info(
                    'You are using PyTorch {}. Tensorboard will use [tensorboardX]'.format(version))
                from tensorboardX import SummaryWriter
            tb_logger = SummaryWriter(log_dir='../tb_logger/' + opt['name'])
    else:
        util.setup_logger('base', opt['path']['log'], 'train', level=logging.INFO, screen=True)
        logger = logging.getLogger('base')

    # convert to NoneDict, which returns None for missing keys
    opt = option.dict_to_nonedict(opt)
    which_model = opt['mode']
    #### random seed
    seed = opt['train']['manual_seed']
    if seed is None:
        seed = random.randint(1, 10000)
    if rank <= 0:
        logger.info('Random seed: {}'.format(seed))
    util.set_random_seed(seed)

    torch.backends.cudnn.benckmark = True
    # torch.backends.cudnn.deterministic = True

    #### create train and val dataloader
    dataset_ratio = 1.0 # dataset_ratio = 200  # enlarge the size of each epoch
    for phase, dataset_opt in opt['datasets'].items():
        if phase == 'train':
            train_set = create_dataset(dataset_opt, train=True)
            train_size = int(math.ceil(len(train_set) / dataset_opt['batch_size'])) # Iter number per epoch
            total_iters = int(opt['train']['niter'])
            total_epochs = int(math.ceil(total_iters / train_size))
            if opt['dist']:
                train_sampler = DistIterSampler(train_set, world_size, rank, dataset_ratio)
                total_epochs = int(math.ceil(total_iters / (train_size * dataset_ratio)))
            else:
                train_sampler = None
            train_loader = create_dataloader(train_set, dataset_opt, opt, train_sampler)
            if rank <= 0:
                logger.info('Number of train images: {:,d}, iters: {:,d}'.format(
                    len(train_set), train_size))
                logger.info('Total epochs needed: {:d} for iters {:,d}'.format(
                    total_epochs, total_iters))
        
        elif phase == 'val':
            val_set = create_dataset(dataset_opt, train=False)
            val_loader = create_dataloader(val_set, dataset_opt, opt, None)
            if rank <= 0:
                logger.info('Number of val images: {:d}'.format(len(val_set)))
        else:
            raise NotImplementedError('Phase [{:s}] is not recognized.'.format(phase))
    assert train_loader is not None
    assert val_loader is not None

    #### create model
    model = create_model(opt)

    #### resume training
    if resume_state:
        logger.info('Resuming training from epoch: {}, iter: {}.'.format(
            resume_state['epoch'], resume_state['iter']))

        start_epoch = resume_state['epoch'] + 1
        current_step = resume_state['iter']
        model.resume_training(resume_state)  # handle optimizers and schedulers
    else:
        current_step = 0
        start_epoch = 0

    #### training
    best_performance = 0
    logger.info('Start training from epoch: {:d}, iter: {:d}'.format(start_epoch, current_step))
    for epoch in range(start_epoch, total_epochs + 1):
        if opt['dist']:
            train_sampler.set_epoch(epoch)
        for _, train_data in enumerate(train_loader):
            current_step += 1
            if current_step > total_iters:
                break
                        #### training
            model.feed_data(train_data)
            model.optimize_parameters()

            #### log
            if current_step % opt['logger']['print_freq'] == 0:
                logs = model.get_current_log()
                message = '<epoch:{:3d}, iter:{:8,d}, lr:('.format(epoch, current_step)
                for v in model.get_current_learning_rate():
                    message += '{:.3e},'.format(v)
                message += ')>'
                for k, v in logs.items():
                    message += '{:s}: {:.4e} '.format(k, v)
                    # tensorboard logger
                    if opt['use_tb_logger'] and 'debug' not in opt['name']:
                        if rank <= 0:
                            tb_logger.add_scalar(k, v, current_step)
                if rank <= 0:
                    logger.info(message)
            #### save models and training states
            if rank <= 0:
                if current_step % opt['logger']['save_checkpoint_freq'] == 0:
                    logger.info('Saving models and training states at the end of step: ' + str(current_step))
                    model.save(current_step)
                    model.save_training_state(epoch, current_step)

            #### update learning rate
            model.update_learning_rate(current_step+1, warmup_iter=opt['train']['warmup_iter'])
            #### validation
            if opt['model'] == 'joint-rec':
                if current_step % opt["train"]["val_freq"] == 0 and rank <= 0:
                    avg_psnr_im1 = 0.0
                    avg_psnr_im2 = 0.0
                    idx = 0
                    for val_data in tqdm(val_loader): 
                        model.feed_data(val_data)
                        model.test()
                        
                        visuals = model.get_current_visuals()
                        img_num = visuals["im1_restore"].shape[0]
                        for i in range(img_num):
                            sr_img_1 = visuals["im1_restore"][i, 0, :, :]  # (1, w, h)
                            gt_img_1 = visuals["im1_GT"][i, 0, :, :]  # (1, w, h) 
                            sr_img_2 = visuals["im2_restore"][i, 0, :, :]  # (1, w, h)
                            gt_img_2 = visuals["im2_GT"][i, 0, :, :]  # (1, w, h)                  
                            # calculate PSNR
                            cur_psnr_im1 = util.calculate_psnr(sr_img_1.numpy()*255., gt_img_1.numpy()*255.)
                            avg_psnr_im1 += cur_psnr_im1
                            cur_psnr_im2 = util.calculate_psnr(sr_img_2.numpy()*255., gt_img_2.numpy()*255.)
                            avg_psnr_im2 += cur_psnr_im2

                        idx += img_num   
                    avg_psnr_im1 = avg_psnr_im1 / idx
                    avg_psnr_im2 = avg_psnr_im2 / idx
                    # log
                    logger.info("# image1 Validation # PSNR: {:.6f}".format(avg_psnr_im1))
                    logger.info("# image2 Validation # PSNR: {:.6f}".format(avg_psnr_im2))
            else:
                if current_step % opt["train"]["val_freq"] == 0 and rank <= 0:
                    avg_psnr_im1 = 0.0
                    idx = 0
                    for val_data in tqdm(val_loader): 
                        model.feed_data(val_data)
                        model.test()
            
                        visuals = model.get_current_visuals()
                        img_num = visuals["im1_GT"].shape[0]
                        # print(visuals["im1_restore"].shape, visuals["im1_GT"].shape)
                        for i in range(img_num):
                            sr_img_1 = visuals["im1_restore"][i, 0, :, :]  # (1, w, h)
                            gt_img_1 = visuals["im1_GT"][i, 0, :, :]  # (1, w, h) 
                                        
                            # calculate PSNR
                            if which_model == 'IXI' or which_model == 'brain' or which_model == 'knee':
                                cur_psnr_im1 = util.calculate_psnr(sr_img_1.numpy()*255., gt_img_1.numpy()*255.)
                            elif which_model == 'fastmri':
                                cur_psnr_im1 = util.calculate_psnr_fastmri(gt_img_1.numpy(), sr_img_1.numpy())
                        
                            avg_psnr_im1 += cur_psnr_im1

                        idx += img_num   
                    avg_psnr_im1 = avg_psnr_im1 / idx
                    # log
                    logger.info("# image1 Validation # PSNR: {:.6f}".format(avg_psnr_im1))


    if rank <= 0:
        logger.info('Saving the final model.')
        model.save('latest')
        logger.info('End of training.')


if __name__ == '__main__':
    main()

train.py
lpcccc-cv/MC-VarNet
# this code is modified from https://github.com/xinntao/EDVR/blob/master/codes/utils/util.py

import os
import sys
import time
import math
import torch.nn.functional as F
from datetime import datetime
import random
import logging
from collections import OrderedDict
import numpy as np
import cv2
import torch
from torchvision.utils import make_grid
from shutil import get_terminal_size
import glob
import re
from skimage.metrics import peak_signal_noise_ratio, structural_similarity

import yaml
try:
    from yaml import CLoader as Loader, CDumper as Dumper
except ImportError:
    from yaml import Loader, Dumper


def OrderedYaml():
    '''yaml orderedDict support'''
    _mapping_tag = yaml.resolver.BaseResolver.DEFAULT_MAPPING_TAG

    def dict_representer(dumper, data):
        return dumper.represent_dict(data.items())

    def dict_constructor(loader, node):
        return OrderedDict(loader.construct_pairs(node))

    Dumper.add_representer(OrderedDict, dict_representer)
    Loader.add_constructor(_mapping_tag, dict_constructor)
    return Loader, Dumper


####################
# miscellaneous
####################
def get_model_total_params(model):
    model_parameters = filter(lambda p: p.requires_grad, model.parameters())
    params = sum([np.prod(p.size()) for p in model_parameters])
    return (1.0*params/(1000*1000))

def get_timestamp():
    return datetime.now().strftime('%y%m%d-%H%M%S')


def mkdir(path):
    if not os.path.exists(path):
        os.makedirs(path)


def mkdirs(paths):
    if isinstance(paths, str):
        mkdir(paths)
    else:
        for path in paths:
            mkdir(path)


def mkdir_and_rename(path):
    if os.path.exists(path):
        new_name = path + '_archived_' + get_timestamp()
        print('Path already exists. Rename it to [{:s}]'.format(new_name))
        logger = logging.getLogger('base')
        logger.info('Path already exists. Rename it to [{:s}]'.format(new_name))
        os.rename(path, new_name)
    os.makedirs(path)


def set_random_seed(seed):
    random.seed(seed)
    np.random.seed(seed)
    torch.manual_seed(seed)
    torch.cuda.manual_seed_all(seed)


def setup_logger(logger_name, root, phase, level=logging.INFO, screen=False, tofile=False):
    '''set up logger'''
    lg = logging.getLogger(logger_name)
    formatter = logging.Formatter('%(asctime)s.%(msecs)03d - %(levelname)s: %(message)s',
                                  datefmt='%y-%m-%d %H:%M:%S')
    lg.setLevel(level)
    if tofile:
        log_file = os.path.join(root, phase + '_{}.log'.format(get_timestamp()))
        fh = logging.FileHandler(log_file, mode='w')
        fh.setFormatter(formatter)
        lg.addHandler(fh)
    if screen:
        sh = logging.StreamHandler()
        sh.setFormatter(formatter)
        lg.addHandler(sh)


####################
# image convert
####################


def tensor2img(tensor, out_type=np.uint8, min_max=(0, 1)):
    '''
    Converts a torch Tensor into an image Numpy array
    Input: 4D(B,(3/1),H,W), 3D(C,H,W), or 2D(H,W), any range, RGB channel order
    Output: 3D(H,W,C) or 2D(H,W), [0,255], np.uint8 (default)
    '''
    tensor = tensor.squeeze().float().cpu().clamp_(*min_max)  # clamp
    tensor = (tensor - min_max[0]) / (min_max[1] - min_max[0])  # to range [0,1]
    n_dim = tensor.dim()
    if n_dim == 4:
        n_img = len(tensor)
        img_np = make_grid(tensor, nrow=int(math.sqrt(n_img)), normalize=False).numpy()
        img_np = np.transpose(img_np[[2, 1, 0], :, :], (1, 2, 0))  # HWC, BGR
    elif n_dim == 3:
        img_np = tensor.numpy()
        img_np = np.transpose(img_np[[2, 1, 0], :, :], (1, 2, 0))  # HWC, BGR
    elif n_dim == 2:
        img_np = tensor.numpy()
    else:
        raise TypeError(
            'Only support 4D, 3D and 2D tensor. But received with dimension: {:d}'.format(n_dim))
    if out_type == np.uint8:
        img_np = (img_np * 255.0).round()
        # Important. Unlike matlab, numpy.unit8() WILL NOT round by default.
    return img_np.astype(out_type)


def save_img(img, img_path, mode='RGB'):
    cv2.imwrite(img_path, img)

####################
# metric
####################
def calculate_rmse(img1, img2):
    # img1 and img2 have range [0, 255]
    img1 = img1.astype(np.float64)
    img2 = img2.astype(np.float64)
    rmse = np.sqrt(np.mean((img1 - img2)**2))
    return rmse

def calculate_psnr(img1, img2):
    # img1 and img2 have range [0, 255]
    img1 = img1.astype(np.float64)
    img2 = img2.astype(np.float64)
    mse = np.mean((img1 - img2)**2)
    if mse == 0:
        return float('inf')
    return 20 * math.log10(255.0 / math.sqrt(mse))
# def calculate_psnr_fastmri(img1, img2):
#     img1 = img1.astype(np.float64)
#     img2 = img2.astype(np.float64)
#     # print(img1.min(), img1.max())
#     img1 = (img1-img1.min())/(img1.max()-img1.min())
#     img2 = (img2-img2.min())/(img2.max()-img2.min())
#     mse = np.mean((img1 - img2)**2)
#     if mse == 0:
#         return float('inf')
#     return 20 * math.log10(img2.max() / math.sqrt(mse))
def calculate_psnr_fastmri(gt, pred):
    """Compute Peak Signal to Noise Ratio metric (PSNR)"""
    return peak_signal_noise_ratio(gt, pred, data_range=gt.max())

def ssim(img1, img2):
    C1 = (0.01 * 255)**2
    C2 = (0.03 * 255)**2

    img1 = img1.astype(np.float64)
    img2 = img2.astype(np.float64)
    kernel = cv2.getGaussianKernel(11, 1.5)
    window = np.outer(kernel, kernel.transpose())

    mu1 = cv2.filter2D(img1, -1, window)[5:-5, 5:-5]  # valid
    mu2 = cv2.filter2D(img2, -1, window)[5:-5, 5:-5]
    mu1_sq = mu1**2
    mu2_sq = mu2**2
    mu1_mu2 = mu1 * mu2
    sigma1_sq = cv2.filter2D(img1**2, -1, window)[5:-5, 5:-5] - mu1_sq
    sigma2_sq = cv2.filter2D(img2**2, -1, window)[5:-5, 5:-5] - mu2_sq
    sigma12 = cv2.filter2D(img1 * img2, -1, window)[5:-5, 5:-5] - mu1_mu2

    ssim_map = ((2 * mu1_mu2 + C1) * (2 * sigma12 + C2)) / ((mu1_sq + mu2_sq + C1) *
                                                            (sigma1_sq + sigma2_sq + C2))
    return ssim_map.mean()


def calculate_ssim(img1, img2):
    '''calculate SSIM
    the same outputs as MATLAB's
    img1, img2: [0, 255]
    '''
    if not img1.shape == img2.shape:
        raise ValueError('Input images must have the same dimensions.')
    if img1.ndim == 2:
        return ssim(img1, img2)
    elif img1.ndim == 3:
        if img1.shape[2] == 3:
            ssims = []
            for i in range(3):
                ssims.append(ssim(img1, img2))
            return np.array(ssims).mean()
        elif img1.shape[2] == 1:
            return ssim(np.squeeze(img1), np.squeeze(img2))
    else:
        raise ValueError('Wrong input image dimensions.')


class ProgressBar(object):
    '''A progress bar which can print the progress
    modified from https://github.com/hellock/cvbase/blob/master/cvbase/progress.py
    '''

    def __init__(self, task_num=0, bar_width=50, start=True):
        self.task_num = task_num
        max_bar_width = self._get_max_bar_width()
        self.bar_width = (bar_width if bar_width <= max_bar_width else max_bar_width)
        self.completed = 0
        if start:
            self.start()

    def _get_max_bar_width(self):
        terminal_width, _ = get_terminal_size()
        max_bar_width = min(int(terminal_width * 0.6), terminal_width - 50)
        if max_bar_width < 10:
            print('terminal width is too small ({}), please consider widen the terminal for better '
                  'progressbar visualization'.format(terminal_width))
            max_bar_width = 10
        return max_bar_width

    def start(self):
        if self.task_num > 0:
            sys.stdout.write('[{}] 0/{}, elapsed: 0s, ETA:\n{}\n'.format(
                ' ' * self.bar_width, self.task_num, 'Start...'))
        else:
            sys.stdout.write('completed: 0, elapsed: 0s')
        sys.stdout.flush()
        self.start_time = time.time()

    def update(self, msg='In progress...'):
        self.completed += 1
        elapsed = time.time() - self.start_time
        fps = self.completed / elapsed
        if self.task_num > 0:
            percentage = self.completed / float(self.task_num)
            eta = int(elapsed * (1 - percentage) / percentage + 0.5)
            mark_width = int(self.bar_width * percentage)
            bar_chars = '>' * mark_width + '-' * (self.bar_width - mark_width)
            sys.stdout.write('\033[2F')  # cursor up 2 lines
            sys.stdout.write('\033[J')  # clean the output (remove extra chars since last display)
            sys.stdout.write('[{}] {}/{}, {:.1f} task/s, elapsed: {}s, ETA: {:5}s\n{}\n'.format(
                bar_chars, self.completed, self.task_num, fps, int(elapsed + 0.5), eta, msg))
        else:
            sys.stdout.write('completed: {}, elapsed: {}s, {:.1f} tasks/s'.format(
                self.completed, int(elapsed + 0.5), fps))
        sys.stdout.flush()

####################
# read image
####################

def read_image(img_path):
    '''read one image from img_path
    Return img: HWC, BGR, [0,1], numpy
    '''
    img_GT = cv2.imread(img_path)
    img = img_GT.astype(np.float32) / 255.

    return img

def read_seq_imgs(img_seq_path):
    '''read a sequence of images'''
    img_path_l = glob.glob(img_seq_path + '/*')
    # img_path_l.sort(key=lambda x: int(os.path.basename(x)[:-4]))
    img_path_l.sort(key=lambda x: int(re.search(r'\d+', os.path.basename(x)).group()))
    img_l = [read_image(v) for v in img_path_l]
    # stack to TCHW, RGB, [0,1], torch
    imgs = np.stack(img_l, axis=0)
    imgs = imgs[:, :, :, [2, 1, 0]]
    imgs = torch.from_numpy(np.ascontiguousarray(np.transpose(imgs, (0, 3, 1, 2)))).float()
    return imgs

def test_index_generation_multiple_frames(skip, N_out, len_in, use_topAndEnd=None):
    # number of input frames for the network
    N_in = 1 + N_out // 2
    # input length should be enough to generate the output frames
    assert N_in <= len_in

    sele_list = []
    if skip: 
        right = N_out # init
        while (right <= len_in):
            h_list = [right-N_out+x for x in range(N_out)]
            if use_topAndEnd == True:
                l_list = []
                l_list.append(h_list[0])
                l_list.append(h_list[-1])
            else:
                l_list = h_list[::2]
            right += (N_out - 1)
            sele_list.append([l_list,h_list])
    else:
        right = N_out # init
        right_in = N_in
        while (right_in <= len_in):
            h_list = [right-N_out+x for x in range(N_out)]
            l_list = [right_in-N_in+x for x in range(N_in)]
            right += (N_out - 1)
            right_in += (N_in - 1)
            sele_list.append([l_list,h_list])
    # check if it covers the last image, if not, we should cover it 
    if (skip) and (right != len_in - 1):
        h_list = [len_in - N_out + x for x in range(N_out)]
        if use_topAndEnd == True:
            l_list = []
            l_list.append(h_list[0])
            l_list.append(h_list[-1])
        else:
            l_list = h_list[::2]
        sele_list.append([l_list,h_list])
    if (not skip) and (right_in != len_in - 1):
        right = len_in * 2 - 1;
        h_list = [right-N_out+x for x in range(N_out)]
        l_list = [len_in - N_in + x for x in range(N_in)]
        sele_list.append([l_list,h_list])        
    return sele_list

def test_index_generation_demo(skip, N_out, len_in):
    N_in = 1 + N_out // 2
    assert N_in <= len_in
    sele_list = []
    right = 2 # right = N_out
    while (right < len_in):
        h_list = []
        l_list = []
        
        l_list.append(right - 2)
        l_list.append(right)
        
        for i in range(2):
            for j in range(int((N_out - 1) / 2)):
                h_list.append(str(int(right - 2 + i)) + '-' + str(int(j)) + '-' + str(int(right - 1 + i)))
        h_list.append(str(int(right)) + '-' + str(int(0)) + '-' + str(int(right + 1)))
        right += 2
        sele_list.append([l_list, h_list])

    return sele_list

def test_index_generation(skip, N_out, len_in):
    '''
    params: 
    skip: if skip even number; 
    N_out: number of frames of the network; 
    len_in: length of input frames

    example:
  len_in | N_out  | times | (no skip)                  |   (skip)
    5    |   3    |  4/2  | [0,1], [1,2], [2,3], [3,4] | [0,2],[2,4]
    7    |   3    |  5/3  | [0,1],[1,2][2,3]...[5,6]   | [0,2],[2,4],[4,6] 
    5    |   5    |  2/1  | [0,1,2] [2,3,4]            | [0,2,4]
    '''
    # number of input frames for the network
    N_in = 1 + N_out // 2
    # input length should be enough to generate the output frames
    assert N_in <= len_in

    sele_list = []
    if skip: 
        right = N_out # init
        while (right <= len_in):
            h_list = [right-N_out+x for x in range(N_out)]
            l_list = h_list[::2]
            right += (N_out - 1)
            sele_list.append([l_list,h_list])
    else:
        right = N_out # init
        right_in = N_in
        while (right_in <= len_in):
            h_list = [right-N_out+x for x in range(N_out)]
            l_list = [right_in-N_in+x for x in range(N_in)]
            right += (N_out - 1)
            right_in += (N_in - 1)
            sele_list.append([l_list,h_list])
    # check if it covers the last image, if not, we should cover it 
    if (skip) and (right < len_in - 1):
        h_list = [len_in - N_out + x for x in range(N_out)]
        l_list = h_list[::2]
        sele_list.append([l_list,h_list])
    if (not skip) and (right_in < len_in - 1):
        right = len_in * 2 - 1;
        h_list = [right-N_out+x for x in range(N_out)]
        l_list = [len_in - N_in + x for x in range(N_in)]
        sele_list.append([l_list,h_list])        
    return sele_list


####################
# video
####################

def extract_frames(ffmpeg_dir, video, outDir):
    """
    Converts the `video` to images.
    Parameters
    ----------
        video : string
            full path to the video file.
        outDir : string
            path to directory to output the extracted images.
    Returns
    -------
        error : string
            Error message if error occurs otherwise blank string.
    """

    error = ""
    print('{} -i {} -vsync 0 {}/%06d.png'.format(os.path.join(ffmpeg_dir, "ffmpeg"), video, outDir))
    retn = os.system('{} -i "{}" -vsync 0 {}/%06d.png'.format(os.path.join(ffmpeg_dir, "ffmpeg"), video, outDir))
    if retn:
        error = "Error converting file:{}. Exiting.".format(video)
    return error

def create_video(ffmpeg_dir, dir, output, fps):
    error = ""
    # print('{} -r {} -i {}/%6d.png -vcodec ffvhuff {}'.format(os.path.join(ffmpeg_dir, "ffmpeg"), fps, dir, output))
    # retn = os.system('{} -r {} -i {}/%6d.png -vcodec ffvhuff "{}"'.format(os.path.join(ffmpeg_dir, "ffmpeg"), fps, dir, output))
    print('{} -r {} -f image2 -i {}/%6d.png {}'.format(os.path.join(ffmpeg_dir, "ffmpeg"), fps, dir, output))
    retn = os.system('{} -r {} -f image2 -i {}/%6d.png {}'.format(os.path.join(ffmpeg_dir, "ffmpeg"), fps, dir, output))
    if retn:
        error = "Error creating output video. Exiting."
    return error

# if __name__ == '__main__':

utils/util.py
lpcccc-cv/MC-VarNet