datasets.py

"""
datasets.py
=======================
Houses the DynamicImageDataset class, also functions to help with image color channel normalization, transformers, etc..
"""

import torch
from torchvision import transforms
import os
import dask
#from dask.distributed import Client; Client()
import dask.array as da, pandas as pd, numpy as np
from pathflowai.utils import *
import pysnooper
import nonechucks as nc
from torch.utils.data import Dataset, DataLoader
import random
import albumentations as alb
import copy
from albumentations import pytorch as albtorch
from sklearn.preprocessing import LabelBinarizer
from sklearn.utils.class_weight import compute_class_weight
from pathflowai.losses import class2one_hot
import cv2
from scipy.ndimage.morphology import generate_binary_structure
from dask_image.ndmorph import binary_dilation
cv2.setNumThreads(0)
cv2.ocl.setUseOpenCL(False)


def RandomRotate90():
	"""Transformer for random 90 degree rotation image.

	Returns
	-------
	function
		Transformer function for operation.

	"""
	return (lambda img: img.rotate(random.sample([0, 90, 180, 270], k=1)[0]))

def get_data_transforms(patch_size = None, mean=[], std=[], resize=False, transform_platform='torch', elastic=True, user_transforms=dict()):
	"""Get data transformers for training test and validation sets.

	Parameters
	----------
	patch_size:int
		Original patch size being transformed.
	mean:list of float
		Mean RGB
	std:list of float
		Std RGB
	resize:int
		Which patch size to resize to.
	transform_platform:str
		Use pytorch or albumentation transforms.
	elastic:bool
		Whether to add elastic deformations from albumentations.

	Returns
	-------
	dict
		Transformers.

	"""
	transform_dict=dict(torch=dict(
									colorjitter=lambda kargs: transforms.ColorJitter(**kargs),
									hflip=lambda kargs: transforms.RandomHorizontalFlip(),
									vflip=lambda kargs: transforms.RandomVerticalFlip(),
									r90= lambda kargs: RandomRotate90()
									),
						albumentations=dict(
							huesaturation=lambda kargs: alb.augmentations.transforms.HueSaturationValue(**kargs),
							flip=lambda kargs: alb.augmentations.transforms.Flip(**kargs),
							transpose=lambda kargs: alb.augmentations.transforms.Transpose(**kargs),
							affine=lambda kargs: alb.augmentations.transforms.ShiftScaleRotate(**kargs),
							r90=lambda kargs: alb.augmentations.transforms.RandomRotate90(**kargs),
							elastic=lambda kargs: alb.augmentations.transforms.ElasticTransform(**kargs)
						))
	if 'normalization' in user_transforms:
		mean=user_transforms['normalization'].pop('mean')
		std=user_transforms['normalization'].pop('std')
		del user_transforms['normalization']
	default_transforms=dict() # add normalization custom
	default_transforms['torch']=dict(
							colorjitter=dict(brightness=0.8, contrast=0.8, saturation=0.8, hue=0.5),
							hflip=dict(),
							vflip=dict(),
							r90=dict())
	default_transforms['albumentations']=dict(
							huesaturation=dict(hue_shift_limit=20, sat_shift_limit=30, val_shift_limit=20, p=0.5),
							r90=dict(p=0.5),
							elastic=dict(p=0.5))
	main_transforms = default_transforms[transform_platform] if not user_transforms else user_transforms
	print(main_transforms)
	train_transforms=[transform_dict[transform_platform][k](v) for k,v in main_transforms.items()]
	torch_init=[transforms.ToPILImage(),transforms.Resize((patch_size,patch_size)),transforms.CenterCrop(patch_size)]
	albu_init=[alb.augmentations.transforms.Resize(patch_size, patch_size),
				alb.augmentations.transforms.CenterCrop(patch_size, patch_size)]
	tensor_norm=[transforms.ToTensor(),transforms.Normalize(mean if mean else [0.7, 0.6, 0.7], std if std is not None else [0.15, 0.15, 0.15])] #mean and standard deviations for lung adenocarcinoma resection slides
	data_transforms = { 'torch': {
		'train': transforms.Compose(torch_init+train_transforms+tensor_norm),
		'val': transforms.Compose([
			transforms.ToPILImage(),
			transforms.Resize((patch_size,patch_size)),
			transforms.CenterCrop(patch_size),
			transforms.ToTensor(),
			transforms.Normalize(mean if mean else [0.7, 0.6, 0.7], std if std is not None else [0.15, 0.15, 0.15])
		]),
		'test': transforms.Compose([
			transforms.ToPILImage(),
			transforms.Resize((patch_size,patch_size)),
			transforms.CenterCrop(patch_size),
			transforms.ToTensor(),
			transforms.Normalize(mean if mean else [0.7, 0.6, 0.7], std if std is not None else [0.15, 0.15, 0.15])
		]),
		'pass': transforms.Compose([
			transforms.ToPILImage(),
			transforms.CenterCrop(patch_size),
			transforms.ToTensor(),
		])
	},
	'albumentations':{
	'train':alb.core.composition.Compose(albu_init+train_transforms),
	'val':alb.core.composition.Compose([
		alb.augmentations.transforms.Resize(patch_size, patch_size),
		alb.augmentations.transforms.CenterCrop(patch_size, patch_size)
	]),
	'test':alb.core.composition.Compose([
		alb.augmentations.transforms.Resize(patch_size, patch_size),
		alb.augmentations.transforms.CenterCrop(patch_size, patch_size)
	]),
	'normalize':transforms.Compose([transforms.Normalize(mean if mean else [0.7, 0.6, 0.7], std if std is not None else [0.15, 0.15, 0.15])])
	}}

	return data_transforms[transform_platform]

def create_transforms(mean, std):
	"""Create transformers.

	Parameters
	----------
	mean:list
		See get_data_transforms.
	std:list
		See get_data_transforms.

	Returns
	-------
	dict
		Transformers.

	"""
	return get_data_transforms(patch_size = 224, mean=mean, std=std, resize=True)



def get_normalizer(normalization_file, dataset_opts):
	"""Find mean and standard deviation of images in batches.

	Parameters
	----------
	normalization_file:str
		File to store normalization information.
	dataset_opts:type
		Dictionary storing information to create DynamicDataset class.

	Returns
	-------
	dict
		Stores RGB mean, stdev.

	"""
	if os.path.exists(normalization_file):
		norm_dict = torch.load(normalization_file)
	else:
		norm_dict = {'normalization_file':normalization_file}

	if 'normalization_file' in norm_dict:

		transformers = get_data_transforms(patch_size = 224, mean=[], std=[], resize=True, transform_platform='torch')

		dataset_opts['transformers']=transformers
		#print(dict(pos_annotation_class=pos_annotation_class, segmentation=segmentation, patch_size=patch_size, fix_names=fix_names, other_annotations=other_annotations))

		dataset = DynamicImageDataset(**dataset_opts)#nc.SafeDataset(DynamicImageDataset(**dataset_opts))

		if dataset_opts['classify_annotations']:
			dataset.binarize_annotations()

		dataloader = DataLoader(dataset, batch_size=128, shuffle=True, num_workers=4)

		all_mean = torch.tensor([0.,0.,0.],dtype=torch.float)#[]

		all_std = torch.tensor([0.,0.,0.],dtype=torch.float)

		if torch.cuda.is_available():
			all_mean=all_mean.cuda()
			all_std=all_std.cuda()

		with torch.no_grad():
			for i,(X,_) in enumerate(dataloader): # x,3,224,224
				if torch.cuda.is_available():
					X=X.cuda()
				all_mean += torch.mean(X, (0,2,3))
				all_std += torch.std(X, (0,2,3))

		N=i+1

		all_mean /= float(N) #(np.array(all_mean).mean(axis=0)).tolist()
		all_std /= float(N) #(np.array(all_std).mean(axis=0)).tolist()

		all_mean = all_mean.detach().cpu().numpy().tolist()
		all_std = all_std.detach().cpu().numpy().tolist()

		torch.save(dict(mean=all_mean,std=all_std),norm_dict['normalization_file'])

		norm_dict = torch.load(norm_dict['normalization_file'])
	return norm_dict

def segmentation_transform(img,mask, transformer, normalizer, alb_reduction):
	"""Run albumentations and return an image and its segmentation mask.

	Parameters
	----------
	img:array
		Image as array
	mask:array
		Categorical pixel by pixel.
	transformer :
		Transformation object.

	Returns
	-------
	tuple arrays
		Image and mask array.

	"""
	res=transformer(True, image=img, mask=mask)
	#res_mask_shape = res['mask'].size()
	return normalizer(torch.tensor(np.transpose(res['image']/alb_reduction,axes=(2,0,1)),dtype=torch.float)).float(), torch.tensor(res['mask']).long()#.view(res_mask_shape[0],res_mask_shape[1],res_mask_shape[2])

class DilationJitter:
	def __init__(self, dilation_jitter=dict(), segmentation=True, train_set=False):
		if dilation_jitter and segmentation and train_set:
			self.run_jitter=True
			self.dilation_jitter=dilation_jitter
			self.struct=generate_binary_structure(2,1) #structure=self.struct,
		else:
			self.run_jitter=False


	def __call__(self, mask):
		if self.run_jitter:
			for k in self.dilation_jitter:
				amount_jitter=int(round(max(np.random.normal(self.dilation_jitter[k]['mean'],
														self.dilation_jitter[k]['std']),1)))
				#print((mask==k).compute())
				mask[binary_dilation(mask==k,structure=self.struct,iterations=amount_jitter)]=k

		return mask


class DynamicImageDataset(Dataset):
	"""Generate image dataset that accesses images and annotations via dask.

	Parameters
	----------
	dataset_df:dataframe
		Dataframe with WSI, which set it is in (train/test/val) and corresponding WSI labels if applicable.
	set:str
		Whether train, test, val or pass (normalization) set.
	patch_info_file:str
		SQL db with positional and annotation information on each slide.
	transformers:dict
		Contains transformers to apply on images.
	input_dir:str
		Directory where images comes from.
	target_names:list/str
		Names of initial targets, which may be modified.
	pos_annotation_class:str
		If selected and predicting on WSI, this class is labeled as a positive from the WSI, while the other classes are not.
	other_annotations:list
		Other annotations to consider from patch info db.
	segmentation:bool
		Conducting segmentation task?
	patch_size:int
		Patch size.
	fix_names:bool
		Whether to change the names of dataset_df.
	target_segmentation_class:list
		Now can be used for classification as well, matched with two below options, samples images only from this class. Can specify this and below two options multiple times.
	target_threshold:list
		Sampled only if above this threshold of occurence in the patches.
	oversampling_factor:list
		Over sample them at this amount.
	n_segmentation_classes:int
		Number classes to segment.
	gdl:bool
		Using generalized dice loss?
	mt_bce:bool
		For multi-target prediction tasks.
	classify_annotations:bool
		For classifying annotations.

	"""
	# when building transformers, need a resize patch size to make patches 224 by 224
	#@pysnooper.snoop('init_data.log')
	def __init__(self,dataset_df, set, patch_info_file, transformers, input_dir, target_names, pos_annotation_class, other_annotations=[], segmentation=False, patch_size=224, fix_names=True, target_segmentation_class=-1, target_threshold=0., oversampling_factor=1., n_segmentation_classes=4, gdl=False, mt_bce=False, classify_annotations=False, dilation_jitter=dict(), modify_patches=True):

		#print('check',classify_annotations)
		reduce_alb=True
		self.patch_size=patch_size
		self.input_dir = input_dir
		self.alb_reduction=255. if reduce_alb else 1.
		self.transformer=transformers[set]
		original_set = copy.deepcopy(set)
		if set=='pass':
			set='train'
		self.targets = target_names
		self.mt_bce=mt_bce
		self.set = set
		self.segmentation = segmentation
		self.alb_normalizer=None
		if 'normalize' in transformers:
			self.alb_normalizer = transformers['normalize']
		if len(self.targets)==1:
			self.targets = self.targets[0]
		if original_set == 'pass':
			self.transform_fn = lambda x,y: (self.transformer(x), torch.tensor(1.,dtype=torch.float))
		else:
			if self.segmentation:
				self.transform_fn = lambda x,y: segmentation_transform(x,y, self.transformer, self.alb_normalizer, self.alb_reduction)
			else:
				if 'p' in dir(self.transformer):
					self.transform_fn = lambda x,y: (self.alb_normalizer(torch.tensor(np.transpose(self.transformer(True, image=x)['image']/self.alb_reduction,axes=(2,0,1)),dtype=torch.float)), torch.from_numpy(y).float())
				else:
					self.transform_fn = lambda x,y: (self.transformer(x), torch.from_numpy(y).float())
		self.image_set = dataset_df[dataset_df['set']==set]
		if self.segmentation:
			self.targets='target'
			self.image_set[self.targets] = 1.
		if not self.segmentation and fix_names:
			self.image_set.loc[:,'ID'] = self.image_set['ID'].map(fix_name)
		self.slide_info = pd.DataFrame(self.image_set.set_index('ID').loc[:,self.targets])
		if self.mt_bce and not self.segmentation:
			if pos_annotation_class:
				self.targets = [pos_annotation_class]+list(other_annotations)
			else:
				self.targets = None
		print(self.targets)
		IDs = self.slide_info.index.tolist()
		pi_dict=dict(input_info_db=patch_info_file,
					slide_labels=self.slide_info,
					pos_annotation_class=pos_annotation_class,
					patch_size=patch_size,
					segmentation=self.segmentation,
					other_annotations=other_annotations,
					target_segmentation_class=target_segmentation_class,
					target_threshold=target_threshold,
					classify_annotations=classify_annotations,
					modify_patches=modify_patches)
		self.patch_info = modify_patch_info(**pi_dict)

		if self.segmentation and original_set!='pass':
			#IDs = self.patch_info['ID'].unique()
			self.segmentation_maps = {slide:npy2da(join(input_dir,'{}_mask.npy'.format(slide))) for slide in IDs}
		self.slides = {slide:load_preprocessed_img(join(input_dir,'{}.zarr'.format(slide))) for slide in IDs}
		#print(self.slide_info)
		if original_set =='pass':
			self.segmentation=False
		#print(self.patch_info[self.targets].unique())
		if oversampling_factor > 1:
			self.patch_info = pd.concat([self.patch_info]*int(oversampling_factor),axis=0).reset_index(drop=True)
		elif oversampling_factor < 1:
			self.patch_info = self.patch_info.sample(frac=oversampling_factor).reset_index(drop=True)
		self.length = self.patch_info.shape[0]
		self.n_segmentation_classes = n_segmentation_classes
		self.gdl=gdl if self.segmentation else False
		self.binarized=False
		self.classify_annotations=classify_annotations
		print(self.targets)
		self.dilation_jitter=DilationJitter(dilation_jitter,self.segmentation,(original_set=='train'))
		if not self.targets:
			self.targets = [pos_annotation_class]+list(other_annotations)

	def concat(self, other_dataset):
		"""Concatenate this dataset with others. Updates its own internal attributes.

		Parameters
		----------
		other_dataset:DynamicImageDataset
			Other image dataset.

		"""
		self.patch_info = pd.concat([self.patch_info, other_dataset.patch_info],axis=0).reset_index(drop=True)
		self.length = self.patch_info.shape[0]
		if self.segmentation:
			self.segmentation_maps.update(other_dataset.segmentation_maps)
			#print(self.segmentation_maps.keys())

	def retain_ID(self, ID):
		"""Reduce the sample set to just images from one ID.

		Parameters
		----------
		ID:str
			Basename/ID to predict on.

		Returns
		-------
		self

		"""
		self.patch_info=self.patch_info.loc[self.patch_info['ID']==ID]
		self.length = self.patch_info.shape[0]
		self.segmentation_maps={ID:self.segmentation_maps[ID]}
		return self

	def split_by_ID(self):
		"""Generator similar to groupby, but splits up by ID, generates (ID,data) using retain_ID.

		Returns
		-------
		generator
			ID, DynamicDataset

		"""
		for ID in self.patch_info['ID'].unique():
			new_dataset = copy.deepcopy(self)
			yield ID, new_dataset.retain_ID(ID)

	def select_IDs(self, IDs):
		for ID in IDs:
			if ID in self.patch_info['ID'].unique():
				new_dataset = copy.deepcopy(self)
				yield ID, new_dataset.retain_ID(ID)


	def get_class_weights(self, i=0):#[0,1]
		"""Weight loss function with weights inversely proportional to the class appearence.

		Parameters
		----------
		i:int
			If multi-target, class used for weighting.

		Returns
		-------
		self
			Dataset.

		"""
		if self.segmentation:
			label_counts=self.patch_info[list(map(str,list(range(self.n_segmentation_classes))))].sum(axis=0).values
			freq = label_counts/sum(label_counts)
			weights=1./(freq)
		elif self.mt_bce:
			weights=1./(self.patch_info.loc[:,self.targets].sum(axis=0).values)
			weights=weights/sum(weights)
		else:
			if self.binarized and len(self.targets)>1:
				y=np.argmax(self.patch_info.loc[:,self.targets].values,axis=1)
			elif (type(self.targets)==type('')):
				y=self.patch_info.loc[:,self.targets]
			else:
				y=self.patch_info.loc[:,self.targets[i]]
			y=y.values.astype(int).flatten()
			weights=compute_class_weight(class_weight='balanced',classes=np.unique(y),y=y)
		return weights

	def binarize_annotations(self, binarizer=None, num_targets=1, binary_threshold=0.):
		"""Label binarize some annotations or threshold them if classifying slide annotations.

		Parameters
		----------
		binarizer:LabelBinarizer
			Binarizes the labels of a column(s)
		num_targets:int
			Number of desired targets to preidict on.
		binary_threshold:float
			Amount of annotation in patch before positive annotation.

		Returns
		-------
		binarizer

		"""

		annotations = self.patch_info['annotation']
		annots=[annot for annot in list(self.patch_info.iloc[:,6:]) if annot !='area']
		if not self.mt_bce and num_targets > 1:
			if binarizer == None:
				self.binarizer = LabelBinarizer().fit(annotations)
			else:
				self.binarizer = copy.deepcopy(binarizer)
			self.targets = self.binarizer.classes_
			annotation_labels = pd.DataFrame(self.binarizer.transform(annotations),index=self.patch_info.index,columns=self.targets).astype(float)
			for col in list(annotation_labels):
				if col in list(self.patch_info):
					self.patch_info.loc[:,col]=annotation_labels[col].values
				else:
					self.patch_info[col]=annotation_labels[col].values
		else:
			self.binarizer=None
			self.targets=annots
			if num_targets == 1:
				self.targets = [self.targets[-1]]
			if binary_threshold>0.:
				self.patch_info.loc[:,self.targets]=(self.patch_info[self.targets]>=binary_threshold).values.astype(np.float32)
			print(self.targets)
			#self.patch_info = pd.concat([self.patch_info,annotation_labels],axis=1)
		self.binarized=True
		return self.binarizer

	def subsample(self, p):
		"""Sample subset of dataset.

		Parameters
		----------
		p:float
			Fraction to subsample.

		"""
		np.random.seed(42)
		self.patch_info = self.patch_info.sample(frac=p)
		self.length = self.patch_info.shape[0]

	def update_dataset(self, input_dir, new_db, prediction_basename=[]):
		"""Experimental. Only use for segmentation for now."""
		self.input_dir=input_dir
		self.patch_info=load_sql_df(new_db, self.patch_size)
		IDs = self.patch_info['ID'].unique()
		self.slides = {slide:load_preprocessed_img(join(self.input_dir,'{}.zarr'.format(slide))) for slide in IDs}
		if self.segmentation:
			if prediction_basename:
				self.segmentation_maps = {slide:npy2da(join(self.input_dir,'{}_mask.npy'.format(slide))) for slide in IDs if slide in prediction_basename}
			else:
				self.segmentation_maps = {slide:npy2da(join(self.input_dir,'{}_mask.npy'.format(slide))) for slide in IDs}
		self.length = self.patch_info.shape[0]

	#@pysnooper.snoop("getitem.log")
	def __getitem__(self, i):
		patch_info = self.patch_info.iloc[i]
		ID = patch_info['ID']
		xs = patch_info['x']
		ys = patch_info['y']
		patch_size = patch_info['patch_size']
		if xs==np.nan:
			entire_image=True
		else:
			entire_image=False
		targets=self.targets
		use_long=False
		if not self.segmentation:
			y = patch_info.loc[list(self.targets) if not isinstance(self.targets,str) else self.targets]
			if isinstance(y,pd.Series):
				y=y.values.astype(float)
				if self.binarized and not self.mt_bce and len(y)>1:
					y=np.array(y.argmax())
					use_long=True
			y=np.array(y)
			if not y.shape:
				y=y.reshape(1)
		if self.segmentation:
			arr=self.segmentation_maps[ID]
			if not entire_image:
				arr=arr[xs:xs+patch_size,ys:ys+patch_size]
			arr=self.dilation_jitter(arr)
		y=(y if not self.segmentation else np.array(arr))
		#print(y)
		arr=self.slides[ID]
		if not entire_image:
			arr=arr[xs:xs+patch_size,ys:ys+patch_size,:3]
		image, y = self.transform_fn(arr.compute().astype(np.uint8), y)#.unsqueeze(0) # transpose .transpose([1,0,2])
		if not self.segmentation and not self.mt_bce and self.classify_annotations and use_long:
			y=y.long()
		#image_size=image.size()
		if self.gdl:
			y=class2one_hot(y, self.n_segmentation_classes)
		#	y=one_hot2dist(y)
		return image, y

	def __len__(self):
		return self.length

class NPYDataset(Dataset):
	def __init__(self, patch_info, patch_size, npy_file, transform, mmap=False):
		self.ID=os.path.basename(npy_file).split('.')[0]
		patch_info=patch_info=load_sql_df(patch_info,patch_size)
		self.patch_info=patch_info.loc[patch_info["ID"]==self.ID].reset_index()
		self.X=np.load(npy_file,mmap_mode=(None if not mmap else 'r+'))
		self.transform=transform

	def __getitem__(self,i):
		x,y,patch_size=self.patch_info.loc[i,["x","y","patch_size"]]
		return self.transform(self.X[x:x+patch_size,y:y+patch_size])

	def __len__(self):
		return self.patch_info.shape[0]

	def embed(self,model,batch_size,out_dir):
		Z=[]
		dataloader=DataLoader(self,batch_size=batch_size,shuffle=False)
		n_batches=len(self)//batch_size
		with torch.no_grad():
			for i,X in enumerate(dataloader):
				if torch.cuda.is_available():
					X=X.cuda()
				z=model(X).detach().cpu().numpy()
				Z.append(z)
				print(f"Processed batch {i}/{n_batches}")
		Z=np.vstack(Z)
		torch.save(dict(embeddings=Z,patch_info=self.patch_info),os.path.join(out_dir,f"{self.ID}.pkl"))
		print("Embeddings saved")
		quit()