diff --git a/datasets/lidc/configs.py b/datasets/lidc/configs.py
index 19b243f..0b828af 100644
--- a/datasets/lidc/configs.py
+++ b/datasets/lidc/configs.py
@@ -1,445 +1,445 @@
#!/usr/bin/env python
# Copyright 2019 Division of Medical Image Computing, German Cancer Research Center (DKFZ).
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
# ==============================================================================
import sys
import os
from collections import namedtuple
sys.path.append(os.path.dirname(os.path.realpath(__file__)))
import numpy as np
sys.path.append(os.path.dirname(os.path.realpath(__file__))+"/../..")
from default_configs import DefaultConfigs
# legends, nested classes are not handled well in multiprocessing! hence, Label class def in outer scope
Label = namedtuple("Label", ['id', 'name', 'color', 'm_scores']) # m_scores = malignancy scores
binLabel = namedtuple("binLabel", ['id', 'name', 'color', 'm_scores', 'bin_vals'])
class Configs(DefaultConfigs):
def __init__(self, server_env=None):
super(Configs, self).__init__(server_env)
#########################
# Preprocessing #
#########################
self.root_dir = '/home/gregor/networkdrives/E130-Personal/Goetz/Datenkollektive/Lungendaten/Nodules_LIDC_IDRI'
self.raw_data_dir = '{}/new_nrrd'.format(self.root_dir)
self.pp_dir = '/media/gregor/HDD2TB/data/lidc/pp_20200309_dev'
# 'merged' for one gt per image, 'single_annotator' for four gts per image.
self.gts_to_produce = ["single_annotator", "merged"]
self.target_spacing = (0.7, 0.7, 1.25)
#########################
# I/O #
#########################
# path to preprocessed data.
self.pp_name = 'pp_20190805'
self.input_df_name = 'info_df.pickle'
self.data_sourcedir = '/media/gregor/HDD2TB/data/lidc/{}/'.format(self.pp_name)
#self.data_sourcedir = '/home/gregor/networkdrives/E132-Cluster-Projects/lidc/data/{}/'.format(self.pp_name)
# settings for deployment on cluster.
if server_env:
# path to preprocessed data.
self.data_sourcedir = '/datasets/datasets_ramien/lidc/data/{}_npz/'.format(self.pp_name)
# one out of ['mrcnn', 'retina_net', 'retina_unet', 'detection_fpn'].
- self.model = 'mrcnn'
+ self.model = 'detection_fpn'
self.model_path = 'models/{}.py'.format(self.model if not 'retina' in self.model else 'retina_net')
self.model_path = os.path.join(self.source_dir, self.model_path)
#########################
# Architecture #
#########################
# dimension the model operates in. one out of [2, 3].
- self.dim = 2
+ self.dim = 3
# 'class': standard object classification per roi, pairwise combinable with each of below tasks.
# if 'class' is omitted from tasks, object classes will be fg/bg (1/0) from RPN.
# 'regression': regress some vector per each roi
# 'regression_ken_gal': use kendall-gal uncertainty sigma
# 'regression_bin': classify each roi into a bin related to a regression scale
self.prediction_tasks = ['class']
self.start_filts = 48 if self.dim == 2 else 18
self.end_filts = self.start_filts * 4 if self.dim == 2 else self.start_filts * 2
self.res_architecture = 'resnet50' # 'resnet101' , 'resnet50'
self.norm = "instance_norm" # one of None, 'instance_norm', 'batch_norm'
# one of 'xavier_uniform', 'xavier_normal', or 'kaiming_normal', None (=default = 'kaiming_uniform')
self.weight_init = None
self.regression_n_features = 1
#########################
# Data Loader #
#########################
# distorted gt experiments: train on single-annotator gts in a random fashion to investigate network's
# handling of noisy gts.
# choose 'merged' for single, merged gt per image, or 'single_annotator' for four gts per image.
# validation is always performed on same gt kind as training, testing always on merged gt.
self.training_gts = "merged"
# select modalities from preprocessed data
self.channels = [0]
self.n_channels = len(self.channels)
# patch_size to be used for training. pre_crop_size is the patch_size before data augmentation.
self.pre_crop_size_2D = [320, 320]
self.patch_size_2D = [320, 320]
self.pre_crop_size_3D = [160, 160, 96]
self.patch_size_3D = [160, 160, 96]
self.patch_size = self.patch_size_2D if self.dim == 2 else self.patch_size_3D
self.pre_crop_size = self.pre_crop_size_2D if self.dim == 2 else self.pre_crop_size_3D
# ratio of free sampled batch elements before class balancing is triggered
self.batch_random_ratio = 0.1
self.balance_target = "class_targets" if 'class' in self.prediction_tasks else 'rg_bin_targets'
# set 2D network to match 3D gt boxes.
self.merge_2D_to_3D_preds = self.dim==2
self.observables_rois = []
#self.rg_map = {1:1, 2:2, 3:3, 4:4, 5:5}
#########################
# Colors and Legends #
#########################
self.plot_frequency = 5
binary_cl_labels = [Label(1, 'benign', (*self.dark_green, 1.), (1, 2)),
Label(2, 'malignant', (*self.red, 1.), (3, 4, 5))]
quintuple_cl_labels = [Label(1, 'MS1', (*self.dark_green, 1.), (1,)),
Label(2, 'MS2', (*self.dark_yellow, 1.), (2,)),
Label(3, 'MS3', (*self.orange, 1.), (3,)),
Label(4, 'MS4', (*self.bright_red, 1.), (4,)),
Label(5, 'MS5', (*self.red, 1.), (5,))]
# choose here if to do 2-way or 5-way regression-bin classification
task_spec_cl_labels = quintuple_cl_labels
self.class_labels = [
# #id #name #color #malignancy score
Label( 0, 'bg', (*self.gray, 0.), (0,))]
if "class" in self.prediction_tasks:
self.class_labels += task_spec_cl_labels
else:
self.class_labels += [Label(1, 'lesion', (*self.orange, 1.), (1,2,3,4,5))]
if any(['regression' in task for task in self.prediction_tasks]):
self.bin_labels = [binLabel(0, 'MS0', (*self.gray, 1.), (0,), (0,))]
self.bin_labels += [binLabel(cll.id, cll.name, cll.color, cll.m_scores,
tuple([ms for ms in cll.m_scores])) for cll in task_spec_cl_labels]
self.bin_id2label = {label.id: label for label in self.bin_labels}
self.ms2bin_label = {ms: label for label in self.bin_labels for ms in label.m_scores}
bins = [(min(label.bin_vals), max(label.bin_vals)) for label in self.bin_labels]
self.bin_id2rg_val = {ix: [np.mean(bin)] for ix, bin in enumerate(bins)}
self.bin_edges = [(bins[i][1] + bins[i + 1][0]) / 2 for i in range(len(bins) - 1)]
if self.class_specific_seg:
self.seg_labels = self.class_labels
else:
self.seg_labels = [ # id #name #color
Label(0, 'bg', (*self.gray, 0.)),
Label(1, 'fg', (*self.orange, 1.))
]
self.class_id2label = {label.id: label for label in self.class_labels}
self.class_dict = {label.id: label.name for label in self.class_labels if label.id != 0}
# class_dict is used in evaluator / ap, auc, etc. statistics, and class 0 (bg) only needs to be
# evaluated in debugging
self.class_cmap = {label.id: label.color for label in self.class_labels}
self.seg_id2label = {label.id: label for label in self.seg_labels}
self.cmap = {label.id: label.color for label in self.seg_labels}
self.plot_prediction_histograms = True
self.plot_stat_curves = False
self.has_colorchannels = False
self.plot_class_ids = True
self.num_classes = len(self.class_dict) # for instance classification (excl background)
self.num_seg_classes = len(self.seg_labels) # incl background
#########################
# Data Augmentation #
#########################
self.da_kwargs={
'mirror': True,
'mirror_axes': tuple(np.arange(0, self.dim, 1)),
'do_elastic_deform': True,
'alpha':(0., 1500.),
'sigma':(30., 50.),
'do_rotation':True,
'angle_x': (0., 2 * np.pi),
'angle_y': (0., 0),
'angle_z': (0., 0),
'do_scale': True,
'scale':(0.8, 1.1),
'random_crop':False,
'rand_crop_dist': (self.patch_size[0] / 2. - 3, self.patch_size[1] / 2. - 3),
'border_mode_data': 'constant',
'border_cval_data': 0,
'order_data': 1}
if self.dim == 3:
self.da_kwargs['do_elastic_deform'] = False
self.da_kwargs['angle_x'] = (0, 0.0)
self.da_kwargs['angle_y'] = (0, 0.0) #must be 0!!
self.da_kwargs['angle_z'] = (0., 2 * np.pi)
#################################
# Schedule / Selection / Optim #
#################################
self.num_epochs = 130 if self.dim == 2 else 150
self.num_train_batches = 200 if self.dim == 2 else 200
self.batch_size = 20 if self.dim == 2 else 8
# decide whether to validate on entire patient volumes (like testing) or sampled patches (like training)
# the former is morge accurate, while the latter is faster (depending on volume size)
self.val_mode = 'val_sampling' # only 'val_sampling', 'val_patient' not implemented
if self.val_mode == 'val_patient':
raise NotImplementedError
if self.val_mode == 'val_sampling':
self.num_val_batches = 70
self.save_n_models = 4
# set a minimum epoch number for saving in case of instabilities in the first phase of training.
self.min_save_thresh = 0 if self.dim == 2 else 0
# criteria to average over for saving epochs, 'criterion':weight.
if "class" in self.prediction_tasks:
# 'criterion': weight
if len(self.class_labels)==3:
self.model_selection_criteria = {"benign_ap": 0.5, "malignant_ap": 0.5}
elif len(self.class_labels)==6:
self.model_selection_criteria = {str(label.name)+"_ap": 1./5 for label in self.class_labels if label.id!=0}
elif any("regression" in task for task in self.prediction_tasks):
self.model_selection_criteria = {"lesion_ap": 0.2, "lesion_avp": 0.8}
self.weight_decay = 1e-5
self.exclude_from_wd = []
self.clip_norm = 200 if 'regression_ken_gal' in self.prediction_tasks else None # number or None
# int in [0, dataset_size]. select n patients from dataset for prototyping. If None, all data is used.
self.select_prototype_subset = None #self.batch_size
#########################
# Testing #
#########################
# set the top-n-epochs to be saved for temporal averaging in testing.
self.test_n_epochs = self.save_n_models
self.test_aug_axes = (0,1,(0,1)) # None or list: choices are 0,1,(0,1) (0==spatial y, 1== spatial x).
self.held_out_test_set = False
self.max_test_patients = "all" # "all" or number
self.report_score_level = ['rois', 'patient'] # choose list from 'patient', 'rois'
self.patient_class_of_interest = 2 if 'class' in self.prediction_tasks else 1
self.metrics = ['ap', 'auc']
if any(['regression' in task for task in self.prediction_tasks]):
self.metrics += ['avp', 'rg_MAE_weighted', 'rg_MAE_weighted_tp',
'rg_bin_accuracy_weighted', 'rg_bin_accuracy_weighted_tp']
if 'aleatoric' in self.model:
self.metrics += ['rg_uncertainty', 'rg_uncertainty_tp', 'rg_uncertainty_tp_weighted']
self.evaluate_fold_means = True
self.ap_match_ious = [0.1] # list of ious to be evaluated for ap-scoring.
self.min_det_thresh = 0.1 # minimum confidence value to select predictions for evaluation.
# aggregation method for test and val_patient predictions.
# wbc = weighted box clustering as in https://arxiv.org/pdf/1811.08661.pdf,
# nms = standard non-maximum suppression, or None = no clustering
self.clustering = 'wbc'
# iou thresh (exclusive!) for regarding two preds as concerning the same ROI
self.clustering_iou = 0.1 # has to be larger than desired possible overlap iou of model predictions
self.plot_prediction_histograms = True
self.plot_stat_curves = False
self.n_test_plots = 1
#########################
# Assertions #
#########################
if not 'class' in self.prediction_tasks:
assert self.num_classes == 1
#########################
# Add model specifics #
#########################
{'detection_fpn': self.add_det_fpn_configs,
'mrcnn': self.add_mrcnn_configs, 'mrcnn_aleatoric': self.add_mrcnn_configs,
'retina_net': self.add_mrcnn_configs,
'retina_unet': self.add_mrcnn_configs,
}[self.model]()
def rg_val_to_bin_id(self, rg_val):
return float(np.digitize(np.mean(rg_val), self.bin_edges))
def add_det_fpn_configs(self):
self.learning_rate = [3e-4] * self.num_epochs
self.dynamic_lr_scheduling = False
# RoI score assigned to aggregation from pixel prediction (connected component). One of ['max', 'median'].
self.score_det = 'max'
# max number of roi candidates to identify per batch element and class.
self.n_roi_candidates = 10 if self.dim == 2 else 30
# loss mode: either weighted cross entropy ('wce'), batch-wise dice loss ('dice), or the sum of both ('dice_wce')
self.seg_loss_mode = 'wce'
# if <1, false positive predictions in foreground are penalized less.
self.fp_dice_weight = 1 if self.dim == 2 else 1
if len(self.class_labels)==3:
self.wce_weights = [1., 1., 1.] if self.seg_loss_mode=="dice_wce" else [0.1, 1., 1.]
elif len(self.class_labels)==6:
self.wce_weights = [1., 1., 1., 1., 1., 1.] if self.seg_loss_mode == "dice_wce" else [0.1, 1., 1., 1., 1., 1.]
else:
raise Exception("mismatch loss weights & nr of classes")
self.detection_min_confidence = self.min_det_thresh
self.head_classes = self.num_seg_classes
def add_mrcnn_configs(self):
# learning rate is a list with one entry per epoch.
self.learning_rate = [3e-4] * self.num_epochs
self.dynamic_lr_scheduling = False
# disable the re-sampling of mask proposals to original size for speed-up.
# since evaluation is detection-driven (box-matching) and not instance segmentation-driven (iou-matching),
# mask-outputs are optional.
self.return_masks_in_train = False
self.return_masks_in_val = True
self.return_masks_in_test = False
# set number of proposal boxes to plot after each epoch.
self.n_plot_rpn_props = 5 if self.dim == 2 else 30
# number of classes for network heads: n_foreground_classes + 1 (background)
self.head_classes = self.num_classes + 1
self.frcnn_mode = False
# feature map strides per pyramid level are inferred from architecture.
self.backbone_strides = {'xy': [4, 8, 16, 32], 'z': [1, 2, 4, 8]}
# anchor scales are chosen according to expected object sizes in data set. Default uses only one anchor scale
# per pyramid level. (outer list are pyramid levels (corresponding to BACKBONE_STRIDES), inner list are scales per level.)
self.rpn_anchor_scales = {'xy': [[8], [16], [32], [64]], 'z': [[2], [4], [8], [16]]}
# choose which pyramid levels to extract features from: P2: 0, P3: 1, P4: 2, P5: 3.
self.pyramid_levels = [0, 1, 2, 3]
# number of feature maps in rpn. typically lowered in 3D to save gpu-memory.
self.n_rpn_features = 512 if self.dim == 2 else 128
# anchor ratios and strides per position in feature maps.
self.rpn_anchor_ratios = [0.5, 1, 2]
self.rpn_anchor_stride = 1
# Threshold for first stage (RPN) non-maximum suppression (NMS): LOWER == HARDER SELECTION
self.rpn_nms_threshold = 0.7 if self.dim == 2 else 0.7
# loss sampling settings.
self.rpn_train_anchors_per_image = 64 #per batch element
self.train_rois_per_image = 6 #per batch element
self.roi_positive_ratio = 0.5
self.anchor_matching_iou = 0.7
# factor of top-k candidates to draw from per negative sample (stochastic-hard-example-mining).
# poolsize to draw top-k candidates from will be shem_poolsize * n_negative_samples.
self.shem_poolsize = 10
self.pool_size = (7, 7) if self.dim == 2 else (7, 7, 3)
self.mask_pool_size = (14, 14) if self.dim == 2 else (14, 14, 5)
self.mask_shape = (28, 28) if self.dim == 2 else (28, 28, 10)
self.rpn_bbox_std_dev = np.array([0.1, 0.1, 0.1, 0.2, 0.2, 0.2])
self.bbox_std_dev = np.array([0.1, 0.1, 0.1, 0.2, 0.2, 0.2])
self.window = np.array([0, 0, self.patch_size[0], self.patch_size[1], 0, self.patch_size_3D[2]])
self.scale = np.array([self.patch_size[0], self.patch_size[1], self.patch_size[0], self.patch_size[1],
self.patch_size_3D[2], self.patch_size_3D[2]])
if self.dim == 2:
self.rpn_bbox_std_dev = self.rpn_bbox_std_dev[:4]
self.bbox_std_dev = self.bbox_std_dev[:4]
self.window = self.window[:4]
self.scale = self.scale[:4]
# pre-selection in proposal-layer (stage 1) for NMS-speedup. applied per batch element.
self.pre_nms_limit = 3000 if self.dim == 2 else 6000
# n_proposals to be selected after NMS per batch element. too high numbers blow up memory if "detect_while_training" is True,
# since proposals of the entire batch are forwarded through second stage in as one "batch".
self.roi_chunk_size = 2500 if self.dim == 2 else 600
self.post_nms_rois_training = 500 if self.dim == 2 else 75
self.post_nms_rois_inference = 500
# Final selection of detections (refine_detections)
self.model_max_instances_per_batch_element = 10 if self.dim == 2 else 30 # per batch element and class.
self.detection_nms_threshold = 1e-5 # needs to be > 0, otherwise all predictions are one cluster.
self.model_min_confidence = 0.1
if self.dim == 2:
self.backbone_shapes = np.array(
[[int(np.ceil(self.patch_size[0] / stride)),
int(np.ceil(self.patch_size[1] / stride))]
for stride in self.backbone_strides['xy']])
else:
self.backbone_shapes = np.array(
[[int(np.ceil(self.patch_size[0] / stride)),
int(np.ceil(self.patch_size[1] / stride)),
int(np.ceil(self.patch_size[2] / stride_z))]
for stride, stride_z in zip(self.backbone_strides['xy'], self.backbone_strides['z']
)])
if self.model == 'retina_net' or self.model == 'retina_unet':
self.focal_loss = True
# implement extra anchor-scales according to retina-net publication.
self.rpn_anchor_scales['xy'] = [[ii[0], ii[0] * (2 ** (1 / 3)), ii[0] * (2 ** (2 / 3))] for ii in
self.rpn_anchor_scales['xy']]
self.rpn_anchor_scales['z'] = [[ii[0], ii[0] * (2 ** (1 / 3)), ii[0] * (2 ** (2 / 3))] for ii in
self.rpn_anchor_scales['z']]
self.n_anchors_per_pos = len(self.rpn_anchor_ratios) * 3
self.n_rpn_features = 256 if self.dim == 2 else 128
# pre-selection of detections for NMS-speedup. per entire batch.
self.pre_nms_limit = (500 if self.dim == 2 else 6250) * self.batch_size
# anchor matching iou is lower than in Mask R-CNN according to https://arxiv.org/abs/1708.02002
self.anchor_matching_iou = 0.5
if self.model == 'retina_unet':
self.operate_stride1 = True
diff --git a/datasets/toy/configs.py b/datasets/toy/configs.py
index 10d5889..6780d22 100644
--- a/datasets/toy/configs.py
+++ b/datasets/toy/configs.py
@@ -1,491 +1,491 @@
#!/usr/bin/env python
# Copyright 2019 Division of Medical Image Computing, German Cancer Research Center (DKFZ).
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
# ==============================================================================
import sys
import os
sys.path.append(os.path.dirname(os.path.realpath(__file__)))
import numpy as np
from default_configs import DefaultConfigs
from collections import namedtuple
boxLabel = namedtuple('boxLabel', ["name", "color"])
Label = namedtuple("Label", ['id', 'name', 'shape', 'radius', 'color', 'regression', 'ambiguities', 'gt_distortion'])
binLabel = namedtuple("binLabel", ['id', 'name', 'color', 'bin_vals'])
class Configs(DefaultConfigs):
def __init__(self, server_env=None):
super(Configs, self).__init__(server_env)
#########################
# Prepro #
#########################
self.pp_rootdir = os.path.join('/home/gregor/datasets/toy', "cyl1ps_dev")
self.pp_npz_dir = self.pp_rootdir+"_npz"
self.pre_crop_size = [320,320,8] #y,x,z; determines pp data shape (2D easily implementable, but only 3D for now)
self.min_2d_radius = 6 #in pixels
self.n_train_samples, self.n_test_samples = 1200, 1000
# not actually real one-hot encoding (ohe) but contains more info: roi-overlap only within classes.
self.pp_create_ohe_seg = False
self.pp_empty_samples_ratio = 0.1
self.pp_place_radii_mid_bin = True
self.pp_only_distort_2d = True
# outer-most intensity of blurred radii, relative to inner-object intensity. <1 for decreasing, > 1 for increasing.
# e.g.: setting 0.1 means blurred edge has min intensity 10% as large as inner-object intensity.
self.pp_blur_min_intensity = 0.2
self.max_instances_per_sample = 1 #how many max instances over all classes per sample (img if 2d, vol if 3d)
self.max_instances_per_class = self.max_instances_per_sample # how many max instances per image per class
self.noise_scale = 0. # std-dev of gaussian noise
self.ambigs_sampling = "gaussian" #"gaussian" or "uniform"
""" radius_calib: gt distort for calibrating uncertainty. Range of gt distortion is inferable from
image by distinguishing it from the rest of the object.
blurring width around edge will be shifted so that symmetric rel to orig radius.
blurring scale: if self.ambigs_sampling is uniform, distribution's non-zero range (b-a) will be sqrt(12)*scale
since uniform dist has variance (b-a)²/12. b,a will be placed symmetrically around unperturbed radius.
if sampling is gaussian, then scale parameter sets one std dev, i.e., blurring width will be orig_radius * std_dev * 2.
"""
self.ambiguities = {
#set which classes to apply which ambs to below in class labels
#choose out of: 'outer_radius', 'inner_radius', 'radii_relations'.
#kind #probability #scale (gaussian std, relative to unperturbed value)
#"outer_radius": (1., 0.5),
#"outer_radius_xy": (1., 0.5),
#"inner_radius": (0.5, 0.1),
#"radii_relations": (0.5, 0.1),
"radius_calib": (1., 1./6)
}
# shape choices: 'cylinder', 'block'
# id, name, shape, radius, color, regression, ambiguities, gt_distortion
self.pp_classes = [Label(1, 'cylinder', 'cylinder', ((6,6,1),(40,40,8)), (*self.blue, 1.), "radius_2d", (), ()),
#Label(2, 'block', 'block', ((6,6,1),(40,40,8)), (*self.aubergine,1.), "radii_2d", (), ('radius_calib',))
]
#########################
# I/O #
#########################
self.data_sourcedir = '/home/gregor/datasets/toy/cyl1ps_dev'
if server_env:
self.data_sourcedir = '/datasets/datasets_ramien/toy/data/cyl1ps_dev_npz'
self.test_data_sourcedir = os.path.join(self.data_sourcedir, 'test')
self.data_sourcedir = os.path.join(self.data_sourcedir, "train")
self.info_df_name = 'info_df.pickle'
# one out of ['mrcnn', 'retina_net', 'retina_unet', 'detection_unet', 'ufrcnn', 'detection_fpn'].
self.model = 'mrcnn'
self.model_path = 'models/{}.py'.format(self.model if not 'retina' in self.model else 'retina_net')
self.model_path = os.path.join(self.source_dir, self.model_path)
#########################
# Architecture #
#########################
# one out of [2, 3]. dimension the model operates in.
self.dim = 2
# 'class', 'regression', 'regression_bin', 'regression_ken_gal'
# currently only tested mode is a single-task at a time (i.e., only one task in below list)
# but, in principle, tasks could be combined (e.g., object classes and regression per class)
self.prediction_tasks = ['class', ]
self.start_filts = 48 if self.dim == 2 else 18
self.end_filts = self.start_filts * 4 if self.dim == 2 else self.start_filts * 2
self.res_architecture = 'resnet50' # 'resnet101' , 'resnet50'
self.norm = 'instance_norm' # one of None, 'instance_norm', 'batch_norm'
self.relu = 'relu'
# one of 'xavier_uniform', 'xavier_normal', or 'kaiming_normal', None (=default = 'kaiming_uniform')
self.weight_init = None
self.regression_n_features = 1 # length of regressor target vector
#########################
# Data Loader #
#########################
self.num_epochs = 24
self.num_train_batches = 100 if self.dim == 2 else 180
self.batch_size = 20 if self.dim == 2 else 8
self.n_cv_splits = 4
# select modalities from preprocessed data
self.channels = [0]
self.n_channels = len(self.channels)
# which channel (mod) to show as bg in plotting, will be extra added to batch if not in self.channels
self.plot_bg_chan = 0
self.crop_margin = [20, 20, 1] # has to be smaller than respective patch_size//2
self.patch_size_2D = self.pre_crop_size[:2]
self.patch_size_3D = self.pre_crop_size[:2]+[8]
# patch_size to be used for training. pre_crop_size is the patch_size before data augmentation.
self.patch_size = self.patch_size_2D if self.dim == 2 else self.patch_size_3D
# ratio of free sampled batch elements before class balancing is triggered
# (>0 to include "empty"/background patches.)
self.batch_random_ratio = 0.2
self.balance_target = "class_targets" if 'class' in self.prediction_tasks else "rg_bin_targets"
self.observables_patient = []
self.observables_rois = []
self.seed = 3 #for generating folds
#############################
# Colors, Classes, Legends #
#############################
self.plot_frequency = 4
binary_bin_labels = [binLabel(1, 'r<=25', (*self.green, 1.), (1,25)),
binLabel(2, 'r>25', (*self.red, 1.), (25,))]
quintuple_bin_labels = [binLabel(1, 'r2-10', (*self.green, 1.), (2,10)),
binLabel(2, 'r10-20', (*self.yellow, 1.), (10,20)),
binLabel(3, 'r20-30', (*self.orange, 1.), (20,30)),
binLabel(4, 'r30-40', (*self.bright_red, 1.), (30,40)),
binLabel(5, 'r>40', (*self.red, 1.), (40,))]
# choose here if to do 2-way or 5-way regression-bin classification
task_spec_bin_labels = quintuple_bin_labels
self.class_labels = [
# regression: regression-task label, either value or "(x,y,z)_radius" or "radii".
# ambiguities: name of above defined ambig to apply to image data (not gt); need to be iterables!
# gt_distortion: name of ambig to apply to gt only; needs to be iterable!
# #id #name #shape #radius #color #regression #ambiguities #gt_distortion
Label( 0, 'bg', None, (0, 0, 0), (*self.white, 0.), (0, 0, 0), (), ())]
if "class" in self.prediction_tasks:
self.class_labels += self.pp_classes
else:
self.class_labels += [Label(1, 'object', 'object', ('various',), (*self.orange, 1.), ('radius_2d',), ("various",), ('various',))]
if any(['regression' in task for task in self.prediction_tasks]):
self.bin_labels = [binLabel(0, 'bg', (*self.white, 1.), (0,))]
self.bin_labels += task_spec_bin_labels
self.bin_id2label = {label.id: label for label in self.bin_labels}
bins = [(min(label.bin_vals), max(label.bin_vals)) for label in self.bin_labels]
self.bin_id2rg_val = {ix: [np.mean(bin)] for ix, bin in enumerate(bins)}
self.bin_edges = [(bins[i][1] + bins[i + 1][0]) / 2 for i in range(len(bins) - 1)]
self.bin_dict = {label.id: label.name for label in self.bin_labels if label.id != 0}
if self.class_specific_seg:
self.seg_labels = self.class_labels
self.box_type2label = {label.name: label for label in self.box_labels}
self.class_id2label = {label.id: label for label in self.class_labels}
self.class_dict = {label.id: label.name for label in self.class_labels if label.id != 0}
self.seg_id2label = {label.id: label for label in self.seg_labels}
self.cmap = {label.id: label.color for label in self.seg_labels}
self.plot_prediction_histograms = True
self.plot_stat_curves = False
self.has_colorchannels = False
self.plot_class_ids = True
self.num_classes = len(self.class_dict)
self.num_seg_classes = len(self.seg_labels)
#########################
# Data Augmentation #
#########################
self.do_aug = True
self.da_kwargs = {
'mirror': True,
'mirror_axes': tuple(np.arange(0, self.dim, 1)),
'do_elastic_deform': False,
'alpha': (500., 1500.),
'sigma': (40., 45.),
'do_rotation': False,
'angle_x': (0., 2 * np.pi),
'angle_y': (0., 0),
'angle_z': (0., 0),
'do_scale': False,
'scale': (0.8, 1.1),
'random_crop': False,
'rand_crop_dist': (self.patch_size[0] / 2. - 3, self.patch_size[1] / 2. - 3),
'border_mode_data': 'constant',
'border_cval_data': 0,
'order_data': 1
}
if self.dim == 3:
self.da_kwargs['do_elastic_deform'] = False
self.da_kwargs['angle_x'] = (0, 0.0)
self.da_kwargs['angle_y'] = (0, 0.0) # must be 0!!
self.da_kwargs['angle_z'] = (0., 2 * np.pi)
#########################
# Schedule / Selection #
#########################
# decide whether to validate on entire patient volumes (like testing) or sampled patches (like training)
# the former is morge accurate, while the latter is faster (depending on volume size)
self.val_mode = 'val_sampling' # one of 'val_sampling' , 'val_patient'
if self.val_mode == 'val_patient':
self.max_val_patients = 220 # if 'all' iterates over entire val_set once.
if self.val_mode == 'val_sampling':
self.num_val_batches = 35 if self.dim==2 else 25
self.save_n_models = 2
self.min_save_thresh = 1 if self.dim == 2 else 1 # =wait time in epochs
if "class" in self.prediction_tasks:
self.model_selection_criteria = {name + "_ap": 1. for name in self.class_dict.values()}
elif any("regression" in task for task in self.prediction_tasks):
self.model_selection_criteria = {name + "_ap": 0.2 for name in self.class_dict.values()}
self.model_selection_criteria.update({name + "_avp": 0.8 for name in self.class_dict.values()})
self.lr_decay_factor = 0.25
- self.scheduling_patience = np.ceil(1800 / (self.num_train_batches * self.batch_size))
+ self.scheduling_patience = np.ceil(3600 / (self.num_train_batches * self.batch_size))
self.weight_decay = 3e-5
self.exclude_from_wd = []
self.clip_norm = None # number or None
#########################
# Testing / Plotting #
#########################
self.test_aug_axes = (0,1,(0,1)) # None or list: choices are 0,1,(0,1)
self.held_out_test_set = True
self.max_test_patients = "all" # number or "all" for all
self.test_against_exact_gt = True # only True implemented
self.val_against_exact_gt = False # True is an unrealistic --> irrelevant scenario.
self.report_score_level = ['rois'] # 'patient' or 'rois' (incl)
self.patient_class_of_interest = 1
self.patient_bin_of_interest = 2
self.eval_bins_separately = False#"additionally" if not 'class' in self.prediction_tasks else False
self.metrics = ['ap', 'auc', 'dice']
if any(['regression' in task for task in self.prediction_tasks]):
self.metrics += ['avp', 'rg_MAE_weighted', 'rg_MAE_weighted_tp',
'rg_bin_accuracy_weighted', 'rg_bin_accuracy_weighted_tp']
if 'aleatoric' in self.model:
self.metrics += ['rg_uncertainty', 'rg_uncertainty_tp', 'rg_uncertainty_tp_weighted']
self.evaluate_fold_means = True
self.ap_match_ious = [0.5] # threshold(s) for considering a prediction as true positive
self.min_det_thresh = 0.3
self.model_max_iou_resolution = 0.2
# aggregation method for test and val_patient predictions.
# wbc = weighted box clustering as in https://arxiv.org/pdf/1811.08661.pdf,
# nms = standard non-maximum suppression, or None = no clustering
self.clustering = 'wbc'
# iou thresh (exclusive!) for regarding two preds as concerning the same ROI
self.clustering_iou = self.model_max_iou_resolution # has to be larger than desired possible overlap iou of model predictions
self.merge_2D_to_3D_preds = self.dim==2
self.merge_3D_iou = self.model_max_iou_resolution
self.n_test_plots = 1 # per fold and rank
self.test_n_epochs = self.save_n_models # should be called n_test_ens, since is number of models to ensemble over during testing
# is multiplied by (1 + nr of test augs)
#########################
# Assertions #
#########################
if not 'class' in self.prediction_tasks:
assert self.num_classes == 1
#########################
# Add model specifics #
#########################
{'mrcnn': self.add_mrcnn_configs, 'mrcnn_aleatoric': self.add_mrcnn_configs,
'retina_net': self.add_mrcnn_configs, 'retina_unet': self.add_mrcnn_configs,
'detection_unet': self.add_det_unet_configs, 'detection_fpn': self.add_det_fpn_configs
}[self.model]()
def rg_val_to_bin_id(self, rg_val):
#only meant for isotropic radii!!
# only 2D radii (x and y dims) or 1D (x or y) are expected
return np.round(np.digitize(rg_val, self.bin_edges).mean())
def add_det_fpn_configs(self):
- self.learning_rate = [3 * 1e-4] * self.num_epochs
+ self.learning_rate = [1 * 1e-4] * self.num_epochs
self.dynamic_lr_scheduling = True
self.scheduling_criterion = 'torch_loss'
self.scheduling_mode = 'min' if "loss" in self.scheduling_criterion else 'max'
self.n_roi_candidates = 4 if self.dim == 2 else 6
# max number of roi candidates to identify per image (slice in 2D, volume in 3D)
# loss mode: either weighted cross entropy ('wce'), batch-wise dice loss ('dice), or the sum of both ('dice_wce')
self.seg_loss_mode = 'wce'
self.wce_weights = [1] * self.num_seg_classes if 'dice' in self.seg_loss_mode else [0.1, 1]
self.fp_dice_weight = 1 if self.dim == 2 else 1
# if <1, false positive predictions in foreground are penalized less.
self.detection_min_confidence = 0.05
# how to determine score of roi: 'max' or 'median'
self.score_det = 'max'
def add_det_unet_configs(self):
self.learning_rate = [3 * 1e-4] * self.num_epochs
self.dynamic_lr_scheduling = True
self.scheduling_criterion = "torch_loss"
self.scheduling_mode = 'min' if "loss" in self.scheduling_criterion else 'max'
# max number of roi candidates to identify per image (slice in 2D, volume in 3D)
self.n_roi_candidates = 4 if self.dim == 2 else 6
# loss mode: either weighted cross entropy ('wce'), batch-wise dice loss ('dice), or the sum of both ('dice_wce')
self.seg_loss_mode = 'wce'
self.wce_weights = [1] * self.num_seg_classes if 'dice' in self.seg_loss_mode else [0.1, 1]
# if <1, false positive predictions in foreground are penalized less.
self.fp_dice_weight = 1 if self.dim == 2 else 1
self.detection_min_confidence = 0.05
# how to determine score of roi: 'max' or 'median'
self.score_det = 'max'
self.init_filts = 32
self.kernel_size = 3 # ks for horizontal, normal convs
self.kernel_size_m = 2 # ks for max pool
self.pad = "same" # "same" or integer, padding of horizontal convs
def add_mrcnn_configs(self):
self.learning_rate = [3e-4] * self.num_epochs
self.dynamic_lr_scheduling = True # with scheduler set in exec
self.scheduling_criterion = max(self.model_selection_criteria, key=self.model_selection_criteria.get)
self.scheduling_mode = 'min' if "loss" in self.scheduling_criterion else 'max'
# number of classes for network heads: n_foreground_classes + 1 (background)
self.head_classes = self.num_classes + 1 if 'class' in self.prediction_tasks else 2
# feed +/- n neighbouring slices into channel dimension. set to None for no context.
self.n_3D_context = None
if self.n_3D_context is not None and self.dim == 2:
self.n_channels *= (self.n_3D_context * 2 + 1)
self.detect_while_training = True
# disable the re-sampling of mask proposals to original size for speed-up.
# since evaluation is detection-driven (box-matching) and not instance segmentation-driven (iou-matching),
# mask outputs are optional.
self.return_masks_in_train = True
self.return_masks_in_val = True
self.return_masks_in_test = True
# feature map strides per pyramid level are inferred from architecture. anchor scales are set accordingly.
self.backbone_strides = {'xy': [4, 8, 16, 32], 'z': [1, 2, 4, 8]}
# anchor scales are chosen according to expected object sizes in data set. Default uses only one anchor scale
# per pyramid level. (outer list are pyramid levels (corresponding to BACKBONE_STRIDES), inner list are scales per level.)
self.rpn_anchor_scales = {'xy': [[4], [8], [16], [32]], 'z': [[1], [2], [4], [8]]}
# choose which pyramid levels to extract features from: P2: 0, P3: 1, P4: 2, P5: 3.
self.pyramid_levels = [0, 1, 2, 3]
# number of feature maps in rpn. typically lowered in 3D to save gpu-memory.
self.n_rpn_features = 512 if self.dim == 2 else 64
# anchor ratios and strides per position in feature maps.
self.rpn_anchor_ratios = [0.5, 1., 2.]
self.rpn_anchor_stride = 1
# Threshold for first stage (RPN) non-maximum suppression (NMS): LOWER == HARDER SELECTION
self.rpn_nms_threshold = max(0.7, self.model_max_iou_resolution)
# loss sampling settings.
self.rpn_train_anchors_per_image = 32
self.train_rois_per_image = 6 # per batch_instance
self.roi_positive_ratio = 0.5
self.anchor_matching_iou = 0.8
# k negative example candidates are drawn from a pool of size k*shem_poolsize (stochastic hard-example mining),
# where k<=#positive examples.
self.shem_poolsize = 6
self.pool_size = (7, 7) if self.dim == 2 else (7, 7, 3)
self.mask_pool_size = (14, 14) if self.dim == 2 else (14, 14, 5)
self.mask_shape = (28, 28) if self.dim == 2 else (28, 28, 10)
self.rpn_bbox_std_dev = np.array([0.1, 0.1, 0.1, 0.2, 0.2, 0.2])
self.bbox_std_dev = np.array([0.1, 0.1, 0.1, 0.2, 0.2, 0.2])
self.window = np.array([0, 0, self.patch_size[0], self.patch_size[1], 0, self.patch_size_3D[2]])
self.scale = np.array([self.patch_size[0], self.patch_size[1], self.patch_size[0], self.patch_size[1],
self.patch_size_3D[2], self.patch_size_3D[2]]) # y1,x1,y2,x2,z1,z2
if self.dim == 2:
self.rpn_bbox_std_dev = self.rpn_bbox_std_dev[:4]
self.bbox_std_dev = self.bbox_std_dev[:4]
self.window = self.window[:4]
self.scale = self.scale[:4]
self.plot_y_max = 1.5
self.n_plot_rpn_props = 5 if self.dim == 2 else 30 # per batch_instance (slice in 2D / patient in 3D)
# pre-selection in proposal-layer (stage 1) for NMS-speedup. applied per batch element.
self.pre_nms_limit = 2000 if self.dim == 2 else 4000
# n_proposals to be selected after NMS per batch element. too high numbers blow up memory if "detect_while_training" is True,
# since proposals of the entire batch are forwarded through second stage as one "batch".
self.roi_chunk_size = 1300 if self.dim == 2 else 500
self.post_nms_rois_training = 200 * (self.head_classes-1) if self.dim == 2 else 400
self.post_nms_rois_inference = 200 * (self.head_classes-1)
# Final selection of detections (refine_detections)
self.model_max_instances_per_batch_element = 9 if self.dim == 2 else 18 # per batch element and class.
self.detection_nms_threshold = self.model_max_iou_resolution # needs to be > 0, otherwise all predictions are one cluster.
self.model_min_confidence = 0.2 # iou for nms in box refining (directly after heads), should be >0 since ths>=x in mrcnn.py
if self.dim == 2:
self.backbone_shapes = np.array(
[[int(np.ceil(self.patch_size[0] / stride)),
int(np.ceil(self.patch_size[1] / stride))]
for stride in self.backbone_strides['xy']])
else:
self.backbone_shapes = np.array(
[[int(np.ceil(self.patch_size[0] / stride)),
int(np.ceil(self.patch_size[1] / stride)),
int(np.ceil(self.patch_size[2] / stride_z))]
for stride, stride_z in zip(self.backbone_strides['xy'], self.backbone_strides['z']
)])
if self.model == 'retina_net' or self.model == 'retina_unet':
# whether to use focal loss or SHEM for loss-sample selection
self.focal_loss = False
# implement extra anchor-scales according to https://arxiv.org/abs/1708.02002
self.rpn_anchor_scales['xy'] = [[ii[0], ii[0] * (2 ** (1 / 3)), ii[0] * (2 ** (2 / 3))] for ii in
self.rpn_anchor_scales['xy']]
self.rpn_anchor_scales['z'] = [[ii[0], ii[0] * (2 ** (1 / 3)), ii[0] * (2 ** (2 / 3))] for ii in
self.rpn_anchor_scales['z']]
self.n_anchors_per_pos = len(self.rpn_anchor_ratios) * 3
# pre-selection of detections for NMS-speedup. per entire batch.
self.pre_nms_limit = (500 if self.dim == 2 else 6250) * self.batch_size
# anchor matching iou is lower than in Mask R-CNN according to https://arxiv.org/abs/1708.02002
self.anchor_matching_iou = 0.7
if self.model == 'retina_unet':
self.operate_stride1 = True
diff --git a/datasets/toy_mdt/configs.py b/datasets/toy_mdt/configs.py
index 6533192..2333921 100644
--- a/datasets/toy_mdt/configs.py
+++ b/datasets/toy_mdt/configs.py
@@ -1,355 +1,355 @@
#!/usr/bin/env python
# Copyright 2018 Division of Medical Image Computing, German Cancer Research Center (DKFZ).
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
# ==============================================================================
import sys
import os
sys.path.append(os.path.dirname(os.path.realpath(__file__)))
import numpy as np
from collections import namedtuple
from default_configs import DefaultConfigs
Label = namedtuple("Label", ['id', 'name', 'color'])
class Configs(DefaultConfigs):
def __init__(self, server_env=None):
#########################
# Preprocessing #
#########################
self.root_dir = '/home/gregor/datasets/toy_mdt'
#########################
# I/O #
#########################
# one out of [2, 3]. dimension the model operates in.
self.dim = 2
DefaultConfigs.__init__(self, server_env, self.dim)
# one out of ['mrcnn', 'retina_net', 'retina_unet', 'detection_unet', 'ufrcnn'].
- self.model = 'mrcnn'
+ self.model = 'retina_unet'
self.model_path = 'models/{}.py'.format(self.model if not 'retina' in self.model else 'retina_net')
self.model_path = os.path.join(self.source_dir, self.model_path)
# int [0 < dataset_size]. select n patients from dataset for prototyping.
self.select_prototype_subset = None
self.held_out_test_set = True
# including val set. will be 3/4 train, 1/4 val.
self.n_train_val_data = 2500
# choose one of the 3 toy experiments described in https://arxiv.org/pdf/1811.08661.pdf
# one of ['donuts_shape', 'donuts_pattern', 'circles_scale'].
toy_mode = 'donuts_shape_noise'
# path to preprocessed data.
self.info_df_name = 'info_df.pickle'
self.pp_name = os.path.join(toy_mode, 'train')
self.data_sourcedir = os.path.join(self.root_dir, self.pp_name)
self.pp_test_name = os.path.join(toy_mode, 'test')
self.test_data_sourcedir = os.path.join(self.root_dir, self.pp_test_name)
# settings for deployment in cloud.
if server_env:
# path to preprocessed data.
pp_root_dir = '/datasets/datasets_ramien/toy_exp/data'
self.pp_name = os.path.join(toy_mode, 'train')
self.data_sourcedir = os.path.join(pp_root_dir, self.pp_name)
self.pp_test_name = os.path.join(toy_mode, 'test')
self.test_data_sourcedir = os.path.join(pp_root_dir, self.pp_test_name)
self.select_prototype_subset = None
#########################
# Data Loader #
#########################
# select modalities from preprocessed data
self.channels = [0]
self.n_channels = len(self.channels)
self.plot_bg_chan = 0
# patch_size to be used for training. pre_crop_size is the patch_size before data augmentation.
self.pre_crop_size_2D = [320, 320]
self.patch_size_2D = [320, 320]
self.patch_size = self.patch_size_2D if self.dim == 2 else self.patch_size_3D
self.pre_crop_size = self.pre_crop_size_2D if self.dim == 2 else self.pre_crop_size_3D
# ratio of free sampled batch elements before class balancing is triggered
# (>0 to include "empty"/background patches.)
self.batch_random_ratio = 0.2
# set 2D network to operate in 3D images.
self.merge_2D_to_3D_preds = False
#########################
# Architecture #
#########################
self.start_filts = 48 if self.dim == 2 else 18
self.end_filts = self.start_filts * 4 if self.dim == 2 else self.start_filts * 2
self.res_architecture = 'resnet50' # 'resnet101' , 'resnet50'
self.norm = "instance_norm" # one of None, 'instance_norm', 'batch_norm'
# one of 'xavier_uniform', 'xavier_normal', or 'kaiming_normal', None (=default = 'kaiming_uniform')
self.weight_init = "xavier_uniform"
# compatibility
self.regression_n_features = 1
self.num_classes = 2 # excluding bg
self.num_seg_classes = 3 # incl bg
#########################
# Schedule / Selection #
#########################
- self.num_epochs = 24
+ self.num_epochs = 26
self.num_train_batches = 100 if self.dim == 2 else 200
self.batch_size = 20 if self.dim == 2 else 8
self.do_validation = True
# decide whether to validate on entire patient volumes (like testing) or sampled patches (like training)
# the former is morge accurate, while the latter is faster (depending on volume size)
self.val_mode = 'val_patient' # one of 'val_sampling' , 'val_patient'
if self.val_mode == 'val_patient':
self.max_val_patients = "all" # if 'None' iterates over entire val_set once.
if self.val_mode == 'val_sampling':
self.num_val_batches = 50
self.optimizer = "ADAMW"
# set dynamic_lr_scheduling to True to apply LR scheduling with below settings.
self.dynamic_lr_scheduling = True
self.lr_decay_factor = 0.25
- self.scheduling_patience = np.ceil(2400 / (self.num_train_batches * self.batch_size))
+ self.scheduling_patience = np.ceil(4800 / (self.num_train_batches * self.batch_size))
self.scheduling_criterion = 'donuts_ap'
self.scheduling_mode = 'min' if "loss" in self.scheduling_criterion else 'max'
- self.weight_decay = 3e-5
- self.exclude_from_wd = []
- self.clip_norm = None
+ self.weight_decay = 1e-5
+ self.exclude_from_wd = ["norm"]
+ self.clip_norm = 200
#########################
# Testing / Plotting #
#########################
-
+ self.eval_test_fold_wise = True
# set the top-n-epochs to be saved for temporal averaging in testing.
self.save_n_models = 5
self.test_n_epochs = 5
self.test_aug_axes = (0, 1, (0, 1))
self.n_test_plots = 2
self.clustering = "wbc"
self.clustering_iou = 1e-5
# set a minimum epoch number for saving in case of instabilities in the first phase of training.
self.min_save_thresh = 0 if self.dim == 2 else 0
self.report_score_level = ['patient', 'rois'] # choose list from 'patient', 'rois'
self.class_labels = [Label(0, 'bg', (*self.white, 0.)),
Label(1, 'circles', (*self.orange, .9)),
Label(2, 'donuts', (*self.blue, .9)),]
if self.class_specific_seg:
self.seg_labels = self.class_labels
self.box_type2label = {label.name: label for label in self.box_labels}
self.class_id2label = {label.id: label for label in self.class_labels}
self.class_dict = {label.id: label.name for label in self.class_labels if label.id != 0}
self.seg_id2label = {label.id: label for label in self.seg_labels}
self.cmap = {label.id: label.color for label in self.seg_labels}
self.metrics = ["ap", "auc", "dice"]
self.patient_class_of_interest = 2 # patient metrics are only plotted for one class.
self.ap_match_ious = [0.1] # list of ious to be evaluated for ap-scoring.
self.model_selection_criteria = {name + "_ap": 1. for name in self.class_dict.values()}# criteria to average over for saving epochs.
self.min_det_thresh = 0.1 # minimum confidence value to select predictions for evaluation.
self.plot_prediction_histograms = True
self.plot_stat_curves = False
self.plot_class_ids = True
#########################
# Data Augmentation #
#########################
self.do_aug = False
self.da_kwargs={
'do_elastic_deform': True,
'alpha':(0., 1500.),
'sigma':(30., 50.),
'do_rotation':True,
'angle_x': (0., 2 * np.pi),
'angle_y': (0., 0),
'angle_z': (0., 0),
'do_scale': True,
'scale':(0.8, 1.1),
'random_crop':False,
'rand_crop_dist': (self.patch_size[0] / 2. - 3, self.patch_size[1] / 2. - 3),
'border_mode_data': 'constant',
'border_cval_data': 0,
'order_data': 1
}
if self.dim == 3:
self.da_kwargs['do_elastic_deform'] = False
self.da_kwargs['angle_x'] = (0, 0.0)
self.da_kwargs['angle_y'] = (0, 0.0) #must be 0!!
self.da_kwargs['angle_z'] = (0., 2 * np.pi)
#########################
# Add model specifics #
#########################
{'detection_fpn': self.add_det_fpn_configs,
'mrcnn': self.add_mrcnn_configs,
'retina_net': self.add_mrcnn_configs,
'retina_unet': self.add_mrcnn_configs,
}[self.model]()
def add_det_fpn_configs(self):
self.learning_rate = [3 * 1e-4] * self.num_epochs
self.dynamic_lr_scheduling = True
self.scheduling_criterion = 'torch_loss'
self.scheduling_mode = 'min' if "loss" in self.scheduling_criterion else 'max'
self.n_roi_candidates = 4 if self.dim == 2 else 6
# max number of roi candidates to identify per image (slice in 2D, volume in 3D)
# loss mode: either weighted cross entropy ('wce'), batch-wise dice loss ('dice), or the sum of both ('dice_wce')
self.seg_loss_mode = 'wce'
self.wce_weights = [1] * self.num_seg_classes if 'dice' in self.seg_loss_mode else [0.1, 1., 1.]
self.fp_dice_weight = 1 if self.dim == 2 else 1
# if <1, false positive predictions in foreground are penalized less.
self.detection_min_confidence = 0.05
# how to determine score of roi: 'max' or 'median'
self.score_det = 'max'
def add_mrcnn_configs(self):
# learning rate is a list with one entry per epoch.
self.learning_rate = [3e-4] * self.num_epochs
# disable mask head loss. (e.g. if no pixelwise annotations available)
self.frcnn_mode = False
# disable the re-sampling of mask proposals to original size for speed-up.
# since evaluation is detection-driven (box-matching) and not instance segmentation-driven (iou-matching),
# mask-outputs are optional.
self.return_masks_in_val = True
self.return_masks_in_test = False
# set number of proposal boxes to plot after each epoch.
- self.n_plot_rpn_props = 0 if self.dim == 2 else 0
+ self.n_plot_rpn_props = 2 if self.dim == 2 else 2
# number of classes for head networks: n_foreground_classes + 1 (background)
self.head_classes = self.num_classes + 1
# feature map strides per pyramid level are inferred from architecture.
self.backbone_strides = {'xy': [4, 8, 16, 32], 'z': [1, 2, 4, 8]}
# anchor scales are chosen according to expected object sizes in data set. Default uses only one anchor scale
# per pyramid level. (outer list are pyramid levels (corresponding to BACKBONE_STRIDES), inner list are scales per level.)
self.rpn_anchor_scales = {'xy': [[8], [16], [32], [64]], 'z': [[2], [4], [8], [16]]}
# choose which pyramid levels to extract features from: P2: 0, P3: 1, P4: 2, P5: 3.
self.pyramid_levels = [0, 1, 2, 3]
# number of feature maps in rpn. typically lowered in 3D to save gpu-memory.
self.n_rpn_features = 512 if self.dim == 2 else 128
# anchor ratios and strides per position in feature maps.
self.rpn_anchor_ratios = [0.5, 1., 2.]
self.rpn_anchor_stride = 1
# Threshold for first stage (RPN) non-maximum suppression (NMS): LOWER == HARDER SELECTION
self.rpn_nms_threshold = 0.7 if self.dim == 2 else 0.7
# loss sampling settings.
- self.rpn_train_anchors_per_image = 32 #per batch element
+ self.rpn_train_anchors_per_image = 64 #per batch element
self.train_rois_per_image = 2 #per batch element
self.roi_positive_ratio = 0.5
self.anchor_matching_iou = 0.7
# factor of top-k candidates to draw from per negative sample (stochastic-hard-example-mining).
# poolsize to draw top-k candidates from will be shem_poolsize * n_negative_samples.
- self.shem_poolsize = 10
+ self.shem_poolsize = 4
self.pool_size = (7, 7) if self.dim == 2 else (7, 7, 3)
self.mask_pool_size = (14, 14) if self.dim == 2 else (14, 14, 5)
self.mask_shape = (28, 28) if self.dim == 2 else (28, 28, 10)
self.rpn_bbox_std_dev = np.array([0.1, 0.1, 0.1, 0.2, 0.2, 0.2])
self.bbox_std_dev = np.array([0.1, 0.1, 0.1, 0.2, 0.2, 0.2])
self.window = np.array([0, 0, self.patch_size[0], self.patch_size[1]])
self.scale = np.array([self.patch_size[0], self.patch_size[1], self.patch_size[0], self.patch_size[1]])
if self.dim == 2:
self.rpn_bbox_std_dev = self.rpn_bbox_std_dev[:4]
self.bbox_std_dev = self.bbox_std_dev[:4]
self.window = self.window[:4]
self.scale = self.scale[:4]
# pre-selection in proposal-layer (stage 1) for NMS-speedup. applied per batch element.
self.pre_nms_limit = 3000 if self.dim == 2 else 6000
# n_proposals to be selected after NMS per batch element. too high numbers blow up memory if "detect_while_training" is True,
# since proposals of the entire batch are forwarded through second stage in as one "batch".
self.roi_chunk_size = 800 if self.dim == 2 else 600
self.post_nms_rois_training = 500 if self.dim == 2 else 75
self.post_nms_rois_inference = 500
# Final selection of detections (refine_detections)
self.model_max_instances_per_batch_element = 10 if self.dim == 2 else 30 # per batch element and class.
self.detection_nms_threshold = 1e-5 # needs to be > 0, otherwise all predictions are one cluster.
self.model_min_confidence = 0.1
if self.dim == 2:
self.backbone_shapes = np.array(
[[int(np.ceil(self.patch_size[0] / stride)),
int(np.ceil(self.patch_size[1] / stride))]
for stride in self.backbone_strides['xy']])
else:
self.backbone_shapes = np.array(
[[int(np.ceil(self.patch_size[0] / stride)),
int(np.ceil(self.patch_size[1] / stride)),
int(np.ceil(self.patch_size[2] / stride_z))]
for stride, stride_z in zip(self.backbone_strides['xy'], self.backbone_strides['z']
)])
if self.model == 'retina_net' or self.model == 'retina_unet':
# whether to use focal loss or SHEM for loss-sample selection
self.focal_loss = False
# implement extra anchor-scales according to retina-net publication.
self.rpn_anchor_scales['xy'] = [[ii[0], ii[0] * (2 ** (1 / 3)), ii[0] * (2 ** (2 / 3))] for ii in
self.rpn_anchor_scales['xy']]
self.rpn_anchor_scales['z'] = [[ii[0], ii[0] * (2 ** (1 / 3)), ii[0] * (2 ** (2 / 3))] for ii in
self.rpn_anchor_scales['z']]
self.n_anchors_per_pos = len(self.rpn_anchor_ratios) * 3
self.n_rpn_features = 256 if self.dim == 2 else 64
# pre-selection of detections for NMS-speedup. per entire batch.
self.pre_nms_limit = 10000 if self.dim == 2 else 50000
# anchor matching iou is lower than in Mask R-CNN according to https://arxiv.org/abs/1708.02002
self.anchor_matching_iou = 0.5
if self.model == 'retina_unet':
self.operate_stride1 = True
diff --git a/evaluator.py b/evaluator.py
index 2ac7fc6..26d2be5 100644
--- a/evaluator.py
+++ b/evaluator.py
@@ -1,980 +1,983 @@
#!/usr/bin/env python
# Copyright 2019 Division of Medical Image Computing, German Cancer Research Center (DKFZ).
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
# ==============================================================================
import os
from multiprocessing import Pool
import pickle
import time
import numpy as np
import pandas as pd
from sklearn.metrics import roc_auc_score, average_precision_score
from sklearn.metrics import roc_curve, precision_recall_curve
from sklearn.metrics import mean_squared_error, mean_absolute_error, accuracy_score
import torch
import utils.model_utils as mutils
import plotting as plg
import warnings
def get_roi_ap_from_df(inputs):
'''
:param df: data frame.
:param det_thresh: min_threshold for filtering out low confidence predictions.
:param per_patient_ap: boolean flag. evaluate average precision per patient id and average over per-pid results,
instead of computing one ap over whole data set.
:return: average_precision (float)
'''
df, det_thresh, per_patient_ap = inputs
if per_patient_ap:
pids_list = df.pid.unique()
aps = []
for match_iou in df.match_iou.unique():
iou_df = df[df.match_iou == match_iou]
for pid in pids_list:
pid_df = iou_df[iou_df.pid == pid]
all_p = len(pid_df[pid_df.class_label == 1])
pid_df = pid_df[(pid_df.det_type == 'det_fp') | (pid_df.det_type == 'det_tp')].sort_values('pred_score', ascending=False)
pid_df = pid_df[pid_df.pred_score > det_thresh]
if (len(pid_df) ==0 and all_p == 0):
pass
elif (len(pid_df) > 0 and all_p == 0):
aps.append(0)
else:
aps.append(compute_roi_ap(pid_df, all_p))
return np.mean(aps)
else:
aps = []
for match_iou in df.match_iou.unique():
iou_df = df[df.match_iou == match_iou]
# it's important to not apply the threshold before counting all_p in order to not lose the fn!
all_p = len(iou_df[(iou_df.det_type == 'det_tp') | (iou_df.det_type == 'det_fn')])
# sorting out all entries that are not fp or tp or have confidence(=pred_score) <= detection_threshold
iou_df = iou_df[(iou_df.det_type == 'det_fp') | (iou_df.det_type == 'det_tp')].sort_values('pred_score', ascending=False)
iou_df = iou_df[iou_df.pred_score > det_thresh]
if all_p>0:
aps.append(compute_roi_ap(iou_df, all_p))
return np.mean(aps)
def compute_roi_ap(df, all_p):
"""
adapted from: https://github.com/cocodataset/cocoapi/blob/master/PythonAPI/pycocotools/cocoeval.py
:param df: dataframe containing class labels of predictions sorted in descending manner by their prediction score.
:param all_p: number of all ground truth objects. (for denominator of recall.)
:return:
"""
tp = df.class_label.values
fp = (tp == 0) * 1
#recall thresholds, where precision will be measured
R = np.linspace(0., 1., np.round((1. - 0.) / .01).astype(int) + 1, endpoint=True)
tp_sum = np.cumsum(tp)
fp_sum = np.cumsum(fp)
n_dets = len(tp)
rc = tp_sum / all_p
pr = tp_sum / (fp_sum + tp_sum)
# initialize precision array over recall steps (q=queries).
q = [0. for _ in range(len(R))]
# numpy is slow without cython optimization for accessing elements
# python array gets significant speed improvement
pr = pr.tolist()
for i in range(n_dets - 1, 0, -1):
if pr[i] > pr[i - 1]:
pr[i - 1] = pr[i]
#--> pr[i]<=pr[i-1] for all i since we want to consider the maximum
#precision value for a queried interval
# discretize empiric recall steps with given bins.
assert np.all(rc[:-1]<=rc[1:]), "recall not sorted ascendingly"
inds = np.searchsorted(rc, R, side='left')
try:
for rc_ix, pr_ix in enumerate(inds):
q[rc_ix] = pr[pr_ix]
except IndexError: #now q is filled with pr values up to first non-available index
pass
return np.mean(q)
def roi_avp(inputs):
'''
:param df: data frame.
:param det_thresh: min_threshold for filtering out low confidence predictions.
:param per_patient_ap: boolean flag. evaluate average precision per patient id and average over per-pid results,
instead of computing one ap over whole data set.
:return: average_precision (float)
'''
df, det_thresh, per_patient_ap = inputs
if per_patient_ap:
pids_list = df.pid.unique()
aps = []
for match_iou in df.match_iou.unique():
iou_df = df[df.match_iou == match_iou]
for pid in pids_list:
pid_df = iou_df[iou_df.pid == pid]
all_p = len(pid_df[pid_df.class_label == 1])
mask = ((pid_df.rg_bins == pid_df.rg_bin_target) & (pid_df.det_type == 'det_tp')) | (pid_df.det_type == 'det_fp')
pid_df = pid_df[mask].sort_values('pred_score', ascending=False)
pid_df = pid_df[pid_df.pred_score > det_thresh]
if (len(pid_df) ==0 and all_p == 0):
pass
elif (len(pid_df) > 0 and all_p == 0):
aps.append(0)
else:
aps.append(compute_roi_ap(pid_df, all_p))
return np.mean(aps)
else:
aps = []
for match_iou in df.match_iou.unique():
iou_df = df[df.match_iou == match_iou]
#it's important to not apply the threshold before counting all_positives!
all_p = len(iou_df[(iou_df.det_type == 'det_tp') | (iou_df.det_type == 'det_fn')])
# filtering out tps which don't match rg_bin target at this point is same as reclassifying them as fn.
# also sorting out all entries that are not fp or have confidence(=pred_score) <= detection_threshold
mask = ((iou_df.rg_bins == iou_df.rg_bin_target) & (iou_df.det_type == 'det_tp')) | (iou_df.det_type == 'det_fp')
iou_df = iou_df[mask].sort_values('pred_score', ascending=False)
iou_df = iou_df[iou_df.pred_score > det_thresh]
if all_p>0:
aps.append(compute_roi_ap(iou_df, all_p))
return np.mean(aps)
def compute_prc(df):
"""compute precision-recall curve with maximum precision per recall interval.
:param df:
:param all_p: # of all positive samples in data.
:return: array: [precisions, recall query values]
"""
assert (df.class_label==1).any(), "cannot compute prc when no positives in data."
all_p = len(df[(df.det_type == 'det_tp') | (df.det_type == 'det_fn')])
df = df[(df.det_type=="det_tp") | (df.det_type=="det_fp")]
df = df.sort_values("pred_score", ascending=False)
# recall thresholds, where precision will be measured
scores = df.pred_score.values
labels = df.class_label.values
n_dets = len(scores)
pr = np.zeros((n_dets,))
rc = pr.copy()
for rank in range(n_dets):
tp = np.count_nonzero(labels[:rank+1]==1)
fp = np.count_nonzero(labels[:rank+1]==0)
pr[rank] = tp/(tp+fp)
rc[rank] = tp/all_p
#after obj detection convention/ coco-dataset template: take maximum pr within intervals:
# --> pr[i]<=pr[i-1] for all i since we want to consider the maximum
# precision value for a queried interval
for i in range(n_dets - 1, 0, -1):
if pr[i] > pr[i - 1]:
pr[i - 1] = pr[i]
R = np.linspace(0., 1., np.round((1. - 0.) / .01).astype(int) + 1, endpoint=True)#precision queried at R points
inds = np.searchsorted(rc, R, side='left')
queries = np.zeros((len(R),))
try:
for q_ix, rank in enumerate(inds):
queries[q_ix] = pr[rank]
except IndexError:
pass
return np.array((queries, R))
def RMSE(y_true, y_pred, weights=None):
if len(y_true)>0:
return np.sqrt(mean_squared_error(y_true, y_pred, sample_weight=weights))
else:
return np.nan
def MAE_w_std(y_true, y_pred, weights=None):
if len(y_true)>0:
y_true, y_pred = np.array(y_true), np.array(y_pred)
deltas = np.abs(y_true-y_pred)
mae = np.average(deltas, weights=weights, axis=0).item()
skmae = mean_absolute_error(y_true, y_pred, sample_weight=weights)
assert np.allclose(mae, skmae, atol=1e-6), "mae {}, sklearn mae {}".format(mae, skmae)
std = np.std(weights*deltas)
return mae, std
else:
return np.nan, np.nan
def MAE(y_true, y_pred, weights=None):
if len(y_true)>0:
return mean_absolute_error(y_true, y_pred, sample_weight=weights)
else:
return np.nan
def accuracy(y_true, y_pred, weights=None):
if len(y_true)>0:
return accuracy_score(y_true, y_pred, sample_weight=weights)
else:
return np.nan
# noinspection PyCallingNonCallable
class Evaluator():
""" Evaluates given results dicts. Can return results as updated monitor_metrics. Can save test data frames to
file.
"""
def __init__(self, cf, logger, mode='test'):
"""
:param mode: either 'train', 'val_sampling', 'val_patient' or 'test'. handles prediction lists of different forms.
"""
self.cf = cf
self.logger = logger
self.mode = mode
self.regress_flag = any(['regression' in task for task in self.cf.prediction_tasks])
self.plot_dir = self.cf.test_dir if self.mode == "test" else self.cf.plot_dir
if self.cf.plot_prediction_histograms:
self.hist_dir = os.path.join(self.plot_dir, 'histograms')
os.makedirs(self.hist_dir, exist_ok=True)
if self.cf.plot_stat_curves:
self.curves_dir = os.path.join(self.plot_dir, 'stat_curves')
os.makedirs(self.curves_dir, exist_ok=True)
def eval_losses(self, batch_res_dicts):
if hasattr(self.cf, "losses_to_monitor"):
loss_names = self.cf.losses_to_monitor
else:
loss_names = {name for b_res_dict in batch_res_dicts for name in b_res_dict if 'loss' in name}
self.epoch_losses = {l_name: torch.tensor([b_res_dict[l_name] for b_res_dict in batch_res_dicts if l_name
in b_res_dict.keys()]).mean().item() for l_name in loss_names}
def eval_segmentations(self, batch_res_dicts, pid_list):
batch_dices = [b_res_dict['batch_dices'] for b_res_dict in batch_res_dicts if
'batch_dices' in b_res_dict.keys()] # shape (n_batches, n_seg_classes)
if len(batch_dices) > 0:
batch_dices = np.array(batch_dices) # dims n_batches x 1 in sampling / n_test_epochs x n_classes
assert batch_dices.shape[1] == self.cf.num_seg_classes, "bdices shp {}, n seg cl {}, pid lst len {}".format(
batch_dices.shape, self.cf.num_seg_classes, len(pid_list))
self.seg_df = pd.DataFrame()
for seg_id in range(batch_dices.shape[1]):
self.seg_df[self.cf.seg_id2label[seg_id].name + "_dice"] = batch_dices[:,
seg_id] # one row== one batch, one column== one class
# self.seg_df[self.cf.seg_id2label[seg_id].name+"_dice"] = np.concatenate(batch_dices[:,:,seg_id])
self.seg_df['fold'] = self.cf.fold
if self.mode == "val_patient" or self.mode == "test":
# need to make it more conform between sampling and patient-mode
self.seg_df["pid"] = [pid for pix, pid in enumerate(pid_list)] # for b_inst in batch_inst_boxes[pix]]
else:
self.seg_df["pid"] = np.nan
def eval_boxes(self, batch_res_dicts, pid_list, obj_cl_dict,
obj_cl_identifiers={"gt":'class_targets', "pred":'box_pred_class_id'}):
"""
:param batch_res_dicts:
:param pid_list: [pid_0, pid_1, ...]
:return:
"""
if self.mode == 'train' or self.mode == 'val_sampling':
# one pid per batch element
# batch_size > 1, with varying patients across batch:
# [[[results_0, ...], [pid_0, ...]], [[results_n, ...], [pid_n, ...]], ...]
# -> [results_0, results_1, ..]
batch_inst_boxes = [b_res_dict['boxes'] for b_res_dict in batch_res_dicts] # len: nr of batches in epoch
batch_inst_boxes = [[b_inst_boxes] for whole_batch_boxes in batch_inst_boxes for b_inst_boxes in
whole_batch_boxes] # len: batch instances of whole epoch
assert np.all(len(b_boxes_list) == self.cf.batch_size for b_boxes_list in batch_inst_boxes)
elif self.mode == "val_patient" or self.mode == "test":
# patient processing, one element per batch = one patient.
# [[results_0, pid_0], [results_1, pid_1], ...] -> [results_0, results_1, ..]
# in patientbatchiterator there is only one pid per batch
batch_inst_boxes = [b_res_dict['boxes'] for b_res_dict in batch_res_dicts]
# in patient mode not actually per batch instance, but per whole batch!
if hasattr(self.cf, "eval_test_separately") and self.cf.eval_test_separately:
""" you could write your own routines to add GTs to raw predictions for evaluation.
implemented standard is: cf.eval_test_separately = False or not set --> GTs are saved at same time
and in same file as raw prediction results.
"""
raise NotImplementedError
assert len(batch_inst_boxes) == len(pid_list)
df_list_preds = []
df_list_labels = []
df_list_class_preds = []
df_list_pids = []
df_list_type = []
df_list_match_iou = []
df_list_n_missing = []
df_list_regressions = []
df_list_rg_targets = []
df_list_rg_bins = []
df_list_rg_bin_targets = []
df_list_rg_uncs = []
for match_iou in self.cf.ap_match_ious:
self.logger.info('evaluating with ap_match_iou: {}'.format(match_iou))
for cl in list(obj_cl_dict.keys()):
for pix, pid in enumerate(pid_list):
len_df_list_before_patient = len(df_list_pids)
# input of each batch element is a list of boxes, where each box is a dictionary.
for b_inst_ix, b_boxes_list in enumerate(batch_inst_boxes[pix]):
b_tar_boxes = []
b_cand_boxes, b_cand_scores, b_cand_n_missing = [], [], []
if self.regress_flag:
b_tar_regs, b_tar_rg_bins = [], []
b_cand_regs, b_cand_rg_bins, b_cand_rg_uncs = [], [], []
for box in b_boxes_list:
# each box is either gt or detection or proposal/anchor
# we need all gts in the same order & all dets in same order
if box['box_type'] == 'gt' and box[obj_cl_identifiers["gt"]] == cl:
b_tar_boxes.append(box["box_coords"])
if self.regress_flag:
b_tar_regs.append(np.array(box['regression_targets'], dtype='float32'))
b_tar_rg_bins.append(box['rg_bin_targets'])
if box['box_type'] == 'det' and box[obj_cl_identifiers["pred"]] == cl:
b_cand_boxes.append(box["box_coords"])
b_cand_scores.append(box["box_score"])
b_cand_n_missing.append(box["cluster_n_missing"] if 'cluster_n_missing' in box.keys() else np.nan)
if self.regress_flag:
b_cand_regs.append(box["regression"])
b_cand_rg_bins.append(box["rg_bin"])
b_cand_rg_uncs.append(box["rg_uncertainty"] if 'rg_uncertainty' in box.keys() else np.nan)
b_tar_boxes = np.array(b_tar_boxes)
b_cand_boxes, b_cand_scores, b_cand_n_missing = np.array(b_cand_boxes), np.array(b_cand_scores), np.array(b_cand_n_missing)
if self.regress_flag:
b_tar_regs, b_tar_rg_bins = np.array(b_tar_regs), np.array(b_tar_rg_bins)
b_cand_regs, b_cand_rg_bins, b_cand_rg_uncs = np.array(b_cand_regs), np.array(b_cand_rg_bins), np.array(b_cand_rg_uncs)
# check if predictions and ground truth boxes exist and match them according to match_iou.
if not 0 in b_cand_boxes.shape and not 0 in b_tar_boxes.shape:
assert np.all(np.round(b_cand_scores,6) <= 1.), "there is a box score>1: {}".format(b_cand_scores[~(b_cand_scores<=1.)])
#coords_check = np.array([len(coords)==self.cf.dim*2 for coords in b_cand_boxes])
#assert np.all(coords_check), "cand box with wrong bcoords dim: {}, mode: {}".format(b_cand_boxes[~coords_check], self.mode)
expected_dim = len(b_cand_boxes[0])
assert np.all([len(coords) == expected_dim for coords in b_tar_boxes]), \
"gt/cand box coords mismatch, expected dim: {}.".format(expected_dim)
# overlaps: shape len(cand_boxes) x len(tar_boxes)
overlaps = mutils.compute_overlaps(b_cand_boxes, b_tar_boxes)
# match_cand_ixs: shape (nr_of_matches,)
# theses indices are the indices of b_cand_boxes
match_cand_ixs = np.argwhere(np.max(overlaps, axis=1) > match_iou)[:, 0]
non_match_cand_ixs = np.argwhere(np.max(overlaps, 1) <= match_iou)[:, 0]
# the corresponding gt assigned to the pred boxes by highest iou overlap,
# i.e., match_gt_ixs holds index into b_tar_boxes for each entry in match_cand_ixs,
# i.e., gt_ixs and cand_ixs are paired via their position in their list
# (cand_ixs[j] corresponds to gt_ixs[j])
match_gt_ixs = np.argmax(overlaps[match_cand_ixs, :], axis=1) if \
not 0 in match_cand_ixs.shape else np.array([])
assert len(match_gt_ixs)==len(match_cand_ixs)
#match_gt_ixs: shape (nr_of_matches,) or 0
non_match_gt_ixs = np.array(
[ii for ii in np.arange(b_tar_boxes.shape[0]) if ii not in match_gt_ixs])
unique, counts = np.unique(match_gt_ixs, return_counts=True)
# check for double assignments, i.e. two predictions having been assigned to the same gt.
# according to the COCO-metrics, only one prediction counts as true positive, the rest counts as
# false positive. This case is supposed to be avoided by the model itself by,
# e.g. using a low enough NMS threshold.
if np.any(counts > 1):
double_match_gt_ixs = unique[np.argwhere(counts > 1)[:, 0]]
keep_max = []
double_match_list = []
for dg in double_match_gt_ixs:
double_match_cand_ixs = match_cand_ixs[np.argwhere(match_gt_ixs == dg)]
keep_max.append(double_match_cand_ixs[np.argmax(b_cand_scores[double_match_cand_ixs])])
double_match_list += [ii for ii in double_match_cand_ixs]
fp_ixs = np.array([ii for ii in match_cand_ixs if
(ii in double_match_list and ii not in keep_max)])
# count as fp: boxes that match gt above match_iou threshold but have not highest class confidence score
match_gt_ixs = np.array([gt_ix for ii, gt_ix in enumerate(match_gt_ixs) if match_cand_ixs[ii] not in fp_ixs])
match_cand_ixs = np.array([cand_ix for cand_ix in match_cand_ixs if cand_ix not in fp_ixs])
assert len(match_gt_ixs) == len(match_cand_ixs)
df_list_preds += [ii for ii in b_cand_scores[fp_ixs]]
df_list_labels += [0] * fp_ixs.shape[0] # means label==gt==0==bg for all these fp_ixs
df_list_class_preds += [cl] * fp_ixs.shape[0]
df_list_n_missing += [n for n in b_cand_n_missing[fp_ixs]]
if self.regress_flag:
df_list_regressions += [r for r in b_cand_regs[fp_ixs]]
df_list_rg_bins += [r for r in b_cand_rg_bins[fp_ixs]]
df_list_rg_uncs += [r for r in b_cand_rg_uncs[fp_ixs]]
df_list_rg_targets += [[0.]*self.cf.regression_n_features] * fp_ixs.shape[0]
df_list_rg_bin_targets += [0.] * fp_ixs.shape[0]
df_list_pids += [pid] * fp_ixs.shape[0]
df_list_type += ['det_fp'] * fp_ixs.shape[0]
# matched/tp:
if not 0 in match_cand_ixs.shape:
df_list_preds += list(b_cand_scores[match_cand_ixs])
df_list_labels += [1] * match_cand_ixs.shape[0]
df_list_class_preds += [cl] * match_cand_ixs.shape[0]
df_list_n_missing += list(b_cand_n_missing[match_cand_ixs])
if self.regress_flag:
df_list_regressions += list(b_cand_regs[match_cand_ixs])
df_list_rg_bins += list(b_cand_rg_bins[match_cand_ixs])
df_list_rg_uncs += list(b_cand_rg_uncs[match_cand_ixs])
assert len(match_cand_ixs)==len(match_gt_ixs)
df_list_rg_targets += list(b_tar_regs[match_gt_ixs])
df_list_rg_bin_targets += list(b_tar_rg_bins[match_gt_ixs])
df_list_pids += [pid] * match_cand_ixs.shape[0]
df_list_type += ['det_tp'] * match_cand_ixs.shape[0]
# rest fp:
if not 0 in non_match_cand_ixs.shape:
df_list_preds += list(b_cand_scores[non_match_cand_ixs])
df_list_labels += [0] * non_match_cand_ixs.shape[0]
df_list_class_preds += [cl] * non_match_cand_ixs.shape[0]
df_list_n_missing += list(b_cand_n_missing[non_match_cand_ixs])
if self.regress_flag:
df_list_regressions += list(b_cand_regs[non_match_cand_ixs])
df_list_rg_bins += list(b_cand_rg_bins[non_match_cand_ixs])
df_list_rg_uncs += list(b_cand_rg_uncs[non_match_cand_ixs])
df_list_rg_targets += [[0.]*self.cf.regression_n_features] * non_match_cand_ixs.shape[0]
df_list_rg_bin_targets += [0.] * non_match_cand_ixs.shape[0]
df_list_pids += [pid] * non_match_cand_ixs.shape[0]
df_list_type += ['det_fp'] * non_match_cand_ixs.shape[0]
# fn:
if not 0 in non_match_gt_ixs.shape:
df_list_preds += [0] * non_match_gt_ixs.shape[0]
df_list_labels += [1] * non_match_gt_ixs.shape[0]
df_list_class_preds += [cl] * non_match_gt_ixs.shape[0]
df_list_n_missing += [np.nan] * non_match_gt_ixs.shape[0]
if self.regress_flag:
df_list_regressions += [[0.]*self.cf.regression_n_features] * non_match_gt_ixs.shape[0]
df_list_rg_bins += [0.] * non_match_gt_ixs.shape[0]
df_list_rg_uncs += [np.nan] * non_match_gt_ixs.shape[0]
df_list_rg_targets += list(b_tar_regs[non_match_gt_ixs])
df_list_rg_bin_targets += list(b_tar_rg_bins[non_match_gt_ixs])
df_list_pids += [pid] * non_match_gt_ixs.shape[0]
df_list_type += ['det_fn'] * non_match_gt_ixs.shape[0]
# only fp:
if not 0 in b_cand_boxes.shape and 0 in b_tar_boxes.shape:
# means there is no gt in all samples! any preds have to be fp.
df_list_preds += list(b_cand_scores)
df_list_labels += [0] * b_cand_boxes.shape[0]
df_list_class_preds += [cl] * b_cand_boxes.shape[0]
df_list_n_missing += list(b_cand_n_missing)
if self.regress_flag:
df_list_regressions += list(b_cand_regs)
df_list_rg_bins += list(b_cand_rg_bins)
df_list_rg_uncs += list(b_cand_rg_uncs)
df_list_rg_targets += [[0.]*self.cf.regression_n_features] * b_cand_boxes.shape[0]
df_list_rg_bin_targets += [0.] * b_cand_boxes.shape[0]
df_list_pids += [pid] * b_cand_boxes.shape[0]
df_list_type += ['det_fp'] * b_cand_boxes.shape[0]
# only fn:
if 0 in b_cand_boxes.shape and not 0 in b_tar_boxes.shape:
df_list_preds += [0] * b_tar_boxes.shape[0]
df_list_labels += [1] * b_tar_boxes.shape[0]
df_list_class_preds += [cl] * b_tar_boxes.shape[0]
df_list_n_missing += [np.nan] * b_tar_boxes.shape[0]
if self.regress_flag:
df_list_regressions += [[0.]*self.cf.regression_n_features] * b_tar_boxes.shape[0]
df_list_rg_bins += [0.] * b_tar_boxes.shape[0]
df_list_rg_uncs += [np.nan] * b_tar_boxes.shape[0]
df_list_rg_targets += list(b_tar_regs)
df_list_rg_bin_targets += list(b_tar_rg_bins)
df_list_pids += [pid] * b_tar_boxes.shape[0]
df_list_type += ['det_fn'] * b_tar_boxes.shape[0]
# empty patient with 0 detections needs empty patient score, in order to not disappear from stats.
# filtered out for roi-level evaluation later. During training (and val_sampling),
# tn are assigned per sample independently of associated patients.
# i.e., patient_tn is also meant as sample_tn if a list of samples is evaluated instead of whole patient
if len(df_list_pids) == len_df_list_before_patient:
df_list_preds += [0]
df_list_labels += [0]
df_list_class_preds += [cl]
df_list_n_missing += [np.nan]
if self.regress_flag:
df_list_regressions += [[0.]*self.cf.regression_n_features]
df_list_rg_bins += [0.]
df_list_rg_uncs += [np.nan]
df_list_rg_targets += [[0.]*self.cf.regression_n_features]
df_list_rg_bin_targets += [0.]
df_list_pids += [pid]
df_list_type += ['patient_tn'] # true negative: no ground truth boxes, no detections.
df_list_match_iou += [match_iou] * (len(df_list_preds) - len(df_list_match_iou))
self.test_df = pd.DataFrame()
self.test_df['pred_score'] = df_list_preds
self.test_df['class_label'] = df_list_labels
# class labels are gt, 0,1, only indicate neg/pos (or bg/fg) remapped from all classes
self.test_df['pred_class'] = df_list_class_preds # can be diff than 0,1
self.test_df['pid'] = df_list_pids
self.test_df['det_type'] = df_list_type
self.test_df['fold'] = self.cf.fold
self.test_df['match_iou'] = df_list_match_iou
self.test_df['cluster_n_missing'] = df_list_n_missing
if self.regress_flag:
self.test_df['regressions'] = df_list_regressions
self.test_df['rg_targets'] = df_list_rg_targets
self.test_df['rg_uncertainties'] = df_list_rg_uncs
self.test_df['rg_bins'] = df_list_rg_bins
# super weird error: pandas does not properly add an attribute if column is named "rg_bin_targets" ... ?!?
self.test_df['rg_bin_target'] = df_list_rg_bin_targets
assert hasattr(self.test_df, "rg_bin_target")
#fn_df = self.test_df[self.test_df["det_type"] == "det_fn"]
pass
def evaluate_predictions(self, results_list, monitor_metrics=None):
"""
Performs the matching of predicted boxes and ground truth boxes. Loops over list of matching IoUs and foreground classes.
Resulting info of each prediction is stored as one line in an internal dataframe, with the keys:
det_type: 'tp' (true positive), 'fp' (false positive), 'fn' (false negative), 'tn' (true negative)
pred_class: foreground class which the object predicts.
pid: corresponding patient-id.
pred_score: confidence score [0, 1]
fold: corresponding fold of CV.
match_iou: utilized IoU for matching.
:param results_list: list of model predictions. Either from train/val_sampling (patch processing) for monitoring with form:
[[[results_0, ...], [pid_0, ...]], [[results_n, ...], [pid_n, ...]], ...]
Or from val_patient/testing (patient processing), with form: [[results_0, pid_0], [results_1, pid_1], ...])
:param monitor_metrics (optional): dict of dicts with all metrics of previous epochs.
:return monitor_metrics: if provided (during training), return monitor_metrics now including results of current epoch.
"""
# gets results_list = [[batch_instances_box_lists], [batch_instances_pids]]*n_batches
# we want to evaluate one batch_instance (= 2D or 3D image) at a time.
self.logger.info('evaluating in mode {}'.format(self.mode))
batch_res_dicts = [batch[0] for batch in results_list] # len: nr of batches in epoch
if self.mode == 'train' or self.mode=='val_sampling':
# one pid per batch element
# [[[results_0, ...], [pid_0, ...]], [[results_n, ...], [pid_n, ...]], ...]
# -> [pid_0, pid_1, ...]
# additional list wrapping to make conform with below per-patient batches, where one pid is linked to more than one batch instance
pid_list = [batch_instance_pid for batch in results_list for batch_instance_pid in batch[1]]
elif self.mode == "val_patient" or self.mode=="test":
# [[results_0, pid_0], [results_1, pid_1], ...] -> [pid_0, pid_1, ...]
# in patientbatchiterator there is only one pid per batch
pid_list = [np.unique(batch[1]) for batch in results_list]
assert np.all([len(pid)==1 for pid in pid_list]), "pid list in patient-eval mode, should only contain a single scalar per patient: {}".format(pid_list)
pid_list = [pid[0] for pid in pid_list]
else:
raise Exception("undefined run mode encountered")
self.eval_losses(batch_res_dicts)
self.eval_segmentations(batch_res_dicts, pid_list)
self.eval_boxes(batch_res_dicts, pid_list, self.cf.class_dict)
if monitor_metrics is not None:
# return all_stats, updated monitor_metrics
return self.return_metrics(self.test_df, self.cf.class_dict, monitor_metrics)
def return_metrics(self, df, obj_cl_dict, monitor_metrics=None, boxes_only=False):
"""
Calculates metric scores for internal data frame. Called directly from evaluate_predictions during training for
monitoring, or from score_test_df during inference (for single folds or aggregated test set).
Loops over foreground classes and score_levels ('roi' and/or 'patient'), gets scores and stores them.
Optionally creates plots of prediction histograms and ROC/PR curves.
:param df: Data frame that holds evaluated predictions.
:param obj_cl_dict: Dict linking object-class ids to object-class names. E.g., {1: "bikes", 2 : "cars"}. Set in
configs as cf.class_dict.
:param monitor_metrics: dict of dicts with all metrics of previous epochs. This function adds metrics for
current epoch and returns the same object.
:param boxes_only: whether to produce metrics only for the boxes, not the segmentations.
:return: all_stats: list. Contains dicts with resulting scores for each combination of foreground class and
score_level.
:return: monitor_metrics
"""
# -------------- monitoring independent of class, score level ------------
if monitor_metrics is not None:
for l_name in self.epoch_losses:
monitor_metrics[l_name] = [self.epoch_losses[l_name]]
# -------------- metrics calc dependent on class, score level ------------
all_stats = [] # all_stats: one entry per score_level per class
for cl in list(obj_cl_dict.keys()): # bg eval is neglected
cl_name = obj_cl_dict[cl]
cl_df = df[df.pred_class == cl]
if hasattr(self, "seg_df") and not boxes_only:
dice_col = self.cf.seg_id2label[cl].name+"_dice"
seg_cl_df = self.seg_df.loc[:,['pid', dice_col, 'fold']]
for score_level in self.cf.report_score_level:
stats_dict = {}
stats_dict['name'] = 'fold_{} {} {}'.format(self.cf.fold, score_level, cl_name)
# -------------- RoI-based -----------------
if score_level == 'rois':
stats_dict['auc'] = np.nan
stats_dict['roc'] = np.nan
if monitor_metrics is not None:
tn = len(cl_df[cl_df.det_type == "patient_tn"])
tp = len(cl_df[(cl_df.det_type == "det_tp")&(cl_df.pred_score>self.cf.min_det_thresh)])
fp = len(cl_df[(cl_df.det_type == "det_fp")&(cl_df.pred_score>self.cf.min_det_thresh)])
fn = len(cl_df[cl_df.det_type == "det_fn"])
sens = np.divide(tp, (fn + tp))
monitor_metrics.update({"Bin_Stats/" + cl_name + "_fp": [fp], "Bin_Stats/" + cl_name + "_tp": [tp],
"Bin_Stats/" + cl_name + "_fn": [fn], "Bin_Stats/" + cl_name + "_tn": [tn],
"Bin_Stats/" + cl_name + "_sensitivity": [sens]})
# list wrapping only needed bc other metrics are recorded over all epochs;
spec_df = cl_df[cl_df.det_type != 'patient_tn']
if self.regress_flag:
# filter false negatives out for regression-only eval since regressor didn't predict
truncd_df = spec_df[(((spec_df.det_type == "det_fp") | (
spec_df.det_type == "det_tp")) & spec_df.pred_score > self.cf.min_det_thresh)]
truncd_df_tp = truncd_df[truncd_df.det_type == "det_tp"]
weights, weights_tp = truncd_df.pred_score.tolist(), truncd_df_tp.pred_score.tolist()
y_true, y_pred = truncd_df.rg_targets.tolist(), truncd_df.regressions.tolist()
stats_dict["rg_RMSE"] = RMSE(y_true, y_pred)
stats_dict["rg_MAE"] = MAE(y_true, y_pred)
stats_dict["rg_RMSE_weighted"] = RMSE(y_true, y_pred, weights)
stats_dict["rg_MAE_weighted"] = MAE(y_true, y_pred, weights)
y_true, y_pred = truncd_df_tp.rg_targets.tolist(), truncd_df_tp.regressions.tolist()
stats_dict["rg_MAE_weighted_tp"] = MAE(y_true, y_pred, weights_tp)
stats_dict["rg_MAE_w_std_weighted_tp"] = MAE_w_std(y_true, y_pred, weights_tp)
y_true, y_pred = truncd_df.rg_bin_target.tolist(), truncd_df.rg_bins.tolist()
stats_dict["rg_bin_accuracy"] = accuracy(y_true, y_pred)
stats_dict["rg_bin_accuracy_weighted"] = accuracy(y_true, y_pred, weights)
y_true, y_pred = truncd_df_tp.rg_bin_target.tolist(), truncd_df_tp.rg_bins.tolist()
stats_dict["rg_bin_accuracy_weighted_tp"] = accuracy(y_true, y_pred, weights_tp)
if np.any(~truncd_df.rg_uncertainties.isna()):
# det_fn are expected to be NaN so they drop out in means
stats_dict.update({"rg_uncertainty": truncd_df.rg_uncertainties.mean(),
"rg_uncertainty_tp": truncd_df_tp.rg_uncertainties.mean(),
"rg_uncertainty_tp_weighted": (truncd_df_tp.rg_uncertainties * truncd_df_tp.pred_score).sum()
/ truncd_df_tp.pred_score.sum()
})
if (spec_df.class_label==1).any():
stats_dict['ap'] = get_roi_ap_from_df((spec_df, self.cf.min_det_thresh, self.cf.per_patient_ap))
stats_dict['prc'] = precision_recall_curve(spec_df.class_label.tolist(), spec_df.pred_score.tolist())
if self.regress_flag:
stats_dict['avp'] = roi_avp((spec_df, self.cf.min_det_thresh, self.cf.per_patient_ap))
else:
stats_dict['ap'] = np.nan
stats_dict['prc'] = np.nan
stats_dict['avp'] = np.nan
# np.nan is formattable by __format__ as a float, None-type is not
if hasattr(self, "seg_df") and not boxes_only:
stats_dict["dice"] = seg_cl_df.loc[:,dice_col].mean() # mean per all rois in this epoch
stats_dict["dice_std"] = seg_cl_df.loc[:,dice_col].std()
# for the aggregated test set case, additionally get the scores of averaging over fold results.
if self.cf.evaluate_fold_means and len(df.fold.unique()) > 1:
aps = []
for fold in df.fold.unique():
fold_df = spec_df[spec_df.fold == fold]
if (fold_df.class_label==1).any():
aps.append(get_roi_ap_from_df((fold_df, self.cf.min_det_thresh, self.cf.per_patient_ap)))
stats_dict['ap_folds_mean'] = np.mean(aps) if len(aps)>0 else np.nan
stats_dict['ap_folds_std'] = np.std(aps) if len(aps)>0 else np.nan
stats_dict['auc_folds_mean'] = np.nan
stats_dict['auc_folds_std'] = np.nan
if self.regress_flag:
avps, accuracies, MAEs = [], [], []
for fold in df.fold.unique():
fold_df = spec_df[spec_df.fold == fold]
if (fold_df.class_label == 1).any():
avps.append(roi_avp((fold_df, self.cf.min_det_thresh, self.cf.per_patient_ap)))
truncd_df_tp = fold_df[((fold_df.det_type == "det_tp") & fold_df.pred_score > self.cf.min_det_thresh)]
weights_tp = truncd_df_tp.pred_score.tolist()
y_true, y_pred = truncd_df_tp.rg_bin_target.tolist(), truncd_df_tp.rg_bins.tolist()
accuracies.append(accuracy(y_true, y_pred, weights_tp))
y_true, y_pred = truncd_df_tp.rg_targets.tolist(), truncd_df_tp.regressions.tolist()
MAEs.append(MAE_w_std(y_true, y_pred, weights_tp))
stats_dict['avp_folds_mean'] = np.mean(avps) if len(avps) > 0 else np.nan
stats_dict['avp_folds_std'] = np.std(avps) if len(avps) > 0 else np.nan
stats_dict['rg_bin_accuracy_weighted_tp_folds_mean'] = np.mean(accuracies) if len(accuracies) > 0 else np.nan
stats_dict['rg_bin_accuracy_weighted_tp_folds_std'] = np.std(accuracies) if len(accuracies) > 0 else np.nan
stats_dict['rg_MAE_w_std_weighted_tp_folds_mean'] = np.mean(MAEs, axis=0) if len(MAEs) > 0 else np.nan
stats_dict['rg_MAE_w_std_weighted_tp_folds_std'] = np.std(MAEs, axis=0) if len(MAEs) > 0 else np.nan
if hasattr(self, "seg_df") and not boxes_only and self.cf.evaluate_fold_means and len(seg_cl_df.fold.unique()) > 1:
fold_means = seg_cl_df.groupby(['fold'], as_index=True).agg({dice_col:"mean"})
stats_dict["dice_folds_mean"] = float(fold_means.mean())
stats_dict["dice_folds_std"] = float(fold_means.std())
# -------------- patient-based -----------------
# on patient level, aggregate predictions per patient (pid): The patient predicted score is the highest
# confidence prediction for this class. The patient class label is 1 if roi of this class exists in patient, else 0.
if score_level == 'patient':
#this is the critical part in patient scoring: only the max gt and max pred score are taken per patient!
#--> does mix up values from separate detections
spec_df = cl_df.groupby(['pid'], as_index=False)
agg_args = {'class_label': 'max', 'pred_score': 'max', 'fold': 'first'}
if self.regress_flag:
# pandas throws error if aggregated value is np.array, not if is list.
agg_args.update({'regressions': lambda series: list(series.iloc[np.argmax(series.apply(np.linalg.norm).values)]),
'rg_targets': lambda series: list(series.iloc[np.argmax(series.apply(np.linalg.norm).values)]),
'rg_bins': 'max', 'rg_bin_target': 'max',
'rg_uncertainties': 'max'
})
if hasattr(cl_df, "cluster_n_missing"):
agg_args.update({'cluster_n_missing': 'mean'})
spec_df = spec_df.agg(agg_args)
if len(spec_df.class_label.unique()) > 1:
stats_dict['auc'] = roc_auc_score(spec_df.class_label.tolist(), spec_df.pred_score.tolist())
stats_dict['roc'] = roc_curve(spec_df.class_label.tolist(), spec_df.pred_score.tolist())
else:
stats_dict['auc'] = np.nan
stats_dict['roc'] = np.nan
if (spec_df.class_label == 1).any():
patient_cl_labels = spec_df.class_label.tolist()
stats_dict['ap'] = average_precision_score(patient_cl_labels, spec_df.pred_score.tolist())
stats_dict['prc'] = precision_recall_curve(patient_cl_labels, spec_df.pred_score.tolist())
if self.regress_flag:
avp_scores = spec_df[spec_df.rg_bins == spec_df.rg_bin_target].pred_score.tolist()
avp_scores += [0.] * (len(patient_cl_labels) - len(avp_scores))
stats_dict['avp'] = average_precision_score(patient_cl_labels, avp_scores)
else:
stats_dict['ap'] = np.nan
stats_dict['prc'] = np.nan
stats_dict['avp'] = np.nan
if self.regress_flag:
y_true, y_pred = spec_df.rg_targets.tolist(), spec_df.regressions.tolist()
stats_dict["rg_RMSE"] = RMSE(y_true, y_pred)
stats_dict["rg_MAE"] = MAE(y_true, y_pred)
stats_dict["rg_bin_accuracy"] = accuracy(spec_df.rg_bin_target.tolist(), spec_df.rg_bins.tolist())
stats_dict["rg_uncertainty"] = spec_df.rg_uncertainties.mean()
if hasattr(self, "seg_df") and not boxes_only:
seg_cl_df = seg_cl_df.groupby(['pid'], as_index=False).agg(
{dice_col: "mean", "fold": "first"}) # mean of all rois per patient in this epoch
stats_dict["dice"] = seg_cl_df.loc[:,dice_col].mean() #mean of all patients
stats_dict["dice_std"] = seg_cl_df.loc[:, dice_col].std()
# for the aggregated test set case, additionally get the scores for averaging over fold results.
if self.cf.evaluate_fold_means and len(df.fold.unique()) > 1 and self.mode in ["test", "analysis"]:
aucs = []
aps = []
for fold in df.fold.unique():
fold_df = spec_df[spec_df.fold == fold]
if (fold_df.class_label==1).any():
aps.append(
average_precision_score(fold_df.class_label.tolist(), fold_df.pred_score.tolist()))
if len(fold_df.class_label.unique())>1:
aucs.append(roc_auc_score(fold_df.class_label.tolist(), fold_df.pred_score.tolist()))
stats_dict['auc_folds_mean'] = np.mean(aucs)
stats_dict['auc_folds_std'] = np.std(aucs)
stats_dict['ap_folds_mean'] = np.mean(aps)
stats_dict['ap_folds_std'] = np.std(aps)
if hasattr(self, "seg_df") and not boxes_only and self.cf.evaluate_fold_means and len(seg_cl_df.fold.unique()) > 1:
fold_means = seg_cl_df.groupby(['fold'], as_index=True).agg({dice_col:"mean"})
stats_dict["dice_folds_mean"] = float(fold_means.mean())
stats_dict["dice_folds_std"] = float(fold_means.std())
all_stats.append(stats_dict)
# -------------- monitoring, visualisation -----------------
# fill new results into monitor_metrics dict. for simplicity, only one class (of interest) is monitored on patient level.
patient_interests = [self.cf.class_dict[self.cf.patient_class_of_interest],]
if hasattr(self.cf, "bin_dict"):
patient_interests += [self.cf.bin_dict[self.cf.patient_bin_of_interest]]
if monitor_metrics is not None and (score_level != 'patient' or cl_name in patient_interests):
name = 'patient_'+cl_name if score_level == 'patient' else cl_name
for metric in self.cf.metrics:
if metric in stats_dict.keys():
monitor_metrics[name + '_'+metric].append(stats_dict[metric])
else:
print("WARNING: skipped monitor metric {}_{} since not avail".format(name, metric))
# histograms
if self.cf.plot_prediction_histograms:
out_filename = os.path.join(self.hist_dir, 'pred_hist_{}_{}_{}_{}'.format(
self.cf.fold, self.mode, score_level, cl_name))
plg.plot_prediction_hist(self.cf, spec_df, out_filename)
# analysis of the hyper-parameter cf.min_det_thresh, for optimization on validation set.
if self.cf.scan_det_thresh and "val" in self.mode:
conf_threshs = list(np.arange(0.8, 1, 0.02))
pool = Pool(processes=self.cf.n_workers)
mp_inputs = [[spec_df, ii, self.cf.per_patient_ap] for ii in conf_threshs]
aps = pool.map(get_roi_ap_from_df, mp_inputs, chunksize=1)
pool.close()
pool.join()
self.logger.info('results from scanning over det_threshs: {}'.format([[i, j] for i, j in zip(conf_threshs, aps)]))
class_means = pd.DataFrame(columns=self.cf.report_score_level)
for slevel in self.cf.report_score_level:
level_stats = pd.DataFrame([stats for stats in all_stats if slevel in stats["name"]])[self.cf.metrics]
class_means.loc[:, slevel] = level_stats.mean()
all_stats.extend([{"name": 'fold_{} {} {}'.format(self.cf.fold, slevel, "class_means"), **level_means} for
slevel, level_means in class_means.to_dict().items()])
if self.cf.plot_stat_curves:
out_filename = os.path.join(self.curves_dir, '{}_{}_stat_curves'.format(self.cf.fold, self.mode))
plg.plot_stat_curves(self.cf, all_stats, out_filename)
if self.cf.plot_prediction_histograms and hasattr(df, "cluster_n_missing") and df.cluster_n_missing.notna().any():
out_filename = os.path.join(self.hist_dir, 'n_missing_hist_{}_{}.png'.format(self.cf.fold, self.mode))
plg.plot_wbc_n_missing(self.cf, df, outfile=out_filename)
return all_stats, monitor_metrics
def score_test_df(self, max_fold=None, internal_df=True):
"""
Writes out resulting scores to text files: First checks for class-internal-df (typically current) fold,
gets resulting scores, writes them to a text file and pickles data frame. Also checks if data-frame pickles of
all folds of cross-validation exist in exp_dir. If true, loads all dataframes, aggregates test sets over folds,
and calculates and writes out overall metrics.
"""
# this should maybe be extended to auc, ap stds.
metrics_to_score = self.cf.metrics.copy() # + [ m+ext for m in self.cf.metrics if "dice" in m for ext in ["_std"]]
if internal_df:
self.test_df.to_pickle(os.path.join(self.cf.test_dir, '{}_test_df.pkl'.format(self.cf.fold)))
if hasattr(self, "seg_df"):
self.seg_df.to_pickle(os.path.join(self.cf.test_dir, '{}_test_seg_df.pkl'.format(self.cf.fold)))
stats, _ = self.return_metrics(self.test_df, self.cf.class_dict)
with open(os.path.join(self.cf.test_dir, 'results.txt'), 'a') as handle:
handle.write('\n****************************\n')
handle.write('\nresults for fold {}, {} \n'.format(self.cf.fold, time.strftime("%d/%m/%y %H:%M:%S")))
handle.write('\n****************************\n')
handle.write('\nfold df shape {}\n \n'.format(self.test_df.shape))
for s in stats:
for metric in metrics_to_score:
if metric in s.keys(): #needed as long as no dice on patient level poss
if "accuracy" in metric:
handle.write('{} {:0.4f} '.format(metric, s[metric]))
else:
handle.write('{} {:0.3f} '.format(metric, s[metric]))
else:
print("WARNING: skipped metric {} since not avail".format(metric))
handle.write('{} \n'.format(s['name']))
- fold_df_paths = sorted([ii for ii in os.listdir(self.cf.test_dir) if 'test_df.pkl' in ii])
- fold_seg_df_paths = sorted([ii for ii in os.listdir(self.cf.test_dir) if 'test_seg_df.pkl' in ii])
- for paths in [fold_df_paths, fold_seg_df_paths]:
- assert len(paths)<= self.cf.n_cv_splits, "found {} > nr of cv splits results dfs in {}".format(len(paths), self.cf.test_dir)
+
if max_fold is None:
max_fold = self.cf.n_cv_splits-1
if self.cf.fold == max_fold:
print("max fold/overall stats triggered")
+ self.cf.fold = 'overall'
if self.cf.evaluate_fold_means:
metrics_to_score += [m + ext for m in self.cf.metrics for ext in ("_folds_mean", "_folds_std")]
- with open(os.path.join(self.cf.test_dir, 'results.txt'), 'a') as handle:
+ if not self.cf.held_out_test_set or self.cf.eval_test_fold_wise:
+ fold_df_paths = sorted([ii for ii in os.listdir(self.cf.test_dir) if 'test_df.pkl' in ii])
+ fold_seg_df_paths = sorted([ii for ii in os.listdir(self.cf.test_dir) if 'test_seg_df.pkl' in ii])
+ for paths in [fold_df_paths, fold_seg_df_paths]:
+ assert len(paths) <= self.cf.n_cv_splits, "found {} > nr of cv splits results dfs in {}".format(
+ len(paths), self.cf.test_dir)
+ with open(os.path.join(self.cf.test_dir, 'results.txt'), 'a') as handle:
- self.cf.fold = 'overall'
- dfs_list = [pd.read_pickle(os.path.join(self.cf.test_dir, ii)) for ii in fold_df_paths]
- seg_dfs_list = [pd.read_pickle(os.path.join(self.cf.test_dir, ii)) for ii in fold_seg_df_paths]
- self.test_df = pd.concat(dfs_list, sort=True)
- if len(seg_dfs_list)>0:
- self.seg_df = pd.concat(seg_dfs_list, sort=True)
- stats, _ = self.return_metrics(self.test_df, self.cf.class_dict)
+ dfs_list = [pd.read_pickle(os.path.join(self.cf.test_dir, ii)) for ii in fold_df_paths]
+ seg_dfs_list = [pd.read_pickle(os.path.join(self.cf.test_dir, ii)) for ii in fold_seg_df_paths]
- handle.write('\n****************************\n')
- handle.write('\nOVERALL RESULTS \n')
- handle.write('\n****************************\n')
- handle.write('\ndf shape \n \n'.format(self.test_df.shape))
- for s in stats:
- for metric in metrics_to_score:
- if metric in s.keys():
- handle.write('{} {:0.3f} '.format(metric, s[metric]))
- handle.write('{} \n'.format(s['name']))
+ self.test_df = pd.concat(dfs_list, sort=True)
+ if len(seg_dfs_list)>0:
+ self.seg_df = pd.concat(seg_dfs_list, sort=True)
+ stats, _ = self.return_metrics(self.test_df, self.cf.class_dict)
+
+ handle.write('\n****************************\n')
+ handle.write('\nOVERALL RESULTS \n')
+ handle.write('\n****************************\n')
+ handle.write('\ndf shape \n \n'.format(self.test_df.shape))
+ for s in stats:
+ for metric in metrics_to_score:
+ if metric in s.keys():
+ handle.write('{} {:0.3f} '.format(metric, s[metric]))
+ handle.write('{} \n'.format(s['name']))
results_table_path = os.path.join(self.cf.test_dir,"../../", 'results_table.csv')
with open(results_table_path, 'a') as handle:
#---column headers---
handle.write('\n{},'.format("Experiment Name"))
handle.write('{},'.format("Time Stamp"))
handle.write('{},'.format("Samples Seen"))
handle.write('{},'.format("Spatial Dim"))
handle.write('{},'.format("Patch Size"))
handle.write('{},'.format("CV Folds"))
handle.write('{},'.format("{}-clustering IoU".format(self.cf.clustering)))
handle.write('{},'.format("Merge-2D-to-3D IoU"))
if hasattr(self.cf, "test_against_exact_gt"):
handle.write('{},'.format('Exact GT'))
for s in stats:
- assert "overall" in s['name'].split(" ")[0]
if self.cf.class_dict[self.cf.patient_class_of_interest] in s['name'] or "mean" in s["name"]:
for metric in metrics_to_score:
if metric in s.keys() and not np.isnan(s[metric]):
if metric=='ap':
handle.write('{}_{} : {}_{},'.format(*s['name'].split(" ")[1:], metric, int(np.mean(self.cf.ap_match_ious)*100)))
elif not "folds_std" in metric:
handle.write('{}_{} : {},'.format(*s['name'].split(" ")[1:], metric))
else:
print("WARNING: skipped metric {} since not avail".format(metric))
handle.write('\n')
#--- columns content---
handle.write('{},'.format(self.cf.exp_dir.split(os.sep)[-1]))
handle.write('{},'.format(time.strftime("%d%b%y %H:%M:%S")))
handle.write('{},'.format(self.cf.num_epochs*self.cf.num_train_batches*self.cf.batch_size))
handle.write('{}D,'.format(self.cf.dim))
handle.write('{},'.format("x".join([str(self.cf.patch_size[i]) for i in range(self.cf.dim)])))
handle.write('{},'.format(str(self.test_df.fold.unique().tolist()).replace(",", "")))
handle.write('{},'.format(self.cf.clustering_iou if self.cf.clustering else str("N/A")))
handle.write('{},'.format(self.cf.merge_3D_iou if self.cf.merge_2D_to_3D_preds else str("N/A")))
if hasattr(self.cf, "test_against_exact_gt"):
handle.write('{},'.format(self.cf.test_against_exact_gt))
for s in stats:
if self.cf.class_dict[self.cf.patient_class_of_interest] in s['name'] or "mean" in s["name"]:
for metric in metrics_to_score:
if metric in s.keys() and not np.isnan(s[metric]): # needed as long as no dice on patient level possible
if "folds_mean" in metric:
handle.write('{:0.3f}\u00B1{:0.3f}, '.format(s[metric], s["_".join((*metric.split("_")[:-1], "std"))]))
elif not "folds_std" in metric:
handle.write('{:0.3f}, '.format(s[metric]))
handle.write('\n')
with open(os.path.join(self.cf.test_dir, 'results_extr_scores.txt'), 'w') as handle:
handle.write('\n****************************\n')
handle.write('\nextremal scores for fold {} \n'.format(self.cf.fold))
handle.write('\n****************************\n')
# want: pid & fold (&other) of highest scoring tp & fp in test_df
for cl in self.cf.class_dict.keys():
print("\nClass {}".format(self.cf.class_dict[cl]), file=handle)
cl_df = self.test_df[self.test_df.pred_class == cl] #.dropna(axis=1)
for det_type in ['det_tp', 'det_fp']:
filtered_df = cl_df[cl_df.det_type==det_type]
print("\nHighest scoring {} of class {}".format(det_type, self.cf.class_dict[cl]), file=handle)
if len(filtered_df)>0:
print(filtered_df.loc[filtered_df.pred_score.idxmax()], file=handle)
else:
print("No detections of type {} for class {} in this df".format(det_type, self.cf.class_dict[cl]), file=handle)
handle.write('\n****************************\n')
diff --git a/exec.py b/exec.py
index 85e7716..d6f00e2 100644
--- a/exec.py
+++ b/exec.py
@@ -1,344 +1,344 @@
#!/usr/bin/env python
# Copyright 2019 Division of Medical Image Computing, German Cancer Research Center (DKFZ).
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
# ==============================================================================
""" execution script. this where all routines come together and the only script you need to call.
refer to parse args below to see options for execution.
"""
import plotting as plg
import os
import warnings
import argparse
import time
import torch
import utils.exp_utils as utils
from evaluator import Evaluator
from predictor import Predictor
for msg in ["Attempting to set identical bottom==top results",
"This figure includes Axes that are not compatible with tight_layout",
"Data has no positive values, and therefore cannot be log-scaled.",
".*invalid value encountered in true_divide.*"]:
warnings.filterwarnings("ignore", msg)
def train(cf, logger):
"""
performs the training routine for a given fold. saves plots and selected parameters to the experiment dir
specified in the configs. logs to file and tensorboard.
"""
logger.info('performing training in {}D over fold {} on experiment {} with model {}'.format(
cf.dim, cf.fold, cf.exp_dir, cf.model))
logger.time("train_val")
# -------------- inits and settings -----------------
net = model.net(cf, logger).cuda()
if cf.optimizer == "ADAMW":
optimizer = torch.optim.AdamW(utils.parse_params_for_optim(net, weight_decay=cf.weight_decay,
exclude_from_wd=cf.exclude_from_wd),
lr=cf.learning_rate[0])
elif cf.optimizer == "SGD":
optimizer = torch.optim.SGD(utils.parse_params_for_optim(net, weight_decay=cf.weight_decay),
lr=cf.learning_rate[0], momentum=0.3)
if cf.dynamic_lr_scheduling:
scheduler = torch.optim.lr_scheduler.ReduceLROnPlateau(optimizer, mode=cf.scheduling_mode, factor=cf.lr_decay_factor,
patience=cf.scheduling_patience)
model_selector = utils.ModelSelector(cf, logger)
starting_epoch = 1
if cf.resume:
checkpoint_path = os.path.join(cf.fold_dir, "last_state.pth")
starting_epoch, net, optimizer, model_selector = \
utils.load_checkpoint(checkpoint_path, net, optimizer, model_selector)
logger.info('resumed from checkpoint {} to epoch {}'.format(checkpoint_path, starting_epoch))
# prepare monitoring
monitor_metrics = utils.prepare_monitoring(cf)
logger.info('loading dataset and initializing batch generators...')
batch_gen = data_loader.get_train_generators(cf, logger)
# -------------- training -----------------
for epoch in range(starting_epoch, cf.num_epochs + 1):
logger.info('starting training epoch {}/{}'.format(epoch, cf.num_epochs))
logger.time("train_epoch")
net.train()
train_results_list = []
train_evaluator = Evaluator(cf, logger, mode='train')
for i in range(cf.num_train_batches):
logger.time("train_batch_loadfw")
batch = next(batch_gen['train'])
batch_gen['train'].generator.stats['roi_counts'] += batch['roi_counts']
batch_gen['train'].generator.stats['empty_counts'] += batch['empty_counts']
logger.time("train_batch_loadfw")
logger.time("train_batch_netfw")
results_dict = net.train_forward(batch)
logger.time("train_batch_netfw")
logger.time("train_batch_bw")
optimizer.zero_grad()
results_dict['torch_loss'].backward()
if cf.clip_norm:
torch.nn.utils.clip_grad_norm_(net.parameters(), cf.clip_norm, norm_type=2) # gradient clipping
optimizer.step()
train_results_list.append(({k:v for k,v in results_dict.items() if k != "seg_preds"}, batch["pid"])) # slim res dict
if not cf.server_env:
print("\rFinished training batch " +
"{}/{} in {:.1f}s ({:.2f}/{:.2f} forw load/net, {:.2f} backw).".format(i+1, cf.num_train_batches,
logger.get_time("train_batch_loadfw")+
logger.get_time("train_batch_netfw")
+logger.time("train_batch_bw"),
logger.get_time("train_batch_loadfw",reset=True),
logger.get_time("train_batch_netfw", reset=True),
logger.get_time("train_batch_bw", reset=True)), end="", flush=True)
print()
#--------------- train eval ----------------
if (epoch-1)%cf.plot_frequency==0:
# view an example batch
utils.split_off_process(plg.view_batch, cf, batch, results_dict, has_colorchannels=cf.has_colorchannels,
show_gt_labels=True, get_time="train-example plot",
out_file=os.path.join(cf.plot_dir, 'batch_example_train_{}.png'.format(cf.fold)))
logger.time("evals")
_, monitor_metrics['train'] = train_evaluator.evaluate_predictions(train_results_list, monitor_metrics['train'])
logger.time("evals")
logger.time("train_epoch", toggle=False)
del train_results_list
#----------- validation ------------
logger.info('starting validation in mode {}.'.format(cf.val_mode))
logger.time("val_epoch")
with torch.no_grad():
net.eval()
val_results_list = []
val_evaluator = Evaluator(cf, logger, mode=cf.val_mode)
val_predictor = Predictor(cf, net, logger, mode='val')
for i in range(batch_gen['n_val']):
logger.time("val_batch")
batch = next(batch_gen[cf.val_mode])
if cf.val_mode == 'val_patient':
results_dict = val_predictor.predict_patient(batch)
elif cf.val_mode == 'val_sampling':
results_dict = net.train_forward(batch, is_validation=True)
val_results_list.append([results_dict, batch["pid"]])
if not cf.server_env:
print("\rFinished validation {} {}/{} in {:.1f}s.".format('patient' if cf.val_mode=='val_patient' else 'batch',
i + 1, batch_gen['n_val'],
logger.time("val_batch")), end="", flush=True)
print()
#------------ val eval -------------
if (epoch - 1) % cf.plot_frequency == 0:
utils.split_off_process(plg.view_batch, cf, batch, results_dict, has_colorchannels=cf.has_colorchannels,
show_gt_labels=True, get_time="val-example plot",
out_file=os.path.join(cf.plot_dir, 'batch_example_val_{}.png'.format(cf.fold)))
logger.time("evals")
_, monitor_metrics['val'] = val_evaluator.evaluate_predictions(val_results_list, monitor_metrics['val'])
model_selector.run_model_selection(net, optimizer, monitor_metrics, epoch)
del val_results_list
#----------- monitoring -------------
monitor_metrics.update({"lr":
{str(g) : group['lr'] for (g, group) in enumerate(optimizer.param_groups)}})
logger.metrics2tboard(monitor_metrics, global_step=epoch)
logger.time("evals")
logger.info('finished epoch {}/{}, took {:.2f}s. train total: {:.2f}s, average: {:.2f}s. val total: {:.2f}s, average: {:.2f}s.'.format(
epoch, cf.num_epochs, logger.get_time("train_epoch")+logger.time("val_epoch"), logger.get_time("train_epoch"),
logger.get_time("train_epoch", reset=True)/cf.num_train_batches, logger.get_time("val_epoch"),
logger.get_time("val_epoch", reset=True)/batch_gen["n_val"]))
logger.info("time for evals: {:.2f}s".format(logger.get_time("evals", reset=True)))
#-------------- scheduling -----------------
if cf.dynamic_lr_scheduling:
scheduler.step(monitor_metrics["val"][cf.scheduling_criterion][-1])
else:
for param_group in optimizer.param_groups:
param_group['lr'] = cf.learning_rate[epoch-1]
logger.time("train_val")
logger.info("Training and validating over {} epochs took {}".format(cf.num_epochs, logger.get_time("train_val", format="hms", reset=True)))
batch_gen['train'].generator.print_stats(logger, plot=True)
def test(cf, logger, max_fold=None):
"""performs testing for a given fold (or held out set). saves stats in evaluator.
"""
logger.time("test_fold")
logger.info('starting testing model of fold {} in exp {}'.format(cf.fold, cf.exp_dir))
net = model.net(cf, logger).cuda()
batch_gen = data_loader.get_test_generator(cf, logger)
test_predictor = Predictor(cf, net, logger, mode='test')
test_results_list = test_predictor.predict_test_set(batch_gen, return_results = not hasattr(
cf, "eval_test_separately") or not cf.eval_test_separately)
if test_results_list is not None:
test_evaluator = Evaluator(cf, logger, mode='test')
test_evaluator.evaluate_predictions(test_results_list)
test_evaluator.score_test_df(max_fold=max_fold)
logger.info('Testing of fold {} took {}.\n'.format(cf.fold, logger.get_time("test_fold", reset=True, format="hms")))
if __name__ == '__main__':
stime = time.time()
parser = argparse.ArgumentParser()
parser.add_argument('--dataset_name', type=str, default='toy',
help="path to the dataset-specific code in source_dir/datasets")
parser.add_argument('--exp_dir', type=str, default='/home/gregor/Documents/regrcnn/datasets/toy/experiments/dev',
help='path to experiment dir. will be created if non existent.')
parser.add_argument('-m', '--mode', type=str, default='train_test', help='one out of: create_exp, analysis, train, train_test, or test')
parser.add_argument('-f', '--folds', nargs='+', type=int, default=None, help='None runs over all folds in CV. otherwise specify list of folds.')
parser.add_argument('--server_env', default=False, action='store_true', help='change IO settings to deploy models on a cluster.')
parser.add_argument('--data_dest', type=str, default=None, help="path to final data folder if different from config")
parser.add_argument('--use_stored_settings', default=False, action='store_true',
help='load configs from existing exp_dir instead of source dir. always done for testing, '
'but can be set to true to do the same for training. useful in job scheduler environment, '
'where source code might change before the job actually runs.')
parser.add_argument('--resume', action="store_true", default=False,
help='if given, resume from checkpoint(s) of the specified folds.')
parser.add_argument('-d', '--dev', default=False, action='store_true', help="development mode: shorten everything")
args = parser.parse_args()
args.dataset_name = os.path.join("datasets", args.dataset_name) if not "datasets" in args.dataset_name else args.dataset_name
folds = args.folds
resume = None if args.resume in ['None', 'none'] else args.resume
if args.mode == 'create_exp':
cf = utils.prep_exp(args.dataset_name, args.exp_dir, args.server_env, use_stored_settings=False)
logger = utils.get_logger(cf.exp_dir, cf.server_env, -1)
logger.info('created experiment directory at {}'.format(args.exp_dir))
elif args.mode == 'train' or args.mode == 'train_test':
cf = utils.prep_exp(args.dataset_name, args.exp_dir, args.server_env, args.use_stored_settings)
if args.dev:
folds = [0,1]
cf.batch_size, cf.num_epochs, cf.min_save_thresh, cf.save_n_models = 3 if cf.dim==2 else 1, 2, 0, 2
cf.num_train_batches, cf.num_val_batches, cf.max_val_patients = 5, 1, 1
- cf.test_n_epochs = cf.save_n_models
- cf.max_test_patients = 1
+ cf.test_n_epochs, cf.max_test_patients = cf.save_n_models, 2
torch.backends.cudnn.benchmark = cf.dim==3
else:
torch.backends.cudnn.benchmark = cf.cuda_benchmark
if args.data_dest is not None:
cf.data_dest = args.data_dest
logger = utils.get_logger(cf.exp_dir, cf.server_env, cf.sysmetrics_interval)
data_loader = utils.import_module('data_loader', os.path.join(args.dataset_name, 'data_loader.py'))
model = utils.import_module('model', cf.model_path)
logger.info("loaded model from {}".format(cf.model_path))
if folds is None:
folds = range(cf.n_cv_splits)
for fold in folds:
"""k-fold cross-validation: the dataset is split into k equally-sized folds, one used for validation,
one for testing, the rest for training. This loop iterates k-times over the dataset, cyclically moving the
splits. k==folds, fold in [0,folds) says which split is used for testing.
"""
cf.fold_dir = os.path.join(cf.exp_dir, 'fold_{}'.format(fold)); cf.fold = fold
logger.set_logfile(fold=fold)
cf.resume = resume
if not os.path.exists(cf.fold_dir):
os.mkdir(cf.fold_dir)
train(cf, logger)
cf.resume = None
if args.mode == 'train_test':
test(cf, logger)
elif args.mode == 'test':
cf = utils.prep_exp(args.dataset_name, args.exp_dir, args.server_env, use_stored_settings=True, is_training=False)
if args.data_dest is not None:
cf.data_dest = args.data_dest
logger = utils.get_logger(cf.exp_dir, cf.server_env, cf.sysmetrics_interval)
data_loader = utils.import_module('data_loader', os.path.join(args.dataset_name, 'data_loader.py'))
model = utils.import_module('model', cf.model_path)
logger.info("loaded model from {}".format(cf.model_path))
fold_dirs = sorted([os.path.join(cf.exp_dir, f) for f in os.listdir(cf.exp_dir) if
os.path.isdir(os.path.join(cf.exp_dir, f)) and f.startswith("fold")])
if folds is None:
folds = range(cf.n_cv_splits)
if args.dev:
folds = folds[:2]
- cf.batch_size, cf.max_test_patients, cf.test_n_epochs = 1 if cf.dim==2 else 1, 2, 2
+ cf.max_test_patients, cf.test_n_epochs = 2, 2
else:
torch.backends.cudnn.benchmark = cf.cuda_benchmark
for fold in folds:
cf.fold_dir = os.path.join(cf.exp_dir, 'fold_{}'.format(fold)); cf.fold = fold
logger.set_logfile(fold=fold)
if cf.fold_dir in fold_dirs:
test(cf, logger, max_fold=max([int(f[-1]) for f in fold_dirs]))
else:
logger.info("Skipping fold {} since no model parameters found.".format(fold))
# load raw predictions saved by predictor during testing, run aggregation algorithms and evaluation.
elif args.mode == 'analysis':
""" analyse already saved predictions.
"""
cf = utils.prep_exp(args.dataset_name, args.exp_dir, args.server_env, use_stored_settings=True, is_training=False)
logger = utils.get_logger(cf.exp_dir, cf.server_env, cf.sysmetrics_interval)
if cf.held_out_test_set and not cf.eval_test_fold_wise:
predictor = Predictor(cf, net=None, logger=logger, mode='analysis')
results_list = predictor.load_saved_predictions()
logger.info('starting evaluation...')
- cf.fold = 0
+ cf.fold = "overall"
evaluator = Evaluator(cf, logger, mode='test')
evaluator.evaluate_predictions(results_list)
- evaluator.score_test_df(max_fold=0)
+ evaluator.score_test_df(max_fold=cf.fold)
else:
fold_dirs = sorted([os.path.join(cf.exp_dir, f) for f in os.listdir(cf.exp_dir) if
os.path.isdir(os.path.join(cf.exp_dir, f)) and f.startswith("fold")])
if args.dev:
+ cf.test_n_epochs = 2
fold_dirs = fold_dirs[:1]
if folds is None:
folds = range(cf.n_cv_splits)
for fold in folds:
cf.fold = fold; cf.fold_dir = os.path.join(cf.exp_dir, 'fold_{}'.format(cf.fold))
logger.set_logfile(fold=fold)
if cf.fold_dir in fold_dirs:
predictor = Predictor(cf, net=None, logger=logger, mode='analysis')
results_list = predictor.load_saved_predictions()
# results_list[x][1] is pid, results_list[x][0] is list of len samples-per-patient, each entry hlds
# list of boxes per that sample, i.e., len(results_list[x][y][0]) would be nr of boxes in sample y of patient x
logger.info('starting evaluation...')
evaluator = Evaluator(cf, logger, mode='test')
evaluator.evaluate_predictions(results_list)
max_fold = max([int(f[-1]) for f in fold_dirs])
evaluator.score_test_df(max_fold=max_fold)
else:
logger.info("Skipping fold {} since no model parameters found.".format(fold))
else:
raise ValueError('mode "{}" specified in args is not implemented.'.format(args.mode))
mins, secs = divmod((time.time() - stime), 60)
h, mins = divmod(mins, 60)
t = "{:d}h:{:02d}m:{:02d}s".format(int(h), int(mins), int(secs))
logger.info("{} total runtime: {}".format(os.path.split(__file__)[1], t))
del logger
torch.cuda.empty_cache()
diff --git a/utils/dataloader_utils.py b/utils/dataloader_utils.py
index 45ba333..5428f79 100644
--- a/utils/dataloader_utils.py
+++ b/utils/dataloader_utils.py
@@ -1,727 +1,730 @@
#!/usr/bin/env python
# Copyright 2019 Division of Medical Image Computing, German Cancer Research Center (DKFZ).
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
# ==============================================================================
import plotting as plg
import os
from multiprocessing import Pool, Lock
import pickle
import warnings
import numpy as np
import pandas as pd
from batchgenerators.transforms.abstract_transforms import AbstractTransform
from scipy.ndimage.measurements import label as lb
from torch.utils.data import Dataset as torchDataset
from batchgenerators.dataloading.data_loader import SlimDataLoaderBase
import utils.exp_utils as utils
import data_manager as dmanager
for msg in ["This figure includes Axes that are not compatible with tight_layout",
"Data has no positive values, and therefore cannot be log-scaled."]:
warnings.filterwarnings("ignore", msg)
class AttributeDict(dict):
__getattr__ = dict.__getitem__
__setattr__ = dict.__setitem__
##################################
# data loading, organisation #
##################################
class fold_generator:
"""
generates splits of indices for a given length of a dataset to perform n-fold cross-validation.
splits each fold into 3 subsets for training, validation and testing.
This form of cross validation uses an inner loop test set, which is useful if test scores shall be reported on a
statistically reliable amount of patients, despite limited size of a dataset.
If hold out test set is provided and hence no inner loop test set needed, just add test_idxs to the training data in the dataloader.
This creates straight-forward train-val splits.
:returns names list: list of len n_splits. each element is a list of len 3 for train_ix, val_ix, test_ix.
"""
def __init__(self, seed, n_splits, len_data):
"""
:param seed: Random seed for splits.
:param n_splits: number of splits, e.g. 5 splits for 5-fold cross-validation
:param len_data: number of elements in the dataset.
"""
self.tr_ix = []
self.val_ix = []
self.te_ix = []
self.slicer = None
self.missing = 0
self.fold = 0
self.len_data = len_data
self.n_splits = n_splits
self.myseed = seed
self.boost_val = 0
def init_indices(self):
t = list(np.arange(self.l))
# round up to next splittable data amount.
split_length = int(np.ceil(len(t) / float(self.n_splits)))
self.slicer = split_length
self.mod = len(t) % self.n_splits
if self.mod > 0:
# missing is the number of folds, in which the new splits are reduced to account for missing data.
self.missing = self.n_splits - self.mod
self.te_ix = t[:self.slicer]
self.tr_ix = t[self.slicer:]
self.val_ix = self.tr_ix[:self.slicer]
self.tr_ix = self.tr_ix[self.slicer:]
def new_fold(self):
slicer = self.slicer
if self.fold < self.missing :
slicer = self.slicer - 1
temp = self.te_ix
# catch exception mod == 1: test set collects 1+ data since walk through both roudned up splits.
# account for by reducing last fold split by 1.
if self.fold == self.n_splits-2 and self.mod ==1:
temp += self.val_ix[-1:]
self.val_ix = self.val_ix[:-1]
self.te_ix = self.val_ix
self.val_ix = self.tr_ix[:slicer]
self.tr_ix = self.tr_ix[slicer:] + temp
def get_fold_names(self):
names_list = []
rgen = np.random.RandomState(self.myseed)
cv_names = np.arange(self.len_data)
rgen.shuffle(cv_names)
self.l = len(cv_names)
self.init_indices()
for split in range(self.n_splits):
train_names, val_names, test_names = cv_names[self.tr_ix], cv_names[self.val_ix], cv_names[self.te_ix]
names_list.append([train_names, val_names, test_names, self.fold])
self.new_fold()
self.fold += 1
return names_list
class FoldGenerator():
r"""takes a set of elements (identifiers) and randomly splits them into the specified amt of subsets.
"""
def __init__(self, identifiers, seed, n_splits=5):
self.ids = np.array(identifiers)
self.n_splits = n_splits
self.seed = seed
def generate_splits(self, n_splits=None):
if n_splits is None:
n_splits = self.n_splits
rgen = np.random.RandomState(self.seed)
rgen.shuffle(self.ids)
self.splits = list(np.array_split(self.ids, n_splits, axis=0)) # already returns list, but to be sure
return self.splits
class Dataset(torchDataset):
r"""Parent Class for actual Dataset classes to inherit from!
"""
def __init__(self, cf, data_sourcedir=None):
super(Dataset, self).__init__()
self.cf = cf
self.data_sourcedir = cf.data_sourcedir if data_sourcedir is None else data_sourcedir
self.data_dir = cf.data_dir if hasattr(cf, 'data_dir') else self.data_sourcedir
self.data_dest = cf.data_dest if hasattr(cf, "data_dest") else self.data_sourcedir
self.data = {}
self.set_ids = []
def copy_data(self, cf, file_subset, keep_packed=False, del_after_unpack=False):
if os.path.normpath(self.data_sourcedir) != os.path.normpath(self.data_dest):
self.data_sourcedir = os.path.join(self.data_sourcedir, '')
args = AttributeDict({
"source" : self.data_sourcedir,
"destination" : self.data_dest,
"recursive" : True,
"cp_only_npz" : False,
"keep_packed" : keep_packed,
"del_after_unpack" : del_after_unpack,
"threads" : 16 if self.cf.server_env else os.cpu_count()
})
dmanager.copy(args, file_subset=file_subset)
self.data_dir = self.data_dest
def __len__(self):
return len(self.data)
def __getitem__(self, id):
"""Return a sample of the dataset, i.e.,the dict of the id
"""
return self.data[id]
def __iter__(self):
return self.data.__iter__()
def init_FoldGenerator(self, seed, n_splits):
self.fg = FoldGenerator(self.set_ids, seed=seed, n_splits=n_splits)
def generate_splits(self, check_file):
if not os.path.exists(check_file):
self.fg.generate_splits()
with open(check_file, 'wb') as handle:
pickle.dump(self.fg.splits, handle)
else:
with open(check_file, 'rb') as handle:
self.fg.splits = pickle.load(handle)
def calc_statistics(self, subsets=None, plot_dir=None, overall_stats=True):
if self.df is None:
self.df = pd.DataFrame()
balance_t = self.cf.balance_target if hasattr(self.cf, "balance_target") else "class_targets"
self.df._metadata.append(balance_t)
if balance_t=="class_targets":
mapper = lambda cl_id: self.cf.class_id2label[cl_id]
labels = self.cf.class_id2label.values()
elif balance_t=="rg_bin_targets":
mapper = lambda rg_bin: self.cf.bin_id2label[rg_bin]
labels = self.cf.bin_id2label.values()
# elif balance_t=="regression_targets":
# # todo this wont work
# mapper = lambda rg_val: AttributeDict({"name":rg_val}) #self.cf.bin_id2label[self.cf.rg_val_to_bin_id(rg_val)]
# labels = self.cf.bin_id2label.values()
elif balance_t=="lesion_gleasons":
mapper = lambda gs: self.cf.gs2label[gs]
labels = self.cf.gs2label.values()
else:
mapper = lambda x: AttributeDict({"name":x})
labels = None
for pid, subj_data in self.data.items():
unique_ts, counts = np.unique(subj_data[balance_t], return_counts=True)
self.df = self.df.append(pd.DataFrame({"pid": [pid],
**{mapper(unique_ts[i]).name: [counts[i]] for i in
range(len(unique_ts))}}), ignore_index=True, sort=True)
self.df = self.df.fillna(0)
if overall_stats:
df = self.df.drop("pid", axis=1)
df = df.reindex(sorted(df.columns), axis=1).astype('uint32')
print("Overall dataset roi counts per target kind:"); print(df.sum())
if subsets is not None:
self.df["subset"] = np.nan
self.df["display_order"] = np.nan
for ix, (subset, pids) in enumerate(subsets.items()):
self.df.loc[self.df.pid.isin(pids), "subset"] = subset
self.df.loc[self.df.pid.isin(pids), "display_order"] = ix
df = self.df.groupby("subset").agg("sum").drop("pid", axis=1, errors='ignore').astype('int64')
df = df.sort_values(by=['display_order']).drop('display_order', axis=1)
df = df.reindex(sorted(df.columns), axis=1)
print("Fold {} dataset roi counts per target kind:".format(self.cf.fold)); print(df)
if plot_dir is not None:
os.makedirs(plot_dir, exist_ok=True)
if subsets is not None:
plg.plot_fold_stats(self.cf, df, labels, os.path.join(plot_dir, "data_stats_fold_" + str(self.cf.fold))+".pdf")
if overall_stats:
plg.plot_data_stats(self.cf, df, labels, os.path.join(plot_dir, 'data_stats_overall.pdf'))
return df, labels
def get_class_balanced_patients(all_pids, class_targets, batch_size, num_classes, random_ratio=0):
'''
samples towards equilibrium of classes (on basis of total RoI counts). for highly imbalanced dataset, this might be a too strong requirement.
:param class_targets: dic holding {patient_specifier : ROI class targets}, list position of ROI target corresponds to respective seg label - 1
:param batch_size:
:param num_classes:
:return:
'''
# assert len(all_pids)>=batch_size, "not enough eligible pids {} to form a single batch of size {}".format(len(all_pids), batch_size)
class_counts = {k: 0 for k in range(1,num_classes+1)}
not_picked = np.array(all_pids)
batch_patients = np.empty((batch_size,), dtype=not_picked.dtype)
rarest_class = np.random.randint(1,num_classes+1)
for ix in range(batch_size):
if len(not_picked) == 0:
warnings.warn("Dataset too small to generate batch with unique samples; => recycling.")
not_picked = np.array(all_pids)
np.random.shuffle(not_picked) #this could actually go outside(above) the loop.
pick = not_picked[0]
for cand in not_picked:
if np.count_nonzero(class_targets[cand] == rarest_class) > 0:
pick = cand
cand_rarest_class = np.argmin([np.count_nonzero(class_targets[cand] == cl) for cl in
range(1,num_classes+1)])+1
# if current batch already bigger than the batch random ratio, then
# check that weakest class in this patient is not the weakest in current batch (since needs to be boosted)
# also that at least one roi of this patient belongs to weakest class. If True, keep patient, else keep looking.
if (cand_rarest_class != rarest_class and np.count_nonzero(class_targets[cand] == rarest_class) > 0) \
or ix < int(batch_size * random_ratio):
break
for c in range(1,num_classes+1):
class_counts[c] += np.count_nonzero(class_targets[pick] == c)
if not ix < int(batch_size * random_ratio) and class_counts[rarest_class] == 0: # means searched thru whole set without finding rarest class
print("Class {} not represented in current dataset.".format(rarest_class))
rarest_class = np.argmin(([class_counts[c] for c in range(1,num_classes+1)]))+1
batch_patients[ix] = pick
not_picked = not_picked[not_picked != pick] # removes pick
return batch_patients
class BatchGenerator(SlimDataLoaderBase):
"""
create the training/validation batch generator. Randomly sample batch_size patients
from the data set, (draw a random slice if 2D), pad-crop them to equal sizes and merge to an array.
:param data: data dictionary as provided by 'load_dataset'
:param img_modalities: list of strings ['adc', 'b1500'] from config
:param batch_size: number of patients to sample for the batch
:param pre_crop_size: equal size for merging the patients to a single array (before the final random-crop in data aug.)
:return dictionary containing the batch data / seg / pids as lists; the augmenter will later concatenate them into an array.
"""
def __init__(self, cf, data, sample_pids_w_replace=True, max_batches=None, raise_stop_iteration=False, n_threads=None, seed=0):
if n_threads is None:
n_threads = cf.n_workers
super(BatchGenerator, self).__init__(data, cf.batch_size, number_of_threads_in_multithreaded=n_threads)
self.cf = cf
self.random_count = int(cf.batch_random_ratio * cf.batch_size)
self.plot_dir = os.path.join(self.cf.plot_dir, 'train_generator')
os.makedirs(self.plot_dir, exist_ok=True)
self.max_batches = max_batches
self.raise_stop = raise_stop_iteration
self.thread_id = 0
self.batches_produced = 0
self.dataset_length = len(self._data)
self.dataset_pids = list(self._data.keys())
- self.n_filled_threads = min(int(self.dataset_length/self.batch_size), self.number_of_threads_in_multithreaded)
- if self.n_filled_threads != self.number_of_threads_in_multithreaded:
- print("Adjusting nr of threads from {} to {}.".format(self.number_of_threads_in_multithreaded,
- self.n_filled_threads))
-
- self.rgen = np.random.RandomState(seed=seed)
- self.eligible_pids = self.rgen.permutation(self.dataset_pids.copy())
- self.eligible_pids = np.array_split(self.eligible_pids, self.n_filled_threads)
- self.eligible_pids = sorted(self.eligible_pids, key=len, reverse=True)
self.sample_pids_w_replace = sample_pids_w_replace
+ self.n_filled_threads = min(int(self.dataset_length/self.batch_size), self.number_of_threads_in_multithreaded)
if not self.sample_pids_w_replace:
+ # if not sampling w replace --> iterator-like behaviour but multi-threaded. adjust threads
+ # s.t. each thread has enough patients for at least one batch.
assert len(self.dataset_pids) / self.n_filled_threads >= self.batch_size, \
"at least one batch needed per thread. dataset size: {}, n_threads: {}, batch_size: {}.".format(
len(self.dataset_pids), self.n_filled_threads, self.batch_size)
self.lock = Lock()
+ if self.n_filled_threads != self.number_of_threads_in_multithreaded:
+ print("adjusting nr of threads from {} to {}.".format(self.number_of_threads_in_multithreaded,
+ self.n_filled_threads))
+ self.rgen = np.random.RandomState(seed=seed)
+ self.eligible_pids = self.rgen.permutation(self.dataset_pids.copy())
+ self.eligible_pids = np.array_split(self.eligible_pids, self.n_filled_threads)
+ self.eligible_pids = sorted(self.eligible_pids, key=len, reverse=True)
+
+
if hasattr(cf, "balance_target"):
# WARNING: "balance targets are only implemented for 1-d targets (or 1-component vectors)"
self.balance_target = cf.balance_target
else:
self.balance_target = "class_targets"
self.targets = {k:v[self.balance_target] for (k,v) in self._data.items()}
def set_thread_id(self, thread_id):
self.thread_ids = self.eligible_pids[thread_id]
self.thread_id = thread_id
def reset(self):
self.batches_produced = 0
self.thread_ids = self.rgen.permutation(self.eligible_pids[self.thread_id])
@staticmethod
def sample_targets_to_weights(targets, fg_bg_weights):
weights = targets * fg_bg_weights
return weights
def balance_target_distribution(self, plot=False):
"""Impose a drawing distribution over samples.
Distribution should be designed so that classes' fg and bg examples are (as good as possible) shown in
equal frequency. Since we are dealing with rois, fg/bg weights count a sample (e.g., a patient) with
**at least** one occurrence as fg, otherwise bg. For fg weights among classes, each RoI counts.
:param all_pids:
:param self.targets: dic holding {patient_specifier : patient-wise-unique ROI targets}
:return: probability distribution over all pids. draw without replace from this.
"""
- # oversampling of fg: limit bg weights to 1 s.t. bg is weighted at max as heavily as it occurs.
- max_bg_weight = 1.
+ # oversampling of fg: limit bg weights to anything <= fg weights by setting factor < 1 to overweight fg.
+ bg_weight_factor = 0.1
self.unique_ts = np.unique([v for pat in self.targets.values() for v in pat])
self.sample_stats = pd.DataFrame(columns=[str(ix)+suffix for ix in self.unique_ts for suffix in ["", "_bg"]], index=list(self.targets.keys()))
for pid in self.sample_stats.index:
for targ in self.unique_ts:
fg_count = np.count_nonzero(self.targets[pid] == targ)
self.sample_stats.loc[pid, str(targ)] = int(fg_count > 0)
self.sample_stats.loc[pid, str(targ)+"_bg"] = int(fg_count == 0)
self.targ_stats = self.sample_stats.agg(
("sum", lambda col: col.sum() / len(self._data)), axis=0, sort=False).rename({"<lambda>": "relative"})
anchor = 1. - self.targ_stats.loc["relative"].iloc[0]
self.fg_bg_weights = anchor / self.targ_stats.loc["relative"]
cum_weights = anchor * len(self.fg_bg_weights)
self.fg_bg_weights /= cum_weights
mask = ["_bg" in ix for ix in self.fg_bg_weights.index]
- self.fg_bg_weights.loc[mask] = self.fg_bg_weights.loc[mask].apply(lambda x: min(x, max_bg_weight))
+ self.fg_bg_weights.loc[mask] = self.fg_bg_weights.loc[mask].apply(lambda x: x * bg_weight_factor)
self.p_probs = self.sample_stats.apply(self.sample_targets_to_weights, args=(self.fg_bg_weights,), axis=1).sum(axis=1)
self.p_probs = self.p_probs / self.p_probs.sum()
if plot:
print("Applying class-weights:\n {}".format(self.fg_bg_weights))
self.stats = {"roi_counts": np.zeros(len(self.unique_ts,), dtype='uint32'),
"empty_counts": np.zeros(len(self.unique_ts,), dtype='uint32')}
if plot:
os.makedirs(self.plot_dir, exist_ok=True)
plg.plot_batchgen_distribution(self.cf, self.dataset_pids, self.p_probs, self.balance_target,
out_file=os.path.join(self.plot_dir,
"train_gen_distr_"+str(self.cf.fold)+".png"))
return self.p_probs
def get_batch_pids(self):
if self.max_batches is not None and self.batches_produced * self.n_filled_threads \
+ self.thread_id >= self.max_batches:
self.reset()
raise StopIteration
if self.sample_pids_w_replace:
# fully random patients
batch_pids = list(np.random.choice(self.dataset_pids, size=self.random_count, replace=False))
# target-balanced patients
batch_pids += list(np.random.choice(
self.dataset_pids, size=self.batch_size - self.random_count, replace=False, p=self.p_probs))
else:
with self.lock:
if len(self.thread_ids) == 0:
if self.raise_stop:
self.reset()
raise StopIteration
else:
self.thread_ids = self.rgen.permutation(self.eligible_pids[self.thread_id])
batch_pids = self.thread_ids[:self.batch_size]
# batch_pids = np.random.choice(self.thread_ids, size=self.batch_size, replace=False)
self.thread_ids = [pid for pid in self.thread_ids if pid not in batch_pids]
self.batches_produced += 1
return batch_pids
def generate_train_batch(self):
# to be overriden by child
# everything done in here is per batch
# print statements in here get confusing due to multithreading
raise NotImplementedError
def print_stats(self, logger=None, file=None, plot_file=None, plot=True):
print_f = utils.CombinedPrinter(logger, file)
print_f('\n***Final Training Stats***')
total_count = np.sum(self.stats['roi_counts'])
for tix, count in enumerate(self.stats['roi_counts']):
#name = self.cf.class_dict[tix] if self.balance_target=="class_targets" else str(self.unique_ts[tix])
name=str(self.unique_ts[tix])
print_f('{}: {} rois seen ({:.1f}%).'.format(name, count, count / total_count * 100))
total_samples = self.cf.num_epochs*self.cf.num_train_batches*self.cf.batch_size
empties = [
'{}: {} ({:.1f}%)'.format(str(name), self.stats['empty_counts'][tix],
self.stats['empty_counts'][tix]/total_samples*100)
for tix, name in enumerate(self.unique_ts)
]
empties = ", ".join(empties)
print_f('empty samples seen: {}\n'.format(empties))
if plot:
if plot_file is None:
plot_file = os.path.join(self.plot_dir, "train_gen_stats_{}.png".format(self.cf.fold))
os.makedirs(self.plot_dir, exist_ok=True)
plg.plot_batchgen_stats(self.cf, self.stats, empties, self.balance_target, self.unique_ts, plot_file)
class PatientBatchIterator(SlimDataLoaderBase):
"""
creates a val/test generator. Step through the dataset and return dictionaries per patient.
2D is a special case of 3D patching with patch_size[2] == 1 (slices)
Creates whole Patient batch and targets, and - if necessary - patchwise batch and targets.
Appends patient targets anyway for evaluation.
For Patching, shifts all patches into batch dimension. batch_tiling_forward will take care of exceeding batch dimensions.
This iterator/these batches are not intended to go through MTaugmenter afterwards
"""
def __init__(self, cf, data):
super(PatientBatchIterator, self).__init__(data, 0)
self.cf = cf
self.dataset_length = len(self._data)
self.dataset_pids = list(self._data.keys())
def generate_train_batch(self, pid=None):
# to be overriden by child
return
###################################
# transforms, image manipulation #
###################################
def get_patch_crop_coords(img, patch_size, min_overlap=30):
"""
_:param img (y, x, (z))
_:param patch_size: list of len 2 (2D) or 3 (3D).
_:param min_overlap: minimum required overlap of patches.
If too small, some areas are poorly represented only at edges of single patches.
_:return ndarray: shape (n_patches, 2*dim). crop coordinates for each patch.
"""
crop_coords = []
for dim in range(len(img.shape)):
n_patches = int(np.ceil(img.shape[dim] / patch_size[dim]))
# no crops required in this dimension, add image shape as coordinates.
if n_patches == 1:
crop_coords.append([(0, img.shape[dim])])
continue
# fix the two outside patches to coords patchsize/2 and interpolate.
center_dists = (img.shape[dim] - patch_size[dim]) / (n_patches - 1)
if (patch_size[dim] - center_dists) < min_overlap:
n_patches += 1
center_dists = (img.shape[dim] - patch_size[dim]) / (n_patches - 1)
patch_centers = np.round([(patch_size[dim] / 2 + (center_dists * ii)) for ii in range(n_patches)])
dim_crop_coords = [(center - patch_size[dim] / 2, center + patch_size[dim] / 2) for center in patch_centers]
crop_coords.append(dim_crop_coords)
coords_mesh_grid = []
for ymin, ymax in crop_coords[0]:
for xmin, xmax in crop_coords[1]:
if len(crop_coords) == 3 and patch_size[2] > 1:
for zmin, zmax in crop_coords[2]:
coords_mesh_grid.append([ymin, ymax, xmin, xmax, zmin, zmax])
elif len(crop_coords) == 3 and patch_size[2] == 1:
for zmin in range(img.shape[2]):
coords_mesh_grid.append([ymin, ymax, xmin, xmax, zmin, zmin + 1])
else:
coords_mesh_grid.append([ymin, ymax, xmin, xmax])
return np.array(coords_mesh_grid).astype(int)
def pad_nd_image(image, new_shape=None, mode="edge", kwargs=None, return_slicer=False, shape_must_be_divisible_by=None):
"""
one padder to pad them all. Documentation? Well okay. A little bit. by Fabian Isensee
:param image: nd image. can be anything
:param new_shape: what shape do you want? new_shape does not have to have the same dimensionality as image. If
len(new_shape) < len(image.shape) then the last axes of image will be padded. If new_shape < image.shape in any of
the axes then we will not pad that axis, but also not crop! (interpret new_shape as new_min_shape)
Example:
image.shape = (10, 1, 512, 512); new_shape = (768, 768) -> result: (10, 1, 768, 768). Cool, huh?
image.shape = (10, 1, 512, 512); new_shape = (364, 768) -> result: (10, 1, 512, 768).
:param mode: see np.pad for documentation
:param return_slicer: if True then this function will also return what coords you will need to use when cropping back
to original shape
:param shape_must_be_divisible_by: for network prediction. After applying new_shape, make sure the new shape is
divisibly by that number (can also be a list with an entry for each axis). Whatever is missing to match that will
be padded (so the result may be larger than new_shape if shape_must_be_divisible_by is not None)
:param kwargs: see np.pad for documentation
"""
if kwargs is None:
kwargs = {}
if new_shape is not None:
old_shape = np.array(image.shape[-len(new_shape):])
else:
assert shape_must_be_divisible_by is not None
assert isinstance(shape_must_be_divisible_by, (list, tuple, np.ndarray))
new_shape = image.shape[-len(shape_must_be_divisible_by):]
old_shape = new_shape
num_axes_nopad = len(image.shape) - len(new_shape)
new_shape = [max(new_shape[i], old_shape[i]) for i in range(len(new_shape))]
if not isinstance(new_shape, np.ndarray):
new_shape = np.array(new_shape)
if shape_must_be_divisible_by is not None:
if not isinstance(shape_must_be_divisible_by, (list, tuple, np.ndarray)):
shape_must_be_divisible_by = [shape_must_be_divisible_by] * len(new_shape)
else:
assert len(shape_must_be_divisible_by) == len(new_shape)
for i in range(len(new_shape)):
if new_shape[i] % shape_must_be_divisible_by[i] == 0:
new_shape[i] -= shape_must_be_divisible_by[i]
new_shape = np.array([new_shape[i] + shape_must_be_divisible_by[i] - new_shape[i] % shape_must_be_divisible_by[i] for i in range(len(new_shape))])
difference = new_shape - old_shape
pad_below = difference // 2
pad_above = difference // 2 + difference % 2
pad_list = [[0, 0]]*num_axes_nopad + list([list(i) for i in zip(pad_below, pad_above)])
res = np.pad(image, pad_list, mode, **kwargs)
if not return_slicer:
return res
else:
pad_list = np.array(pad_list)
pad_list[:, 1] = np.array(res.shape) - pad_list[:, 1]
slicer = list(slice(*i) for i in pad_list)
return res, slicer
def convert_seg_to_bounding_box_coordinates(data_dict, dim, roi_item_keys, get_rois_from_seg=False,
class_specific_seg=False):
'''adapted from batchgenerators
:param data_dict: seg: segmentation with labels indicating roi_count (get_rois_from_seg=False) or classes (get_rois_from_seg=True),
class_targets: list where list index corresponds to roi id (roi_count)
:param dim:
:param roi_item_keys: keys of the roi-wise items in data_dict to process
:param n_rg_feats: nr of regression vector features
:param get_rois_from_seg:
:return: coords (y1,x1,y2,x2 (,z1,z2)) where the segmentation GT is framed by +1 voxel, i.e., for an object with
z-extensions z1=0 through z2=5, bbox target coords will be z1=-1, z2=6. (analogically for x,y).
data_dict['roi_masks']: (b, n(b), c, h(n), w(n) (z(n))) list like roi_labels but with arrays (masks) inplace of
integers. c==1 if segmentation not one-hot encoded.
'''
bb_target = []
roi_masks = []
roi_items = {name:[] for name in roi_item_keys}
out_seg = np.copy(data_dict['seg'])
for b in range(data_dict['seg'].shape[0]):
p_coords_list = [] #p for patient?
p_roi_masks_list = []
p_roi_items_lists = {name:[] for name in roi_item_keys}
if np.sum(data_dict['seg'][b] != 0) > 0:
if get_rois_from_seg:
clusters, n_cands = lb(data_dict['seg'][b])
data_dict['class_targets'][b] = [data_dict['class_targets'][b]] * n_cands
else:
n_cands = int(np.max(data_dict['seg'][b]))
rois = np.array(
[(data_dict['seg'][b] == ii) * 1 for ii in range(1, n_cands + 1)], dtype='uint8') # separate clusters
for rix, r in enumerate(rois):
if np.sum(r != 0) > 0: # check if the roi survived slicing (3D->2D) and data augmentation (cropping etc.)
seg_ixs = np.argwhere(r != 0)
coord_list = [np.min(seg_ixs[:, 1]) - 1, np.min(seg_ixs[:, 2]) - 1, np.max(seg_ixs[:, 1]) + 1,
np.max(seg_ixs[:, 2]) + 1]
if dim == 3:
coord_list.extend([np.min(seg_ixs[:, 3]) - 1, np.max(seg_ixs[:, 3]) + 1])
p_coords_list.append(coord_list)
p_roi_masks_list.append(r)
# add background class = 0. rix is a patient wide index of lesions. since 'class_targets' is
# also patient wide, this assignment is not dependent on patch occurrences.
for name in roi_item_keys:
p_roi_items_lists[name].append(data_dict[name][b][rix])
assert data_dict["class_targets"][b][rix]>=1, "convertsegtobbox produced bg roi w cl targ {} and unique roi seg {}".format(data_dict["class_targets"][b][rix], np.unique(r))
if class_specific_seg:
out_seg[b][data_dict['seg'][b] == rix + 1] = data_dict['class_targets'][b][rix]
if not class_specific_seg:
out_seg[b][data_dict['seg'][b] > 0] = 1
bb_target.append(np.array(p_coords_list))
roi_masks.append(np.array(p_roi_masks_list))
for name in roi_item_keys:
roi_items[name].append(np.array(p_roi_items_lists[name]))
else:
bb_target.append([])
roi_masks.append(np.zeros_like(data_dict['seg'][b], dtype='uint8')[None])
for name in roi_item_keys:
roi_items[name].append(np.array([]))
if get_rois_from_seg:
data_dict.pop('class_targets', None)
data_dict['bb_target'] = np.array(bb_target)
data_dict['roi_masks'] = np.array(roi_masks)
data_dict['seg'] = out_seg
for name in roi_item_keys:
data_dict[name] = np.array(roi_items[name])
return data_dict
class ConvertSegToBoundingBoxCoordinates(AbstractTransform):
""" Converts segmentation masks into bounding box coordinates.
"""
def __init__(self, dim, roi_item_keys, get_rois_from_seg=False, class_specific_seg=False):
self.dim = dim
self.roi_item_keys = roi_item_keys
self.get_rois_from_seg = get_rois_from_seg
self.class_specific_seg = class_specific_seg
def __call__(self, **data_dict):
return convert_seg_to_bounding_box_coordinates(data_dict, self.dim, self.roi_item_keys, self.get_rois_from_seg,
self.class_specific_seg)
#############################
# data packing / unpacking # not used, data_manager.py used instead
#############################
def get_case_identifiers(folder):
case_identifiers = [i[:-4] for i in os.listdir(folder) if i.endswith("npz")]
return case_identifiers
def convert_to_npy(npz_file):
if not os.path.isfile(npz_file[:-3] + "npy"):
a = np.load(npz_file)['data']
np.save(npz_file[:-3] + "npy", a)
def unpack_dataset(folder, threads=8):
case_identifiers = get_case_identifiers(folder)
p = Pool(threads)
npz_files = [os.path.join(folder, i + ".npz") for i in case_identifiers]
p.map(convert_to_npy, npz_files)
p.close()
p.join()
def delete_npy(folder):
case_identifiers = get_case_identifiers(folder)
npy_files = [os.path.join(folder, i + ".npy") for i in case_identifiers]
npy_files = [i for i in npy_files if os.path.isfile(i)]
for n in npy_files:
os.remove(n)
\ No newline at end of file
diff --git a/utils/model_utils.py b/utils/model_utils.py
index 745cba2..4a4cfce 100644
--- a/utils/model_utils.py
+++ b/utils/model_utils.py
@@ -1,1563 +1,1567 @@
#!/usr/bin/env python
# Copyright 2019 Division of Medical Image Computing, German Cancer Research Center (DKFZ).
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
# ==============================================================================
"""
Parts are based on https://github.com/multimodallearning/pytorch-mask-rcnn
published under MIT license.
"""
import warnings
warnings.filterwarnings('ignore', '.*From scipy 0.13.0, the output shape of zoom()*')
import numpy as np
import scipy.misc
import scipy.ndimage
import scipy.interpolate
from scipy.ndimage.measurements import label as lb
import torch
import tqdm
from custom_extensions.nms import nms
from custom_extensions.roi_align import roi_align
############################################################
# Segmentation Processing
############################################################
def sum_tensor(input, axes, keepdim=False):
axes = np.unique(axes)
if keepdim:
for ax in axes:
input = input.sum(ax, keepdim=True)
else:
for ax in sorted(axes, reverse=True):
input = input.sum(int(ax))
return input
def get_one_hot_encoding(y, n_classes):
"""
transform a numpy label array to a one-hot array of the same shape.
:param y: array of shape (b, 1, y, x, (z)).
:param n_classes: int, number of classes to unfold in one-hot encoding.
:return y_ohe: array of shape (b, n_classes, y, x, (z))
"""
dim = len(y.shape) - 2
if dim == 2:
y_ohe = np.zeros((y.shape[0], n_classes, y.shape[2], y.shape[3])).astype('int32')
elif dim == 3:
y_ohe = np.zeros((y.shape[0], n_classes, y.shape[2], y.shape[3], y.shape[4])).astype('int32')
else:
raise Exception("invalid dimensions {} encountered".format(y.shape))
for cl in np.arange(n_classes):
y_ohe[:, cl][y[:, 0] == cl] = 1
return y_ohe
def dice_per_batch_inst_and_class(pred, y, n_classes, convert_to_ohe=True, smooth=1e-8):
'''
computes dice scores per batch instance and class.
:param pred: prediction array of shape (b, 1, y, x, (z)) (e.g. softmax prediction with argmax over dim 1)
:param y: ground truth array of shape (b, 1, y, x, (z)) (contains int [0, ..., n_classes]
:param n_classes: int
:return: dice scores of shape (b, c)
'''
if convert_to_ohe:
pred = get_one_hot_encoding(pred, n_classes)
y = get_one_hot_encoding(y, n_classes)
axes = tuple(range(2, len(pred.shape)))
intersect = np.sum(pred*y, axis=axes)
denominator = np.sum(pred, axis=axes)+np.sum(y, axis=axes)
dice = (2.0*intersect + smooth) / (denominator + smooth)
return dice
def dice_per_batch_and_class(pred, targ, n_classes, convert_to_ohe=True, smooth=1e-8):
'''
computes dice scores per batch and class.
:param pred: prediction array of shape (b, 1, y, x, (z)) (e.g. softmax prediction with argmax over dim 1)
:param targ: ground truth array of shape (b, 1, y, x, (z)) (contains int [0, ..., n_classes])
:param n_classes: int
:param smooth: Laplacian smooth, https://en.wikipedia.org/wiki/Additive_smoothing
:return: dice scores of shape (b, c)
'''
if convert_to_ohe:
pred = get_one_hot_encoding(pred, n_classes)
targ = get_one_hot_encoding(targ, n_classes)
axes = (0, *list(range(2, len(pred.shape)))) #(0,2,3(,4))
intersect = np.sum(pred * targ, axis=axes)
denominator = np.sum(pred, axis=axes) + np.sum(targ, axis=axes)
dice = (2.0 * intersect + smooth) / (denominator + smooth)
assert dice.shape==(n_classes,), "dice shp {}".format(dice.shape)
return dice
def batch_dice(pred, y, false_positive_weight=1.0, smooth=1e-6):
'''
compute soft dice over batch. this is a differentiable score and can be used as a loss function.
only dice scores of foreground classes are returned, since training typically
does not benefit from explicit background optimization. Pixels of the entire batch are considered a pseudo-volume to compute dice scores of.
This way, single patches with missing foreground classes can not produce faulty gradients.
:param pred: (b, c, y, x, (z)), softmax probabilities (network output).
:param y: (b, c, y, x, (z)), one hote encoded segmentation mask.
:param false_positive_weight: float [0,1]. For weighting of imbalanced classes,
reduces the penalty for false-positive pixels. Can be beneficial sometimes in data with heavy fg/bg imbalances.
:return: soft dice score (float).This function discards the background score and returns the mena of foreground scores.
'''
if len(pred.size()) == 4:
axes = (0, 2, 3)
intersect = sum_tensor(pred * y, axes, keepdim=False)
denom = sum_tensor(false_positive_weight*pred + y, axes, keepdim=False)
return torch.mean(( (2*intersect + smooth) / (denom + smooth))[1:]) #only fg dice here.
elif len(pred.size()) == 5:
axes = (0, 2, 3, 4)
intersect = sum_tensor(pred * y, axes, keepdim=False)
denom = sum_tensor(false_positive_weight*pred + y, axes, keepdim=False)
return torch.mean(( (2*intersect + smooth) / (denom + smooth))[1:]) #only fg dice here.
else:
raise ValueError('wrong input dimension in dice loss')
############################################################
# Bounding Boxes
############################################################
def compute_iou_2D(box, boxes, box_area, boxes_area):
"""Calculates IoU of the given box with the array of the given boxes.
box: 1D vector [y1, x1, y2, x2] THIS IS THE GT BOX
boxes: [boxes_count, (y1, x1, y2, x2)]
box_area: float. the area of 'box'
boxes_area: array of length boxes_count.
Note: the areas are passed in rather than calculated here for
efficency. Calculate once in the caller to avoid duplicate work.
"""
# Calculate intersection areas
y1 = np.maximum(box[0], boxes[:, 0])
y2 = np.minimum(box[2], boxes[:, 2])
x1 = np.maximum(box[1], boxes[:, 1])
x2 = np.minimum(box[3], boxes[:, 3])
intersection = np.maximum(x2 - x1, 0) * np.maximum(y2 - y1, 0)
union = box_area + boxes_area[:] - intersection[:]
iou = intersection / union
return iou
def compute_iou_3D(box, boxes, box_volume, boxes_volume):
"""Calculates IoU of the given box with the array of the given boxes.
box: 1D vector [y1, x1, y2, x2, z1, z2] (typically gt box)
boxes: [boxes_count, (y1, x1, y2, x2, z1, z2)]
box_area: float. the area of 'box'
boxes_area: array of length boxes_count.
Note: the areas are passed in rather than calculated here for
efficency. Calculate once in the caller to avoid duplicate work.
"""
# Calculate intersection areas
y1 = np.maximum(box[0], boxes[:, 0])
y2 = np.minimum(box[2], boxes[:, 2])
x1 = np.maximum(box[1], boxes[:, 1])
x2 = np.minimum(box[3], boxes[:, 3])
z1 = np.maximum(box[4], boxes[:, 4])
z2 = np.minimum(box[5], boxes[:, 5])
intersection = np.maximum(x2 - x1, 0) * np.maximum(y2 - y1, 0) * np.maximum(z2 - z1, 0)
union = box_volume + boxes_volume[:] - intersection[:]
iou = intersection / union
return iou
def compute_overlaps(boxes1, boxes2):
"""Computes IoU overlaps between two sets of boxes.
boxes1, boxes2: [N, (y1, x1, y2, x2)]. / 3D: (z1, z2))
For better performance, pass the largest set first and the smaller second.
:return: (#boxes1, #boxes2), ious of each box of 1 machted with each of 2
"""
# Areas of anchors and GT boxes
if boxes1.shape[1] == 4:
area1 = (boxes1[:, 2] - boxes1[:, 0]) * (boxes1[:, 3] - boxes1[:, 1])
area2 = (boxes2[:, 2] - boxes2[:, 0]) * (boxes2[:, 3] - boxes2[:, 1])
# Compute overlaps to generate matrix [boxes1 count, boxes2 count]
# Each cell contains the IoU value.
overlaps = np.zeros((boxes1.shape[0], boxes2.shape[0]))
for i in range(overlaps.shape[1]):
box2 = boxes2[i] #this is the gt box
overlaps[:, i] = compute_iou_2D(box2, boxes1, area2[i], area1)
return overlaps
else:
# Areas of anchors and GT boxes
volume1 = (boxes1[:, 2] - boxes1[:, 0]) * (boxes1[:, 3] - boxes1[:, 1]) * (boxes1[:, 5] - boxes1[:, 4])
volume2 = (boxes2[:, 2] - boxes2[:, 0]) * (boxes2[:, 3] - boxes2[:, 1]) * (boxes2[:, 5] - boxes2[:, 4])
# Compute overlaps to generate matrix [boxes1 count, boxes2 count]
# Each cell contains the IoU value.
overlaps = np.zeros((boxes1.shape[0], boxes2.shape[0]))
for i in range(boxes2.shape[0]):
box2 = boxes2[i] # this is the gt box
overlaps[:, i] = compute_iou_3D(box2, boxes1, volume2[i], volume1)
return overlaps
def box_refinement(box, gt_box):
"""Compute refinement needed to transform box to gt_box.
box and gt_box are [N, (y1, x1, y2, x2)] / 3D: (z1, z2))
"""
height = box[:, 2] - box[:, 0]
width = box[:, 3] - box[:, 1]
center_y = box[:, 0] + 0.5 * height
center_x = box[:, 1] + 0.5 * width
gt_height = gt_box[:, 2] - gt_box[:, 0]
gt_width = gt_box[:, 3] - gt_box[:, 1]
gt_center_y = gt_box[:, 0] + 0.5 * gt_height
gt_center_x = gt_box[:, 1] + 0.5 * gt_width
dy = (gt_center_y - center_y) / height
dx = (gt_center_x - center_x) / width
dh = torch.log(gt_height / height)
dw = torch.log(gt_width / width)
result = torch.stack([dy, dx, dh, dw], dim=1)
if box.shape[1] > 4:
depth = box[:, 5] - box[:, 4]
center_z = box[:, 4] + 0.5 * depth
gt_depth = gt_box[:, 5] - gt_box[:, 4]
gt_center_z = gt_box[:, 4] + 0.5 * gt_depth
dz = (gt_center_z - center_z) / depth
dd = torch.log(gt_depth / depth)
result = torch.stack([dy, dx, dz, dh, dw, dd], dim=1)
return result
def unmold_mask_2D(mask, bbox, image_shape):
"""Converts a mask generated by the neural network into a format similar
to it's original shape.
mask: [height, width] of type float. A small, typically 28x28 mask.
bbox: [y1, x1, y2, x2]. The box to fit the mask in.
Returns a binary mask with the same size as the original image.
"""
y1, x1, y2, x2 = bbox
out_zoom = [y2 - y1, x2 - x1]
zoom_factor = [i / j for i, j in zip(out_zoom, mask.shape)]
mask = scipy.ndimage.zoom(mask, zoom_factor, order=1).astype(np.float32)
# Put the mask in the right location.
full_mask = np.zeros(image_shape[:2]) #only y,x
full_mask[y1:y2, x1:x2] = mask
return full_mask
def unmold_mask_2D_torch(mask, bbox, image_shape):
"""Converts a mask generated by the neural network into a format similar
to it's original shape.
mask: [height, width] of type float. A small, typically 28x28 mask.
bbox: [y1, x1, y2, x2]. The box to fit the mask in.
Returns a binary mask with the same size as the original image.
"""
y1, x1, y2, x2 = bbox
out_zoom = [(y2 - y1).float(), (x2 - x1).float()]
zoom_factor = [i / j for i, j in zip(out_zoom, mask.shape)]
mask = mask.unsqueeze(0).unsqueeze(0)
mask = torch.nn.functional.interpolate(mask, scale_factor=zoom_factor)
mask = mask[0][0]
#mask = scipy.ndimage.zoom(mask.cpu().numpy(), zoom_factor, order=1).astype(np.float32)
#mask = torch.from_numpy(mask).cuda()
# Put the mask in the right location.
full_mask = torch.zeros(image_shape[:2]) # only y,x
full_mask[y1:y2, x1:x2] = mask
return full_mask
def unmold_mask_3D(mask, bbox, image_shape):
"""Converts a mask generated by the neural network into a format similar
to it's original shape.
mask: [height, width] of type float. A small, typically 28x28 mask.
bbox: [y1, x1, y2, x2, z1, z2]. The box to fit the mask in.
Returns a binary mask with the same size as the original image.
"""
y1, x1, y2, x2, z1, z2 = bbox
out_zoom = [y2 - y1, x2 - x1, z2 - z1]
zoom_factor = [i/j for i,j in zip(out_zoom, mask.shape)]
mask = scipy.ndimage.zoom(mask, zoom_factor, order=1).astype(np.float32)
# Put the mask in the right location.
full_mask = np.zeros(image_shape[:3])
full_mask[y1:y2, x1:x2, z1:z2] = mask
return full_mask
def nms_numpy(box_coords, scores, thresh):
""" non-maximum suppression on 2D or 3D boxes in numpy.
:param box_coords: [y1,x1,y2,x2 (,z1,z2)] with y1<=y2, x1<=x2, z1<=z2.
:param scores: ranking scores (higher score == higher rank) of boxes.
:param thresh: IoU threshold for clustering.
:return:
"""
y1 = box_coords[:, 0]
x1 = box_coords[:, 1]
y2 = box_coords[:, 2]
x2 = box_coords[:, 3]
assert np.all(y1 <= y2) and np.all(x1 <= x2), """"the definition of the coordinates is crucially important here:
coordinates of which maxima are taken need to be the lower coordinates"""
areas = (x2 - x1) * (y2 - y1)
is_3d = box_coords.shape[1] == 6
if is_3d: # 3-dim case
z1 = box_coords[:, 4]
z2 = box_coords[:, 5]
assert np.all(z1<=z2), """"the definition of the coordinates is crucially important here:
coordinates of which maxima are taken need to be the lower coordinates"""
areas *= (z2 - z1)
order = scores.argsort()[::-1]
keep = []
while order.size > 0: # order is the sorted index. maps order to index: order[1] = 24 means (rank1, ix 24)
i = order[0] # highest scoring element
yy1 = np.maximum(y1[i], y1[order]) # highest scoring element still in >order<, is compared to itself, that is okay.
xx1 = np.maximum(x1[i], x1[order])
yy2 = np.minimum(y2[i], y2[order])
xx2 = np.minimum(x2[i], x2[order])
h = np.maximum(0.0, yy2 - yy1)
w = np.maximum(0.0, xx2 - xx1)
inter = h * w
if is_3d:
zz1 = np.maximum(z1[i], z1[order])
zz2 = np.minimum(z2[i], z2[order])
d = np.maximum(0.0, zz2 - zz1)
inter *= d
iou = inter / (areas[i] + areas[order] - inter)
non_matches = np.nonzero(iou <= thresh)[0] # get all elements that were not matched and discard all others.
order = order[non_matches]
keep.append(i)
return keep
############################################################
# M-RCNN
############################################################
def refine_proposals(rpn_pred_probs, rpn_pred_deltas, proposal_count, batch_anchors, cf):
"""
Receives anchor scores and selects a subset to pass as proposals
to the second stage. Filtering is done based on anchor scores and
non-max suppression to remove overlaps. It also applies bounding
box refinement details to anchors.
:param rpn_pred_probs: (b, n_anchors, 2)
:param rpn_pred_deltas: (b, n_anchors, (y, x, (z), log(h), log(w), (log(d))))
:return: batch_normalized_props: Proposals in normalized coordinates (b, proposal_count, (y1, x1, y2, x2, (z1), (z2), score))
:return: batch_out_proposals: Box coords + RPN foreground scores
for monitoring/plotting (b, proposal_count, (y1, x1, y2, x2, (z1), (z2), score))
"""
std_dev = torch.from_numpy(cf.rpn_bbox_std_dev[None]).float().cuda()
norm = torch.from_numpy(cf.scale).float().cuda()
anchors = batch_anchors.clone()
batch_scores = rpn_pred_probs[:, :, 1]
# norm deltas
batch_deltas = rpn_pred_deltas * std_dev
batch_normalized_props = []
batch_out_proposals = []
# loop over batch dimension.
for ix in range(batch_scores.shape[0]):
scores = batch_scores[ix]
deltas = batch_deltas[ix]
non_nans = deltas == deltas
assert torch.all(non_nans), "deltas have nans: {}".format(deltas[~non_nans])
non_nans = anchors == anchors
assert torch.all(non_nans), "anchors have nans: {}".format(anchors[~non_nans])
# improve performance by trimming to top anchors by score
# and doing the rest on the smaller subset.
pre_nms_limit = min(cf.pre_nms_limit, anchors.size()[0])
scores, order = scores.sort(descending=True)
order = order[:pre_nms_limit]
scores = scores[:pre_nms_limit]
deltas = deltas[order, :]
# apply deltas to anchors to get refined anchors and filter with non-maximum suppression.
if batch_deltas.shape[-1] == 4:
boxes = apply_box_deltas_2D(anchors[order, :], deltas)
non_nans = boxes == boxes
assert torch.all(non_nans), "unnormalized boxes before clip/after delta apply have nans: {}".format(boxes[~non_nans])
boxes = clip_boxes_2D(boxes, cf.window)
else:
boxes = apply_box_deltas_3D(anchors[order, :], deltas)
boxes = clip_boxes_3D(boxes, cf.window)
non_nans = boxes == boxes
assert torch.all(non_nans), "unnormalized boxes before nms/after clip have nans: {}".format(boxes[~non_nans])
# boxes are y1,x1,y2,x2, torchvision-nms requires x1,y1,x2,y2, but consistent swap x<->y is irrelevant.
keep = nms.nms(boxes, scores, cf.rpn_nms_threshold)
keep = keep[:proposal_count]
boxes = boxes[keep, :]
rpn_scores = scores[keep][:, None]
# pad missing boxes with 0.
if boxes.shape[0] < proposal_count:
n_pad_boxes = proposal_count - boxes.shape[0]
zeros = torch.zeros([n_pad_boxes, boxes.shape[1]]).cuda()
boxes = torch.cat([boxes, zeros], dim=0)
zeros = torch.zeros([n_pad_boxes, rpn_scores.shape[1]]).cuda()
rpn_scores = torch.cat([rpn_scores, zeros], dim=0)
# concat box and score info for monitoring/plotting.
batch_out_proposals.append(torch.cat((boxes, rpn_scores), 1).cpu().data.numpy())
# normalize dimensions to range of 0 to 1.
non_nans = boxes == boxes
assert torch.all(non_nans), "unnormalized boxes after nms have nans: {}".format(boxes[~non_nans])
normalized_boxes = boxes / norm
where = normalized_boxes <=1
assert torch.all(where), "normalized box coords >1 found:\n {}\n".format(normalized_boxes[~where])
# add again batch dimension
batch_normalized_props.append(torch.cat((normalized_boxes, rpn_scores), 1).unsqueeze(0))
batch_normalized_props = torch.cat(batch_normalized_props)
batch_out_proposals = np.array(batch_out_proposals)
return batch_normalized_props, batch_out_proposals
def pyramid_roi_align(feature_maps, rois, pool_size, pyramid_levels, dim):
"""
Implements ROI Pooling on multiple levels of the feature pyramid.
:param feature_maps: list of feature maps, each of shape (b, c, y, x , (z))
:param rois: proposals (normalized coords.) as returned by RPN. contain info about original batch element allocation.
(n_proposals, (y1, x1, y2, x2, (z1), (z2), batch_ixs)
:param pool_size: list of poolsizes in dims: [x, y, (z)]
:param pyramid_levels: list. [0, 1, 2, ...]
:return: pooled: pooled feature map rois (n_proposals, c, poolsize_y, poolsize_x, (poolsize_z))
Output:
Pooled regions in the shape: [num_boxes, height, width, channels].
The width and height are those specific in the pool_shape in the layer
constructor.
"""
boxes = rois[:, :dim*2]
batch_ixs = rois[:, dim*2]
# Assign each ROI to a level in the pyramid based on the ROI area.
if dim == 2:
y1, x1, y2, x2 = boxes.chunk(4, dim=1)
else:
y1, x1, y2, x2, z1, z2 = boxes.chunk(6, dim=1)
h = y2 - y1
w = x2 - x1
# Equation 1 in https://arxiv.org/abs/1612.03144. Account for
# the fact that our coordinates are normalized here.
# divide sqrt(h*w) by 1 instead image_area.
roi_level = (4 + torch.log2(torch.sqrt(h*w))).round().int().clamp(pyramid_levels[0], pyramid_levels[-1])
# if Pyramid contains additional level P6, adapt the roi_level assignment accordingly.
if len(pyramid_levels) == 5:
roi_level[h*w > 0.65] = 5
# Loop through levels and apply ROI pooling to each.
pooled = []
box_to_level = []
fmap_shapes = [f.shape for f in feature_maps]
for level_ix, level in enumerate(pyramid_levels):
ix = roi_level == level
if not ix.any():
continue
ix = torch.nonzero(ix)[:, 0]
level_boxes = boxes[ix, :]
# re-assign rois to feature map of original batch element.
ind = batch_ixs[ix].int()
# Keep track of which box is mapped to which level
box_to_level.append(ix)
# Stop gradient propogation to ROI proposals
level_boxes = level_boxes.detach()
if len(pool_size) == 2:
# remap to feature map coordinate system
y_exp, x_exp = fmap_shapes[level_ix][2:] # exp = expansion
level_boxes.mul_(torch.tensor([y_exp, x_exp, y_exp, x_exp], dtype=torch.float32).cuda())
pooled_features = roi_align.roi_align_2d(feature_maps[level_ix],
torch.cat((ind.unsqueeze(1).float(), level_boxes), dim=1),
pool_size)
else:
y_exp, x_exp, z_exp = fmap_shapes[level_ix][2:]
level_boxes.mul_(torch.tensor([y_exp, x_exp, y_exp, x_exp, z_exp, z_exp], dtype=torch.float32).cuda())
pooled_features = roi_align.roi_align_3d(feature_maps[level_ix],
torch.cat((ind.unsqueeze(1).float(), level_boxes), dim=1),
pool_size)
pooled.append(pooled_features)
# Pack pooled features into one tensor
pooled = torch.cat(pooled, dim=0)
# Pack box_to_level mapping into one array and add another
# column representing the order of pooled boxes
box_to_level = torch.cat(box_to_level, dim=0)
# Rearrange pooled features to match the order of the original boxes
_, box_to_level = torch.sort(box_to_level)
pooled = pooled[box_to_level, :, :]
return pooled
def roi_align_3d_numpy(input: np.ndarray, rois, output_size: tuple,
spatial_scale: float = 1., sampling_ratio: int = -1) -> np.ndarray:
""" This fct mainly serves as a verification method for 3D CUDA implementation of RoIAlign, it's highly
inefficient due to the nested loops.
:param input: (ndarray[N, C, H, W, D]): input feature map
:param rois: list (N,K(n), 6), K(n) = nr of rois in batch-element n, single roi of format (y1,x1,y2,x2,z1,z2)
:param output_size:
:param spatial_scale:
:param sampling_ratio:
:return: (List[N, K(n), C, output_size[0], output_size[1], output_size[2]])
"""
out_height, out_width, out_depth = output_size
coord_grid = tuple([np.linspace(0, input.shape[dim] - 1, num=input.shape[dim]) for dim in range(2, 5)])
pooled_rois = [[]] * len(rois)
assert len(rois) == input.shape[0], "batch dim mismatch, rois: {}, input: {}".format(len(rois), input.shape[0])
print("Numpy 3D RoIAlign progress:", end="\n")
for b in range(input.shape[0]):
for roi in tqdm.tqdm(rois[b]):
y1, x1, y2, x2, z1, z2 = np.array(roi) * spatial_scale
roi_height = max(float(y2 - y1), 1.)
roi_width = max(float(x2 - x1), 1.)
roi_depth = max(float(z2 - z1), 1.)
if sampling_ratio <= 0:
sampling_ratio_h = int(np.ceil(roi_height / out_height))
sampling_ratio_w = int(np.ceil(roi_width / out_width))
sampling_ratio_d = int(np.ceil(roi_depth / out_depth))
else:
sampling_ratio_h = sampling_ratio_w = sampling_ratio_d = sampling_ratio # == n points per bin
bin_height = roi_height / out_height
bin_width = roi_width / out_width
bin_depth = roi_depth / out_depth
n_points = sampling_ratio_h * sampling_ratio_w * sampling_ratio_d
pooled_roi = np.empty((input.shape[1], out_height, out_width, out_depth), dtype="float32")
for chan in range(input.shape[1]):
lin_interpolator = scipy.interpolate.RegularGridInterpolator(coord_grid, input[b, chan],
method="linear")
for bin_iy in range(out_height):
for bin_ix in range(out_width):
for bin_iz in range(out_depth):
bin_val = 0.
for i in range(sampling_ratio_h):
for j in range(sampling_ratio_w):
for k in range(sampling_ratio_d):
loc_ijk = [
y1 + bin_iy * bin_height + (i + 0.5) * (bin_height / sampling_ratio_h),
x1 + bin_ix * bin_width + (j + 0.5) * (bin_width / sampling_ratio_w),
z1 + bin_iz * bin_depth + (k + 0.5) * (bin_depth / sampling_ratio_d)]
# print("loc_ijk", loc_ijk)
if not (np.any([c < -1.0 for c in loc_ijk]) or loc_ijk[0] > input.shape[2] or
loc_ijk[1] > input.shape[3] or loc_ijk[2] > input.shape[4]):
for catch_case in range(3):
# catch on-border cases
if int(loc_ijk[catch_case]) == input.shape[catch_case + 2] - 1:
loc_ijk[catch_case] = input.shape[catch_case + 2] - 1
bin_val += lin_interpolator(loc_ijk)
pooled_roi[chan, bin_iy, bin_ix, bin_iz] = bin_val / n_points
pooled_rois[b].append(pooled_roi)
return np.array(pooled_rois)
def refine_detections(cf, batch_ixs, rois, deltas, scores, regressions):
"""
Refine classified proposals (apply deltas to rpn rois), filter overlaps (nms) and return final detections.
:param rois: (n_proposals, 2 * dim) normalized boxes as proposed by RPN. n_proposals = batch_size * POST_NMS_ROIS
:param deltas: (n_proposals, n_classes, 2 * dim) box refinement deltas as predicted by mrcnn bbox regressor.
:param batch_ixs: (n_proposals) batch element assignment info for re-allocation.
:param scores: (n_proposals, n_classes) probabilities for all classes per roi as predicted by mrcnn classifier.
:param regressions: (n_proposals, n_classes, regression_features (+1 for uncertainty if predicted) regression vector
:return: result: (n_final_detections, (y1, x1, y2, x2, (z1), (z2), batch_ix, pred_class_id, pred_score, *regression vector features))
"""
# class IDs per ROI. Since scores of all classes are of interest (not just max class), all are kept at this point.
class_ids = []
fg_classes = cf.head_classes - 1
# repeat vectors to fill in predictions for all foreground classes.
for ii in range(1, fg_classes + 1):
class_ids += [ii] * rois.shape[0]
class_ids = torch.from_numpy(np.array(class_ids)).cuda()
batch_ixs = batch_ixs.repeat(fg_classes)
rois = rois.repeat(fg_classes, 1)
deltas = deltas.repeat(fg_classes, 1, 1)
scores = scores.repeat(fg_classes, 1)
regressions = regressions.repeat(fg_classes, 1, 1)
# get class-specific scores and bounding box deltas
idx = torch.arange(class_ids.size()[0]).long().cuda()
# using idx instead of slice [:,] squashes first dimension.
#len(class_ids)>scores.shape[1] --> probs is broadcasted by expansion from fg_classes-->len(class_ids)
batch_ixs = batch_ixs[idx]
deltas_specific = deltas[idx, class_ids]
class_scores = scores[idx, class_ids]
regressions = regressions[idx, class_ids]
# apply bounding box deltas. re-scale to image coordinates.
std_dev = torch.from_numpy(np.reshape(cf.rpn_bbox_std_dev, [1, cf.dim * 2])).float().cuda()
scale = torch.from_numpy(cf.scale).float().cuda()
refined_rois = apply_box_deltas_2D(rois, deltas_specific * std_dev) * scale if cf.dim == 2 else \
apply_box_deltas_3D(rois, deltas_specific * std_dev) * scale
# round and cast to int since we're dealing with pixels now
refined_rois = clip_to_window(cf.window, refined_rois)
refined_rois = torch.round(refined_rois)
# filter out low confidence boxes
keep = idx
keep_bool = (class_scores >= cf.model_min_confidence)
if not 0 in torch.nonzero(keep_bool).size():
score_keep = torch.nonzero(keep_bool)[:, 0]
pre_nms_class_ids = class_ids[score_keep]
pre_nms_rois = refined_rois[score_keep]
pre_nms_scores = class_scores[score_keep]
pre_nms_batch_ixs = batch_ixs[score_keep]
for j, b in enumerate(unique1d(pre_nms_batch_ixs)):
bixs = torch.nonzero(pre_nms_batch_ixs == b)[:, 0]
bix_class_ids = pre_nms_class_ids[bixs]
bix_rois = pre_nms_rois[bixs]
bix_scores = pre_nms_scores[bixs]
for i, class_id in enumerate(unique1d(bix_class_ids)):
ixs = torch.nonzero(bix_class_ids == class_id)[:, 0]
# nms expects boxes sorted by score.
ix_rois = bix_rois[ixs]
ix_scores = bix_scores[ixs]
ix_scores, order = ix_scores.sort(descending=True)
ix_rois = ix_rois[order, :]
class_keep = nms.nms(ix_rois, ix_scores, cf.detection_nms_threshold)
# map indices back.
class_keep = keep[score_keep[bixs[ixs[order[class_keep]]]]]
# merge indices over classes for current batch element
b_keep = class_keep if i == 0 else unique1d(torch.cat((b_keep, class_keep)))
# only keep top-k boxes of current batch-element
top_ids = class_scores[b_keep].sort(descending=True)[1][:cf.model_max_instances_per_batch_element]
b_keep = b_keep[top_ids]
# merge indices over batch elements.
batch_keep = b_keep if j == 0 else unique1d(torch.cat((batch_keep, b_keep)))
keep = batch_keep
else:
keep = torch.tensor([0]).long().cuda()
# arrange output
output = [refined_rois[keep], batch_ixs[keep].unsqueeze(1)]
output += [class_ids[keep].unsqueeze(1).float(), class_scores[keep].unsqueeze(1)]
output += [regressions[keep]]
result = torch.cat(output, dim=1)
# shape: (n_keeps, catted feats), catted feats: [0:dim*2] are box_coords, [dim*2] are batch_ics,
# [dim*2+1] are class_ids, [dim*2+2] are scores, [dim*2+3:] are regression vector features (incl uncertainty)
return result
def loss_example_mining(cf, batch_proposals, batch_gt_boxes, batch_gt_masks, batch_roi_scores,
batch_gt_class_ids, batch_gt_regressions):
"""
Subsamples proposals for mrcnn losses and generates targets. Sampling is done per batch element, seems to have positive
effects on training, as opposed to sampling over entire batch. Negatives are sampled via stochastic hard-example mining
(SHEM), where a number of negative proposals is drawn from larger pool of highest scoring proposals for stochasticity.
Scoring is obtained here as the max over all foreground probabilities as returned by mrcnn_classifier (worked better than
loss-based class-balancing methods like "online hard-example mining" or "focal loss".)
Classification-regression duality: regressions can be given along with classes (at least fg/bg, only class scores
are used for ranking).
:param batch_proposals: (n_proposals, (y1, x1, y2, x2, (z1), (z2), batch_ixs).
boxes as proposed by RPN. n_proposals here is determined by batch_size * POST_NMS_ROIS.
:param mrcnn_class_logits: (n_proposals, n_classes)
:param batch_gt_boxes: list over batch elements. Each element is a list over the corresponding roi target coordinates.
:param batch_gt_masks: list over batch elements. Each element is binary mask of shape (n_gt_rois, c, y, x, (z))
:param batch_gt_class_ids: list over batch elements. Each element is a list over the corresponding roi target labels.
if no classes predicted (only fg/bg from RPN): expected as pseudo classes [0, 1] for bg, fg.
:param batch_gt_regressions: list over b elements. Each element is a regression target vector. if None--> pseudo
:return: sample_indices: (n_sampled_rois) indices of sampled proposals to be used for loss functions.
:return: target_class_ids: (n_sampled_rois)containing target class labels of sampled proposals.
:return: target_deltas: (n_sampled_rois, 2 * dim) containing target deltas of sampled proposals for box refinement.
:return: target_masks: (n_sampled_rois, y, x, (z)) containing target masks of sampled proposals.
"""
# normalization of target coordinates
#global sample_regressions
if cf.dim == 2:
h, w = cf.patch_size
scale = torch.from_numpy(np.array([h, w, h, w])).float().cuda()
else:
h, w, z = cf.patch_size
scale = torch.from_numpy(np.array([h, w, h, w, z, z])).float().cuda()
positive_count = 0
negative_count = 0
sample_positive_indices = []
sample_negative_indices = []
sample_deltas = []
sample_masks = []
sample_class_ids = []
if batch_gt_regressions is not None:
sample_regressions = []
else:
target_regressions = torch.FloatTensor().cuda()
std_dev = torch.from_numpy(cf.bbox_std_dev).float().cuda()
# loop over batch and get positive and negative sample rois.
for b in range(len(batch_gt_boxes)):
gt_masks = torch.from_numpy(batch_gt_masks[b]).float().cuda()
gt_class_ids = torch.from_numpy(batch_gt_class_ids[b]).int().cuda()
if batch_gt_regressions is not None:
gt_regressions = torch.from_numpy(batch_gt_regressions[b]).float().cuda()
#if np.any(batch_gt_class_ids[b] > 0): # skip roi selection for no gt images.
if np.any([len(coords)>0 for coords in batch_gt_boxes[b]]):
gt_boxes = torch.from_numpy(batch_gt_boxes[b]).float().cuda() / scale
else:
gt_boxes = torch.FloatTensor().cuda()
# get proposals and indices of current batch element.
proposals = batch_proposals[batch_proposals[:, -1] == b][:, :-1]
batch_element_indices = torch.nonzero(batch_proposals[:, -1] == b).squeeze(1)
# Compute overlaps matrix [proposals, gt_boxes]
if not 0 in gt_boxes.size():
if gt_boxes.shape[1] == 4:
assert cf.dim == 2, "gt_boxes shape {} doesnt match cf.dim{}".format(gt_boxes.shape, cf.dim)
overlaps = bbox_overlaps_2D(proposals, gt_boxes)
else:
assert cf.dim == 3, "gt_boxes shape {} doesnt match cf.dim{}".format(gt_boxes.shape, cf.dim)
overlaps = bbox_overlaps_3D(proposals, gt_boxes)
# Determine positive and negative ROIs
roi_iou_max = torch.max(overlaps, dim=1)[0]
# 1. Positive ROIs are those with >= 0.5 IoU with a GT box
positive_roi_bool = roi_iou_max >= (0.5 if cf.dim == 2 else 0.3)
# 2. Negative ROIs are those with < 0.1 with every GT box.
negative_roi_bool = roi_iou_max < (0.1 if cf.dim == 2 else 0.01)
else:
positive_roi_bool = torch.FloatTensor().cuda()
negative_roi_bool = torch.from_numpy(np.array([1]*proposals.shape[0])).cuda()
# Sample Positive ROIs
if not 0 in torch.nonzero(positive_roi_bool).size():
positive_indices = torch.nonzero(positive_roi_bool).squeeze(1)
positive_samples = int(cf.train_rois_per_image * cf.roi_positive_ratio)
rand_idx = torch.randperm(positive_indices.size()[0])
rand_idx = rand_idx[:positive_samples].cuda()
positive_indices = positive_indices[rand_idx]
positive_samples = positive_indices.size()[0]
positive_rois = proposals[positive_indices, :]
# Assign positive ROIs to GT boxes.
positive_overlaps = overlaps[positive_indices, :]
roi_gt_box_assignment = torch.max(positive_overlaps, dim=1)[1]
roi_gt_boxes = gt_boxes[roi_gt_box_assignment, :]
roi_gt_class_ids = gt_class_ids[roi_gt_box_assignment]
if batch_gt_regressions is not None:
roi_gt_regressions = gt_regressions[roi_gt_box_assignment]
# Compute bbox refinement targets for positive ROIs
deltas = box_refinement(positive_rois, roi_gt_boxes)
deltas /= std_dev
roi_masks = gt_masks[roi_gt_box_assignment]
assert roi_masks.shape[1] == 1, "gt masks have more than one channel --> is this desired?"
# Compute mask targets
boxes = positive_rois
box_ids = torch.arange(roi_masks.shape[0]).cuda().unsqueeze(1).float()
if len(cf.mask_shape) == 2:
y_exp, x_exp = roi_masks.shape[2:] # exp = expansion
boxes.mul_(torch.tensor([y_exp, x_exp, y_exp, x_exp], dtype=torch.float32).cuda())
masks = roi_align.roi_align_2d(roi_masks,
torch.cat((box_ids, boxes), dim=1),
cf.mask_shape)
else:
y_exp, x_exp, z_exp = roi_masks.shape[2:] # exp = expansion
boxes.mul_(torch.tensor([y_exp, x_exp, y_exp, x_exp, z_exp, z_exp], dtype=torch.float32).cuda())
masks = roi_align.roi_align_3d(roi_masks,
torch.cat((box_ids, boxes), dim=1),
cf.mask_shape)
masks = masks.squeeze(1)
# Threshold mask pixels at 0.5 to have GT masks be 0 or 1 to use with
# binary cross entropy loss.
masks = torch.round(masks)
sample_positive_indices.append(batch_element_indices[positive_indices])
sample_deltas.append(deltas)
sample_masks.append(masks)
sample_class_ids.append(roi_gt_class_ids)
if batch_gt_regressions is not None:
sample_regressions.append(roi_gt_regressions)
positive_count += positive_samples
else:
positive_samples = 0
# Sample negative ROIs. Add enough to maintain positive:negative ratio, but at least 1. Sample via SHEM.
if not 0 in torch.nonzero(negative_roi_bool).size():
negative_indices = torch.nonzero(negative_roi_bool).squeeze(1)
r = 1.0 / cf.roi_positive_ratio
b_neg_count = np.max((int(r * positive_samples - positive_samples), 1))
roi_scores_neg = batch_roi_scores[batch_element_indices[negative_indices]]
raw_sampled_indices = shem(roi_scores_neg, b_neg_count, cf.shem_poolsize)
sample_negative_indices.append(batch_element_indices[negative_indices[raw_sampled_indices]])
negative_count += raw_sampled_indices.size()[0]
if len(sample_positive_indices) > 0:
target_deltas = torch.cat(sample_deltas)
target_masks = torch.cat(sample_masks)
target_class_ids = torch.cat(sample_class_ids)
if batch_gt_regressions is not None:
target_regressions = torch.cat(sample_regressions)
# Pad target information with zeros for negative ROIs.
if positive_count > 0 and negative_count > 0:
sample_indices = torch.cat((torch.cat(sample_positive_indices), torch.cat(sample_negative_indices)), dim=0)
zeros = torch.zeros(negative_count, cf.dim * 2).cuda()
target_deltas = torch.cat([target_deltas, zeros], dim=0)
zeros = torch.zeros(negative_count, *cf.mask_shape).cuda()
target_masks = torch.cat([target_masks, zeros], dim=0)
zeros = torch.zeros(negative_count).int().cuda()
target_class_ids = torch.cat([target_class_ids, zeros], dim=0)
if batch_gt_regressions is not None:
# regression targets need to have 0 as background/negative with below practice
if 'regression_bin' in cf.prediction_tasks:
zeros = torch.zeros(negative_count, dtype=torch.float).cuda()
else:
zeros = torch.zeros(negative_count, cf.regression_n_features, dtype=torch.float).cuda()
target_regressions = torch.cat([target_regressions, zeros], dim=0)
elif positive_count > 0:
sample_indices = torch.cat(sample_positive_indices)
elif negative_count > 0:
sample_indices = torch.cat(sample_negative_indices)
target_deltas = torch.zeros(negative_count, cf.dim * 2).cuda()
target_masks = torch.zeros(negative_count, *cf.mask_shape).cuda()
target_class_ids = torch.zeros(negative_count).int().cuda()
if batch_gt_regressions is not None:
if 'regression_bin' in cf.prediction_tasks:
target_regressions = torch.zeros(negative_count, dtype=torch.float).cuda()
else:
target_regressions = torch.zeros(negative_count, cf.regression_n_features, dtype=torch.float).cuda()
else:
sample_indices = torch.LongTensor().cuda()
target_class_ids = torch.IntTensor().cuda()
target_deltas = torch.FloatTensor().cuda()
target_masks = torch.FloatTensor().cuda()
target_regressions = torch.FloatTensor().cuda()
return sample_indices, target_deltas, target_masks, target_class_ids, target_regressions
############################################################
# Anchors
############################################################
def generate_anchors(scales, ratios, shape, feature_stride, anchor_stride):
"""
scales: 1D array of anchor sizes in pixels. Example: [32, 64, 128]
ratios: 1D array of anchor ratios of width/height. Example: [0.5, 1, 2]
shape: [height, width] spatial shape of the feature map over which
to generate anchors.
feature_stride: Stride of the feature map relative to the image in pixels.
anchor_stride: Stride of anchors on the feature map. For example, if the
value is 2 then generate anchors for every other feature map pixel.
"""
# Get all combinations of scales and ratios
scales, ratios = np.meshgrid(np.array(scales), np.array(ratios))
scales = scales.flatten()
ratios = ratios.flatten()
# Enumerate heights and widths from scales and ratios
heights = scales / np.sqrt(ratios)
widths = scales * np.sqrt(ratios)
# Enumerate shifts in feature space
shifts_y = np.arange(0, shape[0], anchor_stride) * feature_stride
shifts_x = np.arange(0, shape[1], anchor_stride) * feature_stride
shifts_x, shifts_y = np.meshgrid(shifts_x, shifts_y)
# Enumerate combinations of shifts, widths, and heights
box_widths, box_centers_x = np.meshgrid(widths, shifts_x)
box_heights, box_centers_y = np.meshgrid(heights, shifts_y)
# Reshape to get a list of (y, x) and a list of (h, w)
box_centers = np.stack([box_centers_y, box_centers_x], axis=2).reshape([-1, 2])
box_sizes = np.stack([box_heights, box_widths], axis=2).reshape([-1, 2])
# Convert to corner coordinates (y1, x1, y2, x2)
boxes = np.concatenate([box_centers - 0.5 * box_sizes, box_centers + 0.5 * box_sizes], axis=1)
return boxes
def generate_anchors_3D(scales_xy, scales_z, ratios, shape, feature_stride_xy, feature_stride_z, anchor_stride):
"""
scales: 1D array of anchor sizes in pixels. Example: [32, 64, 128]
ratios: 1D array of anchor ratios of width/height. Example: [0.5, 1, 2]
shape: [height, width] spatial shape of the feature map over which
to generate anchors.
feature_stride: Stride of the feature map relative to the image in pixels.
anchor_stride: Stride of anchors on the feature map. For example, if the
value is 2 then generate anchors for every other feature map pixel.
"""
# Get all combinations of scales and ratios
scales_xy, ratios_meshed = np.meshgrid(np.array(scales_xy), np.array(ratios))
scales_xy = scales_xy.flatten()
ratios_meshed = ratios_meshed.flatten()
# Enumerate heights and widths from scales and ratios
heights = scales_xy / np.sqrt(ratios_meshed)
widths = scales_xy * np.sqrt(ratios_meshed)
depths = np.tile(np.array(scales_z), len(ratios_meshed)//np.array(scales_z)[..., None].shape[0])
# Enumerate shifts in feature space
shifts_y = np.arange(0, shape[0], anchor_stride) * feature_stride_xy #translate from fm positions to input coords.
shifts_x = np.arange(0, shape[1], anchor_stride) * feature_stride_xy
shifts_z = np.arange(0, shape[2], anchor_stride) * (feature_stride_z)
shifts_x, shifts_y, shifts_z = np.meshgrid(shifts_x, shifts_y, shifts_z)
# Enumerate combinations of shifts, widths, and heights
box_widths, box_centers_x = np.meshgrid(widths, shifts_x)
box_heights, box_centers_y = np.meshgrid(heights, shifts_y)
box_depths, box_centers_z = np.meshgrid(depths, shifts_z)
# Reshape to get a list of (y, x, z) and a list of (h, w, d)
box_centers = np.stack(
[box_centers_y, box_centers_x, box_centers_z], axis=2).reshape([-1, 3])
box_sizes = np.stack([box_heights, box_widths, box_depths], axis=2).reshape([-1, 3])
# Convert to corner coordinates (y1, x1, y2, x2, z1, z2)
boxes = np.concatenate([box_centers - 0.5 * box_sizes,
box_centers + 0.5 * box_sizes], axis=1)
boxes = np.transpose(np.array([boxes[:, 0], boxes[:, 1], boxes[:, 3], boxes[:, 4], boxes[:, 2], boxes[:, 5]]), axes=(1, 0))
return boxes
def generate_pyramid_anchors(logger, cf):
"""Generate anchors at different levels of a feature pyramid. Each scale
is associated with a level of the pyramid, but each ratio is used in
all levels of the pyramid.
from configs:
:param scales: cf.RPN_ANCHOR_SCALES , for conformity with retina nets: scale entries need to be list, e.g. [[4], [8], [16], [32]]
:param ratios: cf.RPN_ANCHOR_RATIOS , e.g. [0.5, 1, 2]
:param feature_shapes: cf.BACKBONE_SHAPES , e.g. [array of shapes per feature map] [80, 40, 20, 10, 5]
:param feature_strides: cf.BACKBONE_STRIDES , e.g. [2, 4, 8, 16, 32, 64]
:param anchors_stride: cf.RPN_ANCHOR_STRIDE , e.g. 1
:return anchors: (N, (y1, x1, y2, x2, (z1), (z2)). All generated anchors in one array. Sorted
with the same order of the given scales. So, anchors of scale[0] come first, then anchors of scale[1], and so on.
"""
scales = cf.rpn_anchor_scales
ratios = cf.rpn_anchor_ratios
feature_shapes = cf.backbone_shapes
anchor_stride = cf.rpn_anchor_stride
pyramid_levels = cf.pyramid_levels
feature_strides = cf.backbone_strides
logger.info("anchor scales {} and feature map shapes {}".format(scales, feature_shapes))
expected_anchors = [np.prod(feature_shapes[level]) * len(ratios) * len(scales['xy'][level]) for level in pyramid_levels]
anchors = []
for lix, level in enumerate(pyramid_levels):
if len(feature_shapes[level]) == 2:
anchors.append(generate_anchors(scales['xy'][level], ratios, feature_shapes[level],
feature_strides['xy'][level], anchor_stride))
elif len(feature_shapes[level]) == 3:
anchors.append(generate_anchors_3D(scales['xy'][level], scales['z'][level], ratios, feature_shapes[level],
feature_strides['xy'][level], feature_strides['z'][level], anchor_stride))
else:
raise Exception("invalid feature_shapes[{}] size {}".format(level, feature_shapes[level]))
logger.info("level {}: expected anchors {}, built anchors {}.".format(level, expected_anchors[lix], anchors[-1].shape))
out_anchors = np.concatenate(anchors, axis=0)
logger.info("Total: expected anchors {}, built anchors {}.".format(np.sum(expected_anchors), out_anchors.shape))
return out_anchors
def apply_box_deltas_2D(boxes, deltas):
"""Applies the given deltas to the given boxes.
boxes: [N, 4] where each row is y1, x1, y2, x2
deltas: [N, 4] where each row is [dy, dx, log(dh), log(dw)]
"""
# Convert to y, x, h, w
non_nans = boxes == boxes
assert torch.all(non_nans), "boxes at beginning of delta apply have nans: {}".format(
boxes[~non_nans])
height = boxes[:, 2] - boxes[:, 0]
width = boxes[:, 3] - boxes[:, 1]
center_y = boxes[:, 0] + 0.5 * height
center_x = boxes[:, 1] + 0.5 * width
# Apply deltas
center_y += deltas[:, 0] * height
center_x += deltas[:, 1] * width
- height *= torch.exp(deltas[:, 2])
- width *= torch.exp(deltas[:, 3])
+ # clip delta preds in order to avoid infs and later nans after exponentiation.
+ height *= torch.exp(torch.clamp(deltas[:, 2], max=6.))
+ width *= torch.exp(torch.clamp(deltas[:, 3], max=6.))
+
+ # height *= torch.exp(deltas[:, 2])
+ # width *= torch.exp(deltas[:, 3])
non_nans = width == width
assert torch.all(non_nans), "inside delta apply, width has nans: {}".format(
width[~non_nans])
# 0.*inf results in nan. fix nans to zeros.
- height[height!=height] = 0.
- width[width!=width] = 0.
+ # height[height!=height] = 0.
+ # width[width!=width] = 0.
non_nans = height == height
assert torch.all(non_nans), "inside delta apply, height has nans directly after setting to zero: {}".format(
height[~non_nans])
non_nans = width == width
assert torch.all(non_nans), "inside delta apply, width has nans directly after setting to zero: {}".format(
width[~non_nans])
# Convert back to y1, x1, y2, x2
y1 = center_y - 0.5 * height
x1 = center_x - 0.5 * width
y2 = y1 + height
x2 = x1 + width
result = torch.stack([y1, x1, y2, x2], dim=1)
non_nans = result == result
assert torch.all(non_nans), "inside delta apply, result has nans: {}".format(result[~non_nans])
return result
def apply_box_deltas_3D(boxes, deltas):
"""Applies the given deltas to the given boxes.
boxes: [N, 6] where each row is y1, x1, y2, x2, z1, z2
deltas: [N, 6] where each row is [dy, dx, dz, log(dh), log(dw), log(dd)]
"""
# Convert to y, x, h, w
height = boxes[:, 2] - boxes[:, 0]
width = boxes[:, 3] - boxes[:, 1]
depth = boxes[:, 5] - boxes[:, 4]
center_y = boxes[:, 0] + 0.5 * height
center_x = boxes[:, 1] + 0.5 * width
center_z = boxes[:, 4] + 0.5 * depth
# Apply deltas
center_y += deltas[:, 0] * height
center_x += deltas[:, 1] * width
center_z += deltas[:, 2] * depth
height *= torch.exp(deltas[:, 3])
width *= torch.exp(deltas[:, 4])
depth *= torch.exp(deltas[:, 5])
# Convert back to y1, x1, y2, x2
y1 = center_y - 0.5 * height
x1 = center_x - 0.5 * width
z1 = center_z - 0.5 * depth
y2 = y1 + height
x2 = x1 + width
z2 = z1 + depth
result = torch.stack([y1, x1, y2, x2, z1, z2], dim=1)
return result
def clip_boxes_2D(boxes, window):
"""
boxes: [N, 4] each col is y1, x1, y2, x2
window: [4] in the form y1, x1, y2, x2
"""
boxes = torch.stack( \
[boxes[:, 0].clamp(float(window[0]), float(window[2])),
boxes[:, 1].clamp(float(window[1]), float(window[3])),
boxes[:, 2].clamp(float(window[0]), float(window[2])),
boxes[:, 3].clamp(float(window[1]), float(window[3]))], 1)
return boxes
def clip_boxes_3D(boxes, window):
"""
boxes: [N, 6] each col is y1, x1, y2, x2, z1, z2
window: [6] in the form y1, x1, y2, x2, z1, z2
"""
boxes = torch.stack( \
[boxes[:, 0].clamp(float(window[0]), float(window[2])),
boxes[:, 1].clamp(float(window[1]), float(window[3])),
boxes[:, 2].clamp(float(window[0]), float(window[2])),
boxes[:, 3].clamp(float(window[1]), float(window[3])),
boxes[:, 4].clamp(float(window[4]), float(window[5])),
boxes[:, 5].clamp(float(window[4]), float(window[5]))], 1)
return boxes
from matplotlib import pyplot as plt
def clip_boxes_numpy(boxes, window):
"""
boxes: [N, 4] each col is y1, x1, y2, x2 / [N, 6] in 3D.
window: iamge shape (y, x, (z))
"""
if boxes.shape[1] == 4:
boxes = np.concatenate(
(np.clip(boxes[:, 0], 0, window[0])[:, None],
np.clip(boxes[:, 1], 0, window[0])[:, None],
np.clip(boxes[:, 2], 0, window[1])[:, None],
np.clip(boxes[:, 3], 0, window[1])[:, None]), 1
)
else:
boxes = np.concatenate(
(np.clip(boxes[:, 0], 0, window[0])[:, None],
np.clip(boxes[:, 1], 0, window[0])[:, None],
np.clip(boxes[:, 2], 0, window[1])[:, None],
np.clip(boxes[:, 3], 0, window[1])[:, None],
np.clip(boxes[:, 4], 0, window[2])[:, None],
np.clip(boxes[:, 5], 0, window[2])[:, None]), 1
)
return boxes
def bbox_overlaps_2D(boxes1, boxes2):
"""Computes IoU overlaps between two sets of boxes.
boxes1, boxes2: [N, (y1, x1, y2, x2)].
"""
# 1. Tile boxes2 and repeate boxes1. This allows us to compare
# every boxes1 against every boxes2 without loops.
# TF doesn't have an equivalent to np.repeate() so simulate it
# using tf.tile() and tf.reshape.
boxes1_repeat = boxes2.size()[0]
boxes2_repeat = boxes1.size()[0]
boxes1 = boxes1.repeat(1,boxes1_repeat).view(-1,4)
boxes2 = boxes2.repeat(boxes2_repeat,1)
# 2. Compute intersections
b1_y1, b1_x1, b1_y2, b1_x2 = boxes1.chunk(4, dim=1)
b2_y1, b2_x1, b2_y2, b2_x2 = boxes2.chunk(4, dim=1)
y1 = torch.max(b1_y1, b2_y1)[:, 0]
x1 = torch.max(b1_x1, b2_x1)[:, 0]
y2 = torch.min(b1_y2, b2_y2)[:, 0]
x2 = torch.min(b1_x2, b2_x2)[:, 0]
#--> expects x1<x2 & y1<y2
zeros = torch.zeros(y1.size()[0], requires_grad=False)
if y1.is_cuda:
zeros = zeros.cuda()
intersection = torch.max(x2 - x1, zeros) * torch.max(y2 - y1, zeros)
# 3. Compute unions
b1_area = (b1_y2 - b1_y1) * (b1_x2 - b1_x1)
b2_area = (b2_y2 - b2_y1) * (b2_x2 - b2_x1)
union = b1_area[:,0] + b2_area[:,0] - intersection
# 4. Compute IoU and reshape to [boxes1, boxes2]
iou = intersection / union
assert torch.all(iou<=1), "iou score>1 produced in bbox_overlaps_2D"
overlaps = iou.view(boxes2_repeat, boxes1_repeat) #--> per gt box: ious of all proposal boxes with that gt box
return overlaps
def bbox_overlaps_3D(boxes1, boxes2):
"""Computes IoU overlaps between two sets of boxes.
boxes1, boxes2: [N, (y1, x1, y2, x2, z1, z2)].
"""
# 1. Tile boxes2 and repeate boxes1. This allows us to compare
# every boxes1 against every boxes2 without loops.
# TF doesn't have an equivalent to np.repeate() so simulate it
# using tf.tile() and tf.reshape.
boxes1_repeat = boxes2.size()[0]
boxes2_repeat = boxes1.size()[0]
boxes1 = boxes1.repeat(1,boxes1_repeat).view(-1,6)
boxes2 = boxes2.repeat(boxes2_repeat,1)
# 2. Compute intersections
b1_y1, b1_x1, b1_y2, b1_x2, b1_z1, b1_z2 = boxes1.chunk(6, dim=1)
b2_y1, b2_x1, b2_y2, b2_x2, b2_z1, b2_z2 = boxes2.chunk(6, dim=1)
y1 = torch.max(b1_y1, b2_y1)[:, 0]
x1 = torch.max(b1_x1, b2_x1)[:, 0]
y2 = torch.min(b1_y2, b2_y2)[:, 0]
x2 = torch.min(b1_x2, b2_x2)[:, 0]
z1 = torch.max(b1_z1, b2_z1)[:, 0]
z2 = torch.min(b1_z2, b2_z2)[:, 0]
zeros = torch.zeros(y1.size()[0], requires_grad=False)
if y1.is_cuda:
zeros = zeros.cuda()
intersection = torch.max(x2 - x1, zeros) * torch.max(y2 - y1, zeros) * torch.max(z2 - z1, zeros)
# 3. Compute unions
b1_volume = (b1_y2 - b1_y1) * (b1_x2 - b1_x1) * (b1_z2 - b1_z1)
b2_volume = (b2_y2 - b2_y1) * (b2_x2 - b2_x1) * (b2_z2 - b2_z1)
union = b1_volume[:,0] + b2_volume[:,0] - intersection
# 4. Compute IoU and reshape to [boxes1, boxes2]
iou = intersection / union
overlaps = iou.view(boxes2_repeat, boxes1_repeat)
return overlaps
def gt_anchor_matching(cf, anchors, gt_boxes, gt_class_ids=None):
"""Given the anchors and GT boxes, compute overlaps and identify positive
anchors and deltas to refine them to match their corresponding GT boxes.
anchors: [num_anchors, (y1, x1, y2, x2, (z1), (z2))]
gt_boxes: [num_gt_boxes, (y1, x1, y2, x2, (z1), (z2))]
gt_class_ids (optional): [num_gt_boxes] Integer class IDs for one stage detectors. in RPN case of Mask R-CNN,
set all positive matches to 1 (foreground)
Returns:
anchor_class_matches: [N] (int32) matches between anchors and GT boxes.
1 = positive anchor, -1 = negative anchor, 0 = neutral
anchor_delta_targets: [N, (dy, dx, (dz), log(dh), log(dw), (log(dd)))] Anchor bbox deltas.
"""
anchor_class_matches = np.zeros([anchors.shape[0]], dtype=np.int32)
anchor_delta_targets = np.zeros((cf.rpn_train_anchors_per_image, 2*cf.dim))
anchor_matching_iou = cf.anchor_matching_iou
if gt_boxes is None:
anchor_class_matches = np.full(anchor_class_matches.shape, fill_value=-1)
return anchor_class_matches, anchor_delta_targets
# for mrcnn: anchor matching is done for RPN loss, so positive labels are all 1 (foreground)
if gt_class_ids is None:
gt_class_ids = np.array([1] * len(gt_boxes))
# Compute overlaps [num_anchors, num_gt_boxes]
overlaps = compute_overlaps(anchors, gt_boxes)
# Match anchors to GT Boxes
# If an anchor overlaps a GT box with IoU >= anchor_matching_iou then it's positive.
# If an anchor overlaps a GT box with IoU < 0.1 then it's negative.
# Neutral anchors are those that don't match the conditions above,
# and they don't influence the loss function.
# However, don't keep any GT box unmatched (rare, but happens). Instead,
# match it to the closest anchor (even if its max IoU is < 0.1).
# 1. Set negative anchors first. They get overwritten below if a GT box is
# matched to them. Skip boxes in crowd areas.
anchor_iou_argmax = np.argmax(overlaps, axis=1)
anchor_iou_max = overlaps[np.arange(overlaps.shape[0]), anchor_iou_argmax]
if anchors.shape[1] == 4:
anchor_class_matches[(anchor_iou_max < 0.1)] = -1
elif anchors.shape[1] == 6:
anchor_class_matches[(anchor_iou_max < 0.01)] = -1
else:
raise ValueError('anchor shape wrong {}'.format(anchors.shape))
# 2. Set an anchor for each GT box (regardless of IoU value).
gt_iou_argmax = np.argmax(overlaps, axis=0)
for ix, ii in enumerate(gt_iou_argmax):
anchor_class_matches[ii] = gt_class_ids[ix]
# 3. Set anchors with high overlap as positive.
above_thresh_ixs = np.argwhere(anchor_iou_max >= anchor_matching_iou)
anchor_class_matches[above_thresh_ixs] = gt_class_ids[anchor_iou_argmax[above_thresh_ixs]]
# Subsample to balance positive anchors.
ids = np.where(anchor_class_matches > 0)[0]
extra = len(ids) - (cf.rpn_train_anchors_per_image // 2)
if extra > 0:
# Reset the extra ones to neutral
ids = np.random.choice(ids, extra, replace=False)
anchor_class_matches[ids] = 0
# Leave all negative proposals negative for now and sample from them later in online hard example mining.
# For positive anchors, compute shift and scale needed to transform them to match the corresponding GT boxes.
ids = np.where(anchor_class_matches > 0)[0]
ix = 0 # index into anchor_delta_targets
for i, a in zip(ids, anchors[ids]):
# closest gt box (it might have IoU < anchor_matching_iou)
gt = gt_boxes[anchor_iou_argmax[i]]
# convert coordinates to center plus width/height.
gt_h = gt[2] - gt[0]
gt_w = gt[3] - gt[1]
gt_center_y = gt[0] + 0.5 * gt_h
gt_center_x = gt[1] + 0.5 * gt_w
# Anchor
a_h = a[2] - a[0]
a_w = a[3] - a[1]
a_center_y = a[0] + 0.5 * a_h
a_center_x = a[1] + 0.5 * a_w
if cf.dim == 2:
anchor_delta_targets[ix] = [
(gt_center_y - a_center_y) / a_h,
(gt_center_x - a_center_x) / a_w,
np.log(gt_h / a_h),
np.log(gt_w / a_w),
]
else:
gt_d = gt[5] - gt[4]
gt_center_z = gt[4] + 0.5 * gt_d
a_d = a[5] - a[4]
a_center_z = a[4] + 0.5 * a_d
anchor_delta_targets[ix] = [
(gt_center_y - a_center_y) / a_h,
(gt_center_x - a_center_x) / a_w,
(gt_center_z - a_center_z) / a_d,
np.log(gt_h / a_h),
np.log(gt_w / a_w),
np.log(gt_d / a_d)
]
# normalize.
anchor_delta_targets[ix] /= cf.rpn_bbox_std_dev
ix += 1
return anchor_class_matches, anchor_delta_targets
def clip_to_window(window, boxes):
"""
window: (y1, x1, y2, x2) / 3D: (z1, z2). The window in the image we want to clip to.
boxes: [N, (y1, x1, y2, x2)] / 3D: (z1, z2)
"""
boxes[:, 0] = boxes[:, 0].clamp(float(window[0]), float(window[2]))
boxes[:, 1] = boxes[:, 1].clamp(float(window[1]), float(window[3]))
boxes[:, 2] = boxes[:, 2].clamp(float(window[0]), float(window[2]))
boxes[:, 3] = boxes[:, 3].clamp(float(window[1]), float(window[3]))
if boxes.shape[1] > 5:
boxes[:, 4] = boxes[:, 4].clamp(float(window[4]), float(window[5]))
boxes[:, 5] = boxes[:, 5].clamp(float(window[4]), float(window[5]))
return boxes
############################################################
# Connected Componenent Analysis
############################################################
def get_coords(binary_mask, n_components, dim):
"""
loops over batch to perform connected component analysis on binary input mask. computes box coordinates around
n_components - biggest components (rois).
:param binary_mask: (b, y, x, (z)). binary mask for one specific foreground class.
:param n_components: int. number of components to extract per batch element and class.
:return: coords (b, n, (y1, x1, y2, x2 (,z1, z2))
:return: batch_components (b, n, (y1, x1, y2, x2, (z1), (z2))
"""
assert len(binary_mask.shape)==dim+1
binary_mask = binary_mask.astype('uint8')
batch_coords = []
batch_components = []
for ix,b in enumerate(binary_mask):
clusters, n_cands = lb(b) # performs connected component analysis.
uniques, counts = np.unique(clusters, return_counts=True)
keep_uniques = uniques[1:][np.argsort(counts[1:])[::-1]][:n_components] #only keep n_components largest components
p_components = np.array([(clusters == ii) * 1 for ii in keep_uniques]) # separate clusters and concat
p_coords = []
if p_components.shape[0] > 0:
for roi in p_components:
mask_ixs = np.argwhere(roi != 0)
# get coordinates around component.
roi_coords = [np.min(mask_ixs[:, 0]) - 1, np.min(mask_ixs[:, 1]) - 1, np.max(mask_ixs[:, 0]) + 1,
np.max(mask_ixs[:, 1]) + 1]
if dim == 3:
roi_coords += [np.min(mask_ixs[:, 2]), np.max(mask_ixs[:, 2])+1]
p_coords.append(roi_coords)
p_coords = np.array(p_coords)
#clip coords.
p_coords[p_coords < 0] = 0
p_coords[:, :4][p_coords[:, :4] > binary_mask.shape[-2]] = binary_mask.shape[-2]
if dim == 3:
p_coords[:, 4:][p_coords[:, 4:] > binary_mask.shape[-1]] = binary_mask.shape[-1]
batch_coords.append(p_coords)
batch_components.append(p_components)
return batch_coords, batch_components
# noinspection PyCallingNonCallable
def get_coords_gpu(binary_mask, n_components, dim):
"""
loops over batch to perform connected component analysis on binary input mask. computes box coordiantes around
n_components - biggest components (rois).
:param binary_mask: (b, y, x, (z)). binary mask for one specific foreground class.
:param n_components: int. number of components to extract per batch element and class.
:return: coords (b, n, (y1, x1, y2, x2 (,z1, z2))
:return: batch_components (b, n, (y1, x1, y2, x2, (z1), (z2))
"""
raise Exception("throws floating point exception")
assert len(binary_mask.shape)==dim+1
binary_mask = binary_mask.type(torch.uint8)
batch_coords = []
batch_components = []
for ix,b in enumerate(binary_mask):
clusters, n_cands = lb(b.cpu().data.numpy()) # peforms connected component analysis.
clusters = torch.from_numpy(clusters).cuda()
uniques = torch.unique(clusters)
counts = torch.stack([(clusters==unique).sum() for unique in uniques])
keep_uniques = uniques[1:][torch.sort(counts[1:])[1].flip(0)][:n_components] #only keep n_components largest components
p_components = torch.cat([(clusters == ii).unsqueeze(0) for ii in keep_uniques]).cuda() # separate clusters and concat
p_coords = []
if p_components.shape[0] > 0:
for roi in p_components:
mask_ixs = torch.nonzero(roi)
# get coordinates around component.
roi_coords = [torch.min(mask_ixs[:, 0]) - 1, torch.min(mask_ixs[:, 1]) - 1,
torch.max(mask_ixs[:, 0]) + 1,
torch.max(mask_ixs[:, 1]) + 1]
if dim == 3:
roi_coords += [torch.min(mask_ixs[:, 2]), torch.max(mask_ixs[:, 2])+1]
p_coords.append(roi_coords)
p_coords = torch.tensor(p_coords)
#clip coords.
p_coords[p_coords < 0] = 0
p_coords[:, :4][p_coords[:, :4] > binary_mask.shape[-2]] = binary_mask.shape[-2]
if dim == 3:
p_coords[:, 4:][p_coords[:, 4:] > binary_mask.shape[-1]] = binary_mask.shape[-1]
batch_coords.append(p_coords)
batch_components.append(p_components)
return batch_coords, batch_components
############################################################
# Pytorch Utility Functions
############################################################
def unique1d(tensor):
"""discard all elements of tensor that occur more than once; make tensor unique.
:param tensor:
:return:
"""
if tensor.size()[0] == 0 or tensor.size()[0] == 1:
return tensor
tensor = tensor.sort()[0]
unique_bool = tensor[1:] != tensor[:-1]
first_element = torch.tensor([True], dtype=torch.bool, requires_grad=False)
if tensor.is_cuda:
first_element = first_element.cuda()
unique_bool = torch.cat((first_element, unique_bool), dim=0)
return tensor[unique_bool.data]
def intersect1d(tensor1, tensor2):
aux = torch.cat((tensor1, tensor2), dim=0)
aux = aux.sort(descending=True)[0]
return aux[:-1][(aux[1:] == aux[:-1]).data]
def shem(roi_probs_neg, negative_count, poolsize):
"""
stochastic hard example mining: from a list of indices (referring to non-matched predictions),
determine a pool of highest scoring (worst false positives) of size negative_count*poolsize.
Then, sample n (= negative_count) predictions of this pool as negative examples for loss.
:param roi_probs_neg: tensor of shape (n_predictions, n_classes).
:param negative_count: int.
:param poolsize: int.
:return: (negative_count). indices refer to the positions in roi_probs_neg. If pool smaller than expected due to
limited negative proposals availabel, this function will return sampled indices of number < negative_count without
throwing an error.
"""
# sort according to higehst foreground score.
probs, order = roi_probs_neg[:, 1:].max(1)[0].sort(descending=True)
select = torch.tensor((poolsize * int(negative_count), order.size()[0])).min().int()
pool_indices = order[:select]
rand_idx = torch.randperm(pool_indices.size()[0])
return pool_indices[rand_idx[:negative_count].cuda()]
############################################################
# Weight Init
############################################################
def initialize_weights(net):
"""Initialize model weights. Current Default in Pytorch (version 0.4.1) is initialization from a uniform distriubtion.
Will expectably be changed to kaiming_uniform in future versions.
"""
init_type = net.cf.weight_init
for m in [module for module in net.modules() if type(module) in [torch.nn.Conv2d, torch.nn.Conv3d,
torch.nn.ConvTranspose2d,
torch.nn.ConvTranspose3d,
torch.nn.Linear]]:
if init_type == 'xavier_uniform':
torch.nn.init.xavier_uniform_(m.weight.data)
if m.bias is not None:
m.bias.data.zero_()
elif init_type == 'xavier_normal':
torch.nn.init.xavier_normal_(m.weight.data)
if m.bias is not None:
m.bias.data.zero_()
elif init_type == "kaiming_uniform":
torch.nn.init.kaiming_uniform_(m.weight.data, mode='fan_out', nonlinearity=net.cf.relu, a=0)
if m.bias is not None:
fan_in, fan_out = torch.nn.init._calculate_fan_in_and_fan_out(m.weight.data)
bound = 1 / np.sqrt(fan_out)
torch.nn.init.uniform_(m.bias, -bound, bound)
elif init_type == "kaiming_normal":
torch.nn.init.kaiming_normal_(m.weight.data, mode='fan_out', nonlinearity=net.cf.relu, a=0)
if m.bias is not None:
fan_in, fan_out = torch.nn.init._calculate_fan_in_and_fan_out(m.weight.data)
bound = 1 / np.sqrt(fan_out)
torch.nn.init.normal_(m.bias, -bound, bound)
net.logger.info("applied {} weight init.".format(init_type))
\ No newline at end of file