MRI_recon/code/Frequency-Diffusion/dataset/basic.py

# Dataloader for abdominal images
import glob
import numpy as np
from .m4_utils import niftiio as nio
from .m4_utils import transform_utils as trans
from .m4_utils.abd_dataset_utils import get_normalize_op
from .m4_utils.transform_albu import get_albu_transforms, get_resize_transforms
import copy
import random, cv2, os
import torch.utils.data as torch_data
import math
import itertools
from pdb import set_trace
from multiprocessing import Process
import albumentations as A
from tqdm import tqdm


def get_basedir(data_dir):
    return os.path.join(data_dir, "Abdominal")


class BasicDataset(torch_data.Dataset):
    def __init__(self, fineSize, mode, transforms, base_dir, domains: list, aux_modality, pseudo = False, 
                 idx_pct = [0.7, 0.1, 0.2], tile_z_dim = 3, extern_norm_fn = None, 
                 LABEL_NAME=["bg", "fore"], debug=False, nclass=4, num_channels=3,
                 filter_non_labeled=False, use_diff_axis_view=False, chunksize=200):
        """
        Args:
            mode:               'train', 'val', 'test', 'test_all'
            transforms:         naive data augmentations used by default. Photometric transformations slightly better than those configured by Zhang et al. (bigaug)
            idx_pct:            train-val-test split for source domain
            extern_norm_fn:     feeding data normalization functions from external, only concerns CT-MR cross domain scenario
        """
        super(BasicDataset, self).__init__()

        self.fineSize = fineSize
        self.transforms = transforms
        self.nclass = nclass
        self.debug = debug
        self.is_train = True if mode == 'train' else False
        self.phase = mode
        self.domains = domains
        self.num_channels = num_channels
        
        # Modality
        self.main_modality = domains[-1].split("-")[-1]
        self.aux_modality = aux_modality.upper()
        
        print(f"=== Donmain: {domains}, Main modality: {self.main_modality}, Aux modality: {self.aux_modality}")
        
        self.pseudo = pseudo
        self.all_label_names = LABEL_NAME
        self.nclass = len(LABEL_NAME)
        self.tile_z_dim = tile_z_dim
        self._base_dir = base_dir
        self.idx_pct = idx_pct
        # self.albu_transform = get_albu_transforms((fineSize, fineSize))
        self.test_resizer   = get_resize_transforms(fineSize)
        self.fake_interpolate    = True # True
        self.use_diff_axis_view  = use_diff_axis_view
        self.filter_non_labeled  = filter_non_labeled
        self.input_window = 1 

        self.resizer = A.Compose([
            A.Resize(fineSize[0], fineSize[1], interpolation=cv2.INTER_NEAREST)
        ], p=1.0, additional_targets={'image2': 'image', "mask2": "mask"})

        self.img_pids = {} 
        for _domain in self.domains: # load file names
            if "BraTS" in _domain:
                self.img_pids[_domain] = sorted([ fid.split("-")[-2] for fid in
                                     glob.glob(self._base_dir + "/" + _domain + "/img/*.nii.gz") ],
                                     key = lambda x: int(x))
                
            else:
                self.img_pids[_domain] = sorted([fid.split("_")[-1].split(".nii.gz")[0] for fid in
                                                 glob.glob(self._base_dir + "/" + _domain + "/img/*.nii.gz")],
                                                key=lambda x: int(x))
                
        self.scan_ids = self.__get_scanids(mode, idx_pct) # train val test split in terms of patient ids
        try:
            print(f'For {self.phase} on {[_dm for _dm in self.domains]} using scan ids len = ' + \
              f'{[len(self.scan_ids[_dm]) for _dm in self.scan_ids.keys()]}')
        except:
            print("Errors of self.scan_ids")
            print(self.scan_ids)


        self.info_by_scan = None
        self.sample_list = self.__search_samples(self.scan_ids) # image files names according to self.scan_ids
        if self.is_train:
            
            self.pid_curr_load = self.scan_ids
        elif mode == 'val':
            self.pid_curr_load = self.scan_ids
        elif mode == 'test': # Source domain test
            self.pid_curr_load = self.scan_ids
        elif mode == 'test_all':
            # Choose this when being used as a target domain testing set. Liu et al.
            self.pid_curr_load = self.scan_ids
            
        if extern_norm_fn is None:
            self.normalize_op = get_normalize_op(self.domains[0], [itm['img_fid'] for _, itm in
                                                                   self.sample_list[self.domains[0]].items() ])
            print(f'{self.phase}_{self.domains[0]}: Using fold data statistics for normalization')
            
        else:
            # assert len(self.domains) == 1, 'for now we only support one normalization function for the entire set'
            self.normalize_op = extern_norm_fn


        # load to memory
        # self.sample_list All
        self.actual_dataset = None
        self.chunksize  = chunksize if not debug else 3
        
        self.chunk_id = 0
        self.chunk_pool, self.current_chunk = {}, {}
        for _domain, item in self.sample_list.items():
            self.chunk_pool[_domain] = list(item.keys())

        chunk, status = self.next_chunk(self.sample_list)
        self.actual_dataset = self.__read_dataset(chunk, status)
        self.size = len(self.actual_dataset)     # 2D
        
        print("----- Set up dataset for", self.phase, "with chunksize=", chunksize)

    def update_chunk(self):
        chunk, status = self.next_chunk(self.sample_list)
        self.actual_dataset = self.__read_dataset(chunk, status)

    def __get_scanids(self, mode, idx_pct):
        """
        index by domains given that we might need to load multi-domain data
        idx_pct: [0.7 0.1 0.2] for train val test. with order te val tr
        """
        tr_ids      = {}
        val_ids     = {}
        te_ids      = {}
        te_all_ids  = {}

        for _domain in self.domains:
            dset_size   = len(self.img_pids[_domain])
            tr_size     = round(dset_size * idx_pct[0])
            val_size    = math.floor(dset_size * idx_pct[1])
            te_size     = dset_size - tr_size - val_size
            # print('te_size = ', te_size)

            te_ids[_domain]     = self.img_pids[_domain][: te_size]
            val_ids[_domain]    = self.img_pids[_domain][te_size: te_size + val_size]
            tr_ids[_domain]     = self.img_pids[_domain][te_size + val_size: ]
            te_all_ids[_domain] = list(itertools.chain(tr_ids[_domain], te_ids[_domain], val_ids[_domain]   ))

        print(" self.phase = ", self.phase)
        if self.phase == 'train':
            return tr_ids
        elif self.phase == 'val':
            return val_ids
        elif self.phase == 'test':
            return te_ids
        elif self.phase == 'test_all':
            return te_all_ids

    def __search_samples(self, scan_ids):
        """search for filenames for images and masks
        """
        out_list = {}
        for _domain, id_list in scan_ids.items():
            domain_dir = os.path.join(self._base_dir, _domain)
            print("=== reading domains from:", domain_dir)
            out_list[_domain] = {}
            for curr_id in id_list:
                curr_dict = {}
                if "BraTS" in _domain:

                    _img_fid = os.path.join(domain_dir, 'img', f'{_domain[:-4]}-{curr_id}-000.nii.gz')
                    if not self.pseudo:
                        _lb_fid  = os.path.join(domain_dir, 'seg', f'{_domain[:-4]}-{curr_id}-000.nii.gz')
                    else:
                        _lb_fid  = os.path.join(domain_dir, 'seg', f'{_domain[:-4]}-{curr_id}-000.nii.gz.npy')  # npy

                    _aux_fid = _img_fid.replace(self.main_modality, self.aux_modality)

                

                curr_dict["img_fid"] = _img_fid
                curr_dict["lbs_fid"] = _lb_fid
                curr_dict["aux_fid"] = _aux_fid
                out_list[_domain][str(curr_id)] = curr_dict

        print("=== search sample num:", len(out_list))
        return out_list
    
    
    def filter_with_label(self, img, lb, aux):
        # H, W, C, filter zero
        if self.phase == "train":
        
            filter = np.any(np.any(img, axis=0), axis=0)
            img, lb, aux = img[..., filter], lb[..., filter], aux[..., filter]

            
            if self.filter_non_labeled:
            
                if self.dataset_key == "knee":
                    filter2 = np.any(np.any(lb == 2, axis=0), axis=0)
                    filter4 = np.any(np.any(lb == 4, axis=0), axis=0)
                    filter = filter2 + filter4
                else:
                    filter = np.any(np.any(lb, axis=0), axis=0)
                
                filter_right = np.roll(filter, 3)
                filter_left  = np.roll(filter, -3)
                filter = filter + filter_right + filter_left
                filter = filter > 0
                
                # HWC
                img, lb, aux = img[..., filter], lb[..., filter], aux[..., filter]
                
        return img, lb, aux

    def __read_dataset(self, chunk, status):
        """
        Read the dataset into memory
        """

        out_list = []
        self.info_by_scan = {} # meta data of each scan
        glb_idx = 0 # global index of a certain slice in a certain scan in entire dataset
        for _domain, _curr_chunk in tqdm(chunk.items()):  # .items()
            domain_ids = 0
            if status[_domain] != 3:
                print(f"==== UPDATE dataset for: {_domain} w/ status = {status[_domain]}")
            
            for scan_id in _curr_chunk:
                domain_ids += 1
                if domain_ids > self.chunksize:
                    print(f"=== UPDATE finished")
                    break
            
                itm = self.sample_list[_domain][scan_id]
                if scan_id not in self.pid_curr_load[_domain]:
                    continue 
                
                # Keep the original dataset
                if (status[_domain] == 0) or (status[_domain] == 2 and domain_ids <= self.chunksize // 2):
                    size = self.actual_dataset[glb_idx]['size']
                    out_list.extend(self.actual_dataset[glb_idx: glb_idx + size])  # Original dataset
                    glb_idx += size
                    continue
                
                if (status[_domain] == 1 and domain_ids > self.chunksize // 2):
                    try:
                        size = self.actual_dataset[glb_idx]['size']
                        out_list.extend(self.actual_dataset[glb_idx: glb_idx + size])  # Original dataset
                        glb_idx += size
                        continue
                    except:
                        print(f"=== Warning (domain_ids={domain_ids}) getting glb_idx={glb_idx} from actual_dataset length={len(self.actual_dataset)}") 
                        # print(self.actual_dataset)
                        
                img, _info = nio.read_nii_bysitk(itm["img_fid"], peel_info = True) # get the meta information out
                self.info_by_scan[_domain + '_' + scan_id] = _info

                img_original = np.float32(img)
                img = img_original.copy()
                
                aux = nio.read_nii_bysitk(itm["aux_fid"])
                aux_original = np.float32(aux)
                aux = aux_original.copy()
                

                # img, self.mean, self.std = self.normalize_op(img)
                _, mean, std = self.normalize_op(img)
                _, aux_mean, aux_std = self.normalize_op(aux)

                if not self.pseudo:
                    lb = nio.read_nii_bysitk(itm["lbs_fid"])
                else:
                    uncertainty_thr = 0.05  #  0.05
                    lb_cache = np.load(itm["lbs_fid"], allow_pickle=True).item()
                    lb = lb_cache['pseudo'].cpu().numpy()              # "pseudo": curr_pred, "score":curr_score , "uncertainty"
                    uncertainty = lb_cache['uncertainty'].cpu().numpy()   # Z, C, H, W
                    uncertainty = np.float32(uncertainty)
                    
                    new_lb = np.zeros_like(lb)
                    for cls in range(self.nclass - 1):
                        un_mask = (uncertainty[:, cls+1] < uncertainty_thr ) * (cls+1)
                        new_lb[lb == (cls+1)] = un_mask[lb == (cls+1)]

                    lb = new_lb

                lb_original = np.float32(lb)
                lb = lb_original.copy()

                # -> H, W, C
                img, lb, aux = map(lambda arr: np.transpose(arr, (1, 2, 0)), [img, lb, aux])
                assert img.shape[-1] == lb.shape[-1], f"ASSERT {img.shape} = {lb.shape}"

                # Resize:
                if img.shape[1] != self.fineSize[1]:
                    # H, W, C
                    res = self.resizer(image=img, mask=lb, image2=aux)
                    img, lb, aux = res['image'], res['mask'], res['image2']

                prt_cache = f" {_domain} stat ({domain_ids}/{len(_curr_chunk)}): shape={img.shape}, max={img.max()}, min={img.min()}"

                # Filter vacant slices
                if self.phase == "train":
                    filter = np.any(np.any(img, axis=0), axis=0)
                    img, lb, aux = img[..., filter], lb[..., filter], aux[..., filter]
            
                img, lb, aux = self.filter_with_label(img, lb, aux)

                out_list, glb_idx = self.add_to_list(glb_idx, out_list, img, lb,
                                                     aux, mean, aux_mean, aux_std, std, _domain,
                                                     scan_id, itm["img_fid"])

                if (domain_ids) % (len(_curr_chunk) // 2) == 0:
                    print(prt_cache + f", filtered shape={img.shape}, mask max={lb.max()}")


                # Add various axis view !!!
                if self.phase == "train" and self.use_diff_axis_view:
                    # C, W, H
                    img, lb, aux = img_original, lb_original, aux_original
                    # Resize:
                    if img.shape[1] != self.fineSize[1]:
                        res = self.resizer(image=img, mask=lb, image2=aux) # assume H, W, (C)<-
                        img, lb, aux = res['image'], res['mask'], res['image2']
                    
                    img, lb, aux = self.filter_with_label(img, lb, aux)
                    
                    out_list, glb_idx = self.add_to_list(glb_idx, out_list, img, lb,
                                                     aux, mean, aux_mean, aux_std, std, _domain,
                                                     scan_id, itm["img_fid"])

                del img, lb, aux, img_original, lb_original, aux_original

        del self.actual_dataset
        return out_list
    
    def next_chunk(self, all_samples):
        # 0 No update, 1 First half, 2 Second half, 3 All updates Chunk
        status = {}
        self.last_chunk = copy.deepcopy(self.current_chunk)
        for _domain, _sample_list in tqdm(all_samples.items()):
            # Default value
            status[_domain] = 3 
            
            # Put all in - validation or small dataset
            if ((not self.is_train) or len(_sample_list) < self.chunksize) and not self.debug:
                self.current_chunk[_domain] = _sample_list
                if _domain not in self.last_chunk:
                    status[_domain] = 3  # all
                else:
                    status[_domain] = 0  # not updates
                print("=== Put all data in for", _domain)
                continue
            
            # chunksize
            random.shuffle(self.chunk_pool[_domain])
            if _domain not in self.last_chunk:
                status[_domain] = 3
                self.current_chunk[_domain] = self.chunk_pool[_domain][:self.chunksize]
                self.chunk_pool[_domain]    = self.chunk_pool[_domain][self.chunksize:]

            else:
                status[_domain] = self.chunk_id//2 + 1  # 1, 2
                candidate                = self.chunk_pool[_domain][:self.chunksize//2]
                self.chunk_pool[_domain] = self.chunk_pool[_domain][self.chunksize //2:]
                if status[_domain] == 1:
                    self.current_chunk[_domain][:self.chunksize // 2] = candidate
                else:
                    self.current_chunk[_domain][self.chunksize  // 2:] = candidate

            if _domain in self.last_chunk:
                self.chunk_pool[_domain] = self.chunk_pool[_domain] + self.last_chunk[_domain]
        
        self.chunk_id += 1

        return self.current_chunk, status


    def add_to_list(self, glb_idx, out_list, img, lb, aux, mean, std, aux_mean, aux_std, _domain, scan_id, file_id):
        # now start writing everthing in
        c = 3
        
        for ii in range(img.shape[-1]):
            is_end   = False
            is_start   = False
            if ii == 0:
                is_start = True
                # write the beginning frame
                if self.input_window == 3:
                    _img = img[..., 0: c].copy()
                    _img[..., 1] = _img[..., 0]
                elif self.input_window == 1:
                    _img = img[..., 0: 0 + 1].copy()
               
                
            elif ii < img.shape[-1] - 1:
                if self.input_window == 3:
                    _img = img[..., ii -1: ii + 2].copy()
                elif self.input_window == 1:
                    _img = img[..., ii: ii + 1].copy()
                    
            else:
                is_end = True
                if self.input_window == 3:
                    _img = img[..., ii-2: ii + 1].copy()
                    _img[..., 0] = _img[..., 1]
                elif self.input_window == 1:
                    _img = img[..., ii: ii+ 1].copy()
            
            _lb  = lb[..., ii: ii + 1]
            _aux = aux[..., ii: ii + 1]
            
            out_list.append(
                       {"img": _img, "lb":_lb, "aux":_aux, "size": img.shape[-1],
                        "mean":mean, "std":std,
                        "aux_mean": aux_mean, "aux_std": aux_std,
                        "is_start": is_start, "is_end": is_end,
                        "domain": _domain, "nframe": img.shape[-1],
                        "scan_id": _domain + "_" + scan_id,
                        "pid": scan_id, "file_id": file_id, "z_id":ii})
            glb_idx += 1
        
        return out_list, glb_idx


    def get_patch_from_img(self, img_H, img_L, img_L2, crop_size=[320, 320], zslice_dim=2):
        # --------------------------------
        # randomly crop the patch
        # --------------------------------

        H, W, _ = img_H.shape
        rnd_h = random.randint(0, max(0, H - crop_size[0]))
        rnd_w = random.randint(0, max(0, W - crop_size[1]))

        # image = torch.index_select(image, 0, torch.tensor([1]))
        if zslice_dim == 2:
            patch_H = img_H[rnd_h:rnd_h + crop_size[0], rnd_w:rnd_w + crop_size[1], :]
            patch_L = img_L[rnd_h:rnd_h + crop_size[0], rnd_w:rnd_w + crop_size[1], :]
            patch_L2 = img_L2[rnd_h:rnd_h + crop_size[0], rnd_w:rnd_w + crop_size[1], :]
        elif zslice_dim == 0:
            patch_H = img_H[:, rnd_h:rnd_h + crop_size[0], rnd_w:rnd_w + crop_size[1]]
            patch_L = img_L[:, rnd_h:rnd_h + crop_size[0], rnd_w:rnd_w + crop_size[1]]
            patch_L2 = img_L2[:, rnd_h:rnd_h + crop_size[0], rnd_w:rnd_w + crop_size[1]]

        return patch_H, patch_L, patch_L2


    def __len__(self):
        """
        copy-paste from basic naive dataset configuration
        """
        return len(self.actual_dataset)