"""
@package medpy.metric.surface
Holds a metrics class computing surface metrics over two 3D-images contain each a binary object.
Classes:
- Surface: Computes different surface metrics between two 3D-images contain each an object.
@author Oskar Maier
@version r0.4.1
@since 2011-12-01
@status Release
"""
# build-in modules
import math
# third-party modules
import scipy.spatial
import scipy.ndimage.morphology
# own modules
# code
class Surface(object):
"""
Computes different surface metrics between two 3D-images contain each an object.
The surface of the objects is computed using a 18-neighbourhood edge detection.
The distance metrics are computed over all points of the surfaces using the nearest
neighbour approach.
Beside this provides a number of statistics of the two images.
During the initialization the edge detection is run for both images, taking up to
5 min (on 512^3 images). The first call to one of the metric measures triggers the
computation of the nearest neighbours, taking up to 7 minutes (based on 250.000 edge
point for each of the objects, which corresponds to a typical liver mask). All
subsequent calls to one of the metrics measures can be expected be in the
sub-millisecond area.
Metrics defined in:
Heimann, T.; van Ginneken, B.; Styner, M.A.; Arzhaeva, Y.; Aurich, V.; Bauer, C.; Beck, A.; Becker, C.; Beichel, R.; Bekes, G.; Bello, F.; Binnig, G.; Bischof, H.; Bornik, A.; Cashman, P.; Ying Chi; Cordova, A.; Dawant, B.M.; Fidrich, M.; Furst, J.D.; Furukawa, D.; Grenacher, L.; Hornegger, J.; Kainmuller, D.; Kitney, R.I.; Kobatake, H.; Lamecker, H.; Lange, T.; Jeongjin Lee; Lennon, B.; Rui Li; Senhu Li; Meinzer, H.-P.; Nemeth, G.; Raicu, D.S.; Rau, A.-M.; van Rikxoort, E.M.; Rousson, M.; Rusko, L.; Saddi, K.A.; Schmidt, G.; Seghers, D.; Shimizu, A.; Slagmolen, P.; Sorantin, E.; Soza, G.; Susomboon, R.; Waite, J.M.; Wimmer, A.; Wolf, I.; , "Comparison and Evaluation of Methods for Liver Segmentation From CT Datasets," Medical Imaging, IEEE Transactions on , vol.28, no.8, pp.1251-1265, Aug. 2009
doi: 10.1109/TMI.2009.2013851
"""
# The edge points of the mask object.
__mask_edge_points = None
# The edge points of the reference object.
__reference_edge_points = None
# The nearest neighbours distances between mask and reference edge points.
__mask_reference_nn = None
# The nearest neighbours distances between reference and mask edge points.
__reference_mask_nn = None
# Distances of the two objects surface points.
__distance_matrix = None
def __init__(self, mask, reference, physical_voxel_spacing = [1,1,1], mask_offset = [0,0,0], reference_offset = [0,0,0]):
"""
Initialize the class with two binary images, each containing a single object.
Assumes the input to be a representation of a 3D image, that fits one of the
following formats:
- 1. all 0 values denoting background, all others the foreground/object
- 2. all False values denoting the background, all others the foreground/object
The first image passed is referred to as 'mask', the second as 'reference'. This
is only important for some metrics that are not symmetric (and therefore not
really metrics).
@param mask binary mask as an scipy array (3D image)
@param reference binary reference as an scipy array (3D image)
@param physical_voxel_spacing The physical voxel spacing of the two images
(must be the same for both)
@param mask_offset offset of the mask array to 0,0,0-origin
@param reference_offset offset of the reference array to 0,0,0-origin
"""
# compute edge images
mask_edge_image = Surface.compute_contour(mask)
reference_edge_image = Surface.compute_contour(reference)
# collect the object edge voxel positions
# !TODO: When the distance matrix is already calculated here
# these points don't have to be actually stored, only their number.
# But there might be some later metric implementation that requires the
# points and then it would be good to have them. What is better?
mask_pts = mask_edge_image.nonzero()
mask_edge_points = zip(mask_pts[0], mask_pts[1], mask_pts[2])
reference_pts = reference_edge_image.nonzero()
reference_edge_points = zip(reference_pts[0], reference_pts[1], reference_pts[2])
# check if there is actually an object present
if 0 >= len(mask_edge_points):
raise Exception('The mask image does not seem to contain an object.')
if 0 >= len(reference_edge_points):
raise Exception('The reference image does not seem to contain an object.')
# add offsets to the voxels positions and multiply with physical voxel spacing
# to get the real positions in millimeters
physical_voxel_spacing = scipy.array(physical_voxel_spacing)
mask_edge_points += scipy.array(mask_offset)
mask_edge_points *= physical_voxel_spacing
reference_edge_points += scipy.array(reference_offset)
reference_edge_points *= physical_voxel_spacing
# set member vars
self.__mask_edge_points = mask_edge_points
self.__reference_edge_points = reference_edge_points
def get_maximum_symmetric_surface_distance(self):
"""
Computes the maximum symmetric surface distance, also known as Hausdorff
distance, between the two objects surfaces.
@return the maximum symmetric surface distance in millimeters
For a perfect segmentation this distance is 0. This metric is sensitive to
outliers and returns the true maximum error.
Metric definition:
Let \f$S(A)\f$ denote the set of surface voxels of \f$A\f$. The shortest
distance of an arbitrary voxel \f$v\f$ to \f$S(A)\f$ is defined as:
\f[
d(v,S(A)) = \min_{s_A\in S(A)} ||v-s_A||
\f]
where \f$||.||\f$ denotes the Euclidean distance. The maximum symmetric
surface distance is then given by:
\f[
MSD(A,B) = \max
\left\{
\max_{s_A\in S(A)} d(s_A,S(B)),
\max_{s_B\in S(B)} d(s_B,S(A)),
\right\}
\f]
"""
# Get the maximum of the nearest neighbour distances
A_B_distance = self.get_mask_reference_nn().max()
B_A_distance = self.get_reference_mask_nn().max()
# compute result and return
return max(A_B_distance, B_A_distance)
def get_root_mean_square_symmetric_surface_distance(self):
"""
Computes the root mean square symmetric surface distance between the
two objects surfaces.
@return root mean square symmetric surface distance in millimeters
For a perfect segmentation this distance is 0. This metric punishes large
deviations from the true contour stronger than the average symmetric surface
distance.
Metric definition:
Let \f$S(A)\f$ denote the set of surface voxels of \f$A\f$. The shortest
distance of an arbitrary voxel \f$v\f$ to \f$S(A)\f$ is defined as:
\f[
d(v,S(A)) = \min_{s_A\in S(A)} ||v-s_A||
\f]
where \f$||.||\f$ denotes the Euclidean distance. The root mean square
symmetric surface distance is then given by:
\f[
RMSD(A,B) =
\sqrt{\frac{1}{|S(A)|+|S(B)|}}
\times
\sqrt{
\sum_{s_A\in S(A)} d^2(s_A,S(B))
+
\sum_{s_B\in S(B)} d^2(s_B,S(A))
}
\f]
"""
# get object sizes
mask_surface_size = len(self.get_mask_edge_points())
reference_surface_sice = len(self.get_reference_edge_points())
# get minimal nearest neighbours distances
A_B_distances = self.get_mask_reference_nn()
B_A_distances = self.get_reference_mask_nn()
# square the distances
A_B_distances_sqrt = A_B_distances * A_B_distances
B_A_distances_sqrt = B_A_distances * B_A_distances
# sum the minimal distances
A_B_distances_sum = A_B_distances_sqrt.sum()
B_A_distances_sum = B_A_distances_sqrt.sum()
# compute result and return
return math.sqrt(1. / (mask_surface_size + reference_surface_sice)) * math.sqrt(A_B_distances_sum + B_A_distances_sum)
def get_average_symmetric_surface_distance(self):
"""
Computes the average symmetric surface distance between the
two objects surfaces.
@return average symmetric surface distance in millimeters
For a perfect segmentation this distance is 0.
Metric definition:
Let \f$S(A)\f$ denote the set of surface voxels of \f$A\f$. The shortest
distance of an arbitrary voxel \f$v\f$ to \f$S(A)\f$ is defined as:
\f[
d(v,S(A)) = \min_{s_A\in S(A)} ||v-s_A||
\f]
where \f$||.||\f$ denotes the Euclidean distance. The average symmetric
surface distance is then given by:
\f[
ASD(A,B) =
\frac{1}{|S(A)|+|S(B)|}
\left(
\sum_{s_A\in S(A)} d(s_A,S(B))
+
\sum_{s_B\in S(B)} d(s_B,S(A))
\right)
\f]
"""
# get object sizes
mask_surface_size = len(self.get_mask_edge_points())
reference_surface_sice = len(self.get_reference_edge_points())
# get minimal nearest neighbours distances
A_B_distances = self.get_mask_reference_nn()
B_A_distances = self.get_reference_mask_nn()
# sum the minimal distances
A_B_distances = A_B_distances.sum()
B_A_distances = B_A_distances.sum()
# compute result and return
return 1. / (mask_surface_size + reference_surface_sice) * (A_B_distances + B_A_distances)
def get_mask_reference_nn(self):
"""
@return The distances of the nearest neighbours of all mask edge points to all
reference edge points.
"""
# Note: see note for @see get_reference_mask_nn
if None == self.__mask_reference_nn:
tree = scipy.spatial.cKDTree(self.get_mask_edge_points())
self.__mask_reference_nn, _ = tree.query(self.get_reference_edge_points())
return self.__mask_reference_nn
def get_reference_mask_nn(self):
"""
@return The distances of the nearest neighbours of all reference edge points
to all mask edge points.
The underlying algorithm used for the scipy.spatial.KDTree implementation is
based on:
Sunil Arya, David M. Mount, Nathan S. Netanyahu, Ruth Silverman, and
Angela Y. Wu. 1998. An optimal algorithm for approximate nearest neighbor
searching fixed dimensions. J. ACM 45, 6 (November 1998), 891-923
"""
# Note: KDTree is faster than scipy.spatial.distance.cdist when the number of
# voxels exceeds 10.000 (computationally tested). The maximum complexity is
# O(D*N^2) vs. O(D*N*log(N), where D=3 and N=number of voxels
if None == self.__reference_mask_nn:
tree = scipy.spatial.cKDTree(self.get_reference_edge_points())
self.__reference_mask_nn, _ = tree.query(self.get_mask_edge_points())
return self.__reference_mask_nn
def get_mask_edge_points(self):
"""
@return The edge points of the mask object.
"""
return self.__mask_edge_points
def get_reference_edge_points(self):
"""
@return The edge points of the reference object.
"""
return self.__reference_edge_points
@staticmethod
def compute_contour(array):
"""
Uses a 18-neighbourhood filter to create an edge image of the input object.
Assumes the input to be a representation of a 3D image, that fits one of the
following formats:
- 1. all 0 values denoting background, all others the foreground/object
- 2. all False values denoting the background, all others the foreground/object
The area outside the array is assumed to contain background voxels. The method
does not ensure that the object voxels are actually connected, this is silently
assumed.
@param array a numpy array with only 0/N\{0} or False/True values.
@return a boolean numpy array with the input objects edges
"""
# set 18-neighbourhood/conectivity (for 3D images) alias face-and-edge kernel
# all values covered by 1/True passed to the function
# as a 1D array in order left-right, top-down
# Note: all in all 19 ones, as the center value
# also has to be checked (if it is a masked pixel)
# [[[0, 1, 0], [[1, 1, 1], [[0, 1, 0],
# [1, 1, 1], [1, 1, 1], [1, 1, 1],
# [0, 1, 0]], [1, 1, 1]], [0, 1, 0]]]
footprint = scipy.ndimage.morphology.generate_binary_structure(3, 2)
# create an erode version of the array
erode_array = scipy.ndimage.morphology.binary_erosion(array, footprint)
# xor the erode_array with the original and return
return array ^ erode_array
验证4
'''
Contains common functions for reading data out of leveldb
@author: Mohamed.Ezz
'''
import plyvel, lmdb
import numpy as np
from caffe.proto import caffe_pb2
IMG_DTYPE = np.float
SEG_DTYPE = np.uint8
def denormalize_img_255(arr):
""" Denormalizes a nparray to 0-255 values """
min = arr.min()
max = arr.max()
new = (arr - min) * (255.0 / (max - min))
return new.astype(np.uint8)
def leveldb_arrays(leveldbdir):
""" Generator. Given leveldb directory, iterate the stored data as numpy arrays. Yields (Key, NumpyArray) """
db = CaffeDatabase(leveldbdir)
for k, v in db.iterator():
yield k, to_numpy_matrix(v)
def nth_datum(caffedb, n):
""" Returns nth datum. 0-based index"""
n += 1
it = caffedb.iterator()
for _ in range(n):
_, v = it.next()
datum = caffe_pb2.Datum()
datum.ParseFromString(v)
return datum
def get_data_type(datum):
""" By simple calculations, conclude the size of integers stored in datum.data """
n_values = datum.height * datum.width * datum.channels
n_bytes = len(datum.data)
int_size = float(n_bytes) / n_values
if int_size != int(int_size) or int_size not in [1, 2, 4, 8]:
raise ValueError("Can't find int size. n_values : %i , n_bytes : %i" % (n_values, n_bytes))
types = {1: np.int8, 2: np.int16, 4: np.int32, 8: np.int64}
type_ = types[int(int_size)]
return type_
def find_keycount(caffedb, count_values=None):
""" Takes a CaffeDatabase or plyvel.DB instance and returns number of keys found and count of each value.
count_values is a list of values to count, e.g. count_values=[0,1,2] will return [count of 1s, count of 2s, count of 3s]
if count_values is None, return value of this function is [],key_count"""
count = 0
total_value_counts = np.array([0] * len(count_values or []))
for _, v in caffedb.iterator():
count += 1
if count_values is not None:
array = to_numpy_matrix(v)
current_count = np.array([0] * len(count_values))
for i, val in enumerate(count_values):
current_count[i] = np.sum(array == val)
total_value_counts += current_count
return total_value_counts, count
def to_numpy_matrix(v):
""" Convert leveldb/lmdb value to numpy matrix of shape N x N """
datum = caffe_pb2.Datum()
datum.ParseFromString(v)
# Three cases
# 1- int imgs in data,
# 2- int8 labels in data
if len(datum.data) > 0:
type_ = get_data_type(datum)
matrix = np.fromstring(datum.data, dtype=type_)
# 3- float imgs in float_data,
elif len(datum.float_data) > 0:
matrix = np.array(datum.float_data)
else:
raise ValueError("Serialized datum have empty data and float_data.")
matrix = matrix.reshape((datum.height, datum.width))
return matrix
def norm_hounsfield_dyn(arr, c_min=0.1, c_max=0.3):
""" Converts from hounsfield units to float64 image with range 0.0 to 1.0 """
# calc min and max
min, max = np.amin(arr), np.amax(arr)
arr = arr.astype(IMG_DTYPE)
if min <= 0:
arr = np.clip(arr, min * c_min, max * c_max)
# right shift to zero
arr = np.abs(min * c_min) + arr
else:
arr = np.clip(arr, min, max * c_max)
# left shift to zero
arr = arr - min
# normalization
norm_fac = np.amax(arr)
if norm_fac != 0:
norm = np.divide(
np.multiply(arr, 255),
np.amax(arr))
else: # don't divide through 0
norm = np.multiply(arr, 255)
norm = np.clip(np.multiply(norm, 0.00390625), 0, 1)
return norm
def norm_hounsfield_stat(arr, c_min=-100, c_max=200):
min = np.amin(arr)
arr = np.array(arr, dtype=IMG_DTYPE)
if min <= 0:
# clip
c_arr = np.clip(arr, c_min, c_max)
# right shift to zero
slc_0 = np.add(np.abs(min), c_arr)
else:
# clip
c_arr = np.clip(arr, c_min, c_max)
# left shift to zero
slc_0 = np.subtract(c_arr, min)
# normalization
norm_fac = np.amax(slc_0)
if norm_fac != 0:
norm = np.divide(
np.multiply(
slc_0,
255
),
np.amax(slc_0)
)
else: # don't divide through 0
norm = np.multiply(slc_0, 255)
norm = np.clip(np.multiply(norm, 0.00390625), 0, 1)
return norm
class CaffeDatabase():
""" Abstraction layer over lmdb and leveldb """
def __init__(self, path, backend='lmdb'):
self.backend = backend
assert backend in ['lmdb', 'leveldb'], "Database backend not known :%s" % backend
if backend == 'lmdb':
self.db = lmdb.open(path)
elif backend == 'leveldb':
self.db = plyvel.DB(path)
def iterator(self):
if self.backend == 'lmdb':
txn = self.db.begin()
cursor = txn.cursor()
it = cursor.iternext()
elif self.backend == 'leveldb':
it = self.db.iterator()
return it
验证5
import config
import logging
import scipy as sp
import scipy.misc, scipy.ndimage.interpolation
import caffe
caffe.set_mode_gpu()
import matplotlib
from matplotlib import pyplot as plt
matplotlib.use('Agg')
import os
from denseinference import CRFProcessor
from medpy import metric
import nibabel as nib
import numpy as np
import IPython
# this should actually be part of medpy. Apparently it isn't (anymore). So the surface.py file from http://pydoc.net/Python/MedPy/0.2.2/medpy.metric._surface/ should be manually imported
from surface import Surface
from utils import norm_hounsfield_stat, norm_hounsfield_dyn
IMG_DTYPE = np.float
SEG_DTYPE = np.uint8
def miccaiimshow(img, seg, preds, fname, titles=None, plot_separate_img=True):
"""Takes raw image img, seg in range 0-2, list of predictions in range 0-2"""
plt.figure(figsize=(25, 25))
ALPHA = 1
n_plots = len(preds)
subplot_offset = 0
plt.set_cmap('gray')
if plot_separate_img:
n_plots += 1
subplot_offset = 1
plt.subplot(1, n_plots, 1)
plt.subplots_adjust(wspace=0, hspace=0)
plt.title("Image")
plt.axis('off')
plt.imshow(img, cmap="gray")
if type(preds) != list:
preds = [preds]
for i, pred in enumerate(preds):
# Order of overaly
########## OLD
# lesion= pred==2
# difflesion = set_minus(seg==2,lesion)
# liver = set_minus(pred==1, [lesion, difflesion])
# diffliver = set_minus(seg==1, [liver,lesion,difflesion])
##########
lesion = pred == 2
difflesion = np.logical_xor(seg == 2, lesion)
liver = pred == 1
diffliver = np.logical_xor(seg == 1, liver)
plt.subplot(1, n_plots, i + 1 + subplot_offset)
title = titles[i] if titles is not None and i < len(titles) else ""
plt.title(title)
plt.axis('off')
plt.imshow(img);
plt.hold(True)
# Liver prediction
plt.imshow(np.ma.masked_where(liver == 0, liver), cmap="Greens", vmin=0.1, vmax=1.2, alpha=ALPHA);
plt.hold(True)
# Liver : Pixels in ground truth, not in prediction
plt.imshow(np.ma.masked_where(diffliver == 0, diffliver), cmap="Spectral", vmin=0.1, vmax=2.2, alpha=ALPHA);
plt.hold(True)
# Lesion prediction
plt.imshow(np.ma.masked_where(lesion == 0, lesion), cmap="Blues", vmin=0.1, vmax=1.2, alpha=ALPHA);
plt.hold(True)
# Lesion : Pixels in ground truth, not in prediction
plt.imshow(np.ma.masked_where(difflesion == 0, difflesion), cmap="Reds", vmin=0.1, vmax=1.5, alpha=ALPHA)
plt.savefig(fname)
plt.close()
def to_scale(img, shape=None):
if shape is None:
shape = config.slice_shape
height, width = shape
if img.dtype == SEG_DTYPE:
return scipy.misc.imresize(img, (height, width), interp="nearest").astype(SEG_DTYPE)
elif img.dtype == IMG_DTYPE:
max_ = np.max(img)
factor = 256.0 / max_ if max_ != 0 else 1
return (scipy.misc.imresize(img, (height, width), interp="nearest") / factor).astype(IMG_DTYPE)
else:
raise TypeError(
'Error. To scale the image array, its type must be np.uint8 or np.float64. (' + str(img.dtype) + ')')
def histeq_processor(img):
"""Histogram equalization"""
nbr_bins = 256
# get image histogram
imhist, bins = np.histogram(img.flatten(), nbr_bins, normed=True)
cdf = imhist.cumsum() # cumulative distribution function
cdf = 255 * cdf / cdf[-1] # normalize
# use linear interpolation of cdf to find new pixel values
original_shape = img.shape
img = np.interp(img.flatten(), bins[:-1], cdf)
img = img / 255.0
return img.reshape(original_shape)
def downscale_img_label(imgvol, label_vol):
"""
Downscales an image volume and an label volume. Normalizes the hounsfield units of the image volume
:param imgvol:
:param label_vol:
:return:
"""
imgvol = imgvol.astype(IMG_DTYPE)
label_vol = label_vol.astype(SEG_DTYPE)
slc=None
imgvol_downscaled = np.zeros((config.slice_shape[0], config.slice_shape[1], imgvol.shape[2]))
label_vol_downscaled = np.zeros((config.slice_shape[0], config.slice_shape[1], imgvol.shape[2]))
# Copy image volume
# copy_imgvol = np.copy(imgvol)
# Truncate metal and high absorbative objects
logging.info('Found' + str(np.sum(imgvol > 1200)) + 'values > 1200 !!')
imgvol[imgvol > 1200] = 0
for i in range(imgvol.shape[2]):
# Get the current slc, normalize and downscale
slc = imgvol[:, :, i]
if config.ct_window_type == 'dyn':
slc = norm_hounsfield_dyn(slc, c_min=config.ct_window_type_min, c_max=config.ct_window_type_max)
elif config.ct_window_type == 'stat':
slc = norm_hounsfield_stat(slc, c_min=config.ct_window_type_min, c_max=config.ct_window_type_max)
else:
print "CT Windowing did not work."
slc = to_scale(slc, config.slice_shape)
# slc = histeq_processor(slc)
imgvol_downscaled[:, :, i] = slc
# downscale the label slc for the crf
label_vol_downscaled[:, :, i] = to_scale(label_vol[:, :, i], config.slice_shape)
return [imgvol_downscaled, label_vol_downscaled]
def scorer(pred, label):
"""
:param pred:
:param label:
:param voxelspacing:
:return:
"""
volscores = {}
volscores['dice'] = metric.dc(pred, label)
volscores['jaccard'] = metric.binary.jc(pred, label)
volscores['voe'] = 1. - volscores['jaccard']
volscores['rvd'] = metric.ravd(label, pred)
if np.count_nonzero(pred) == 0 or np.count_nonzero(label) == 0:
volscores['assd'] = 0
volscores['msd'] = 0
# else:
# evalsurf = Surface(pred, label, physical_voxel_spacing=vxlspacing, mask_offset=[0., 0., 0.],
# reference_offset=[0., 0., 0.])
# volscores['assd'] = evalsurf.get_average_symmetric_surface_distance()
#
# volscores['msd'] = metric.hd(label, pred, voxelspacing=vxlspacing)
logging.info("\tDice " + str(volscores['dice']))
logging.info("\tJaccard " + str(volscores['jaccard']))
logging.info("\tVOE " + str(volscores['voe']))
logging.info("\tRVD " + str(volscores['rvd']))
# logging.info("\tASSD " + str(volscores['assd']))
# logging.info("\tMSD " + str(volscores['msd']))
return volscores
def get_average_score(scorelist, scorename, mode=None):
"""
:param scorelist:
:param scorename:
:return:
"""
score = 0.
for e in scorelist:
if mode == 'abs':
score += np.abs(e[scorename])
else:
score += e[scorename]
score /= float(len(scorelist))
return score
def zoomliver_UNET_processor(img, seg):
""" Custom preprocessing of img,seg for UNET architecture:
Crops the background and upsamples the found patch."""
# Remove background !
img = np.multiply(img, np.clip(seg, 0, 1))
# get patch size
col_maxes = np.max(seg, axis=0) # a row
row_maxes = np.max(seg, axis=1) # a column
nonzero_colmaxes = np.nonzero(col_maxes)[0]
nonzero_rowmaxes = np.nonzero(row_maxes)[0]
x1, x2 = nonzero_colmaxes[0], nonzero_colmaxes[-1]
y1, y2 = nonzero_rowmaxes[0], nonzero_rowmaxes[-1]
width = x2 - x1
height = y2 - y1
MIN_WIDTH = 60
MIN_HEIGHT = 60
x_pad = (MIN_WIDTH - width) / 2 if width < MIN_WIDTH else 0
y_pad = (MIN_HEIGHT - height) / 2 if height < MIN_HEIGHT else 0
x1 = max(0, x1 - x_pad)
x2 = min(img.shape[1], x2 + x_pad)
y1 = max(0, y1 - y_pad)
y2 = min(img.shape[0], y2 + y_pad)
img = img[y1:y2 + 1, x1:x2 + 1]
seg = seg[y1:y2 + 1, x1:x2 + 1]
img = to_scale(img, (388, 388))
seg = to_scale(seg, (388, 388))
# All non-lesion is background
seg[seg == 1] = 0
# Lesion label becomes 1
seg[seg == 2] = 1
# Now do padding for UNET, which takes 572x572
# seg=np.pad(seg,((92,92),(92,92)),mode='reflect')
img = np.pad(img, 92, mode='reflect')
return img, (x1, x2, y1, y2)
if __name__ == '__main__':
try:
logging.basicConfig(filename=os.path.join(config.output_dir, config.logfile), filemode='w',
level=config.log_level, format='%(asctime)s %(levelname)s:%(message)s',
datefmt='%d-%m-%Y %I:%M:%S %p')
# lists to calculate the overall score over all folds from, i.e. holds scores of all volumes
overall_score_liver = []
overall_score_lesion_crf = []
overall_score_liver_crf = []
overall_score_lesion = []
# Iterate folds and corresponding models
for fold, model, deployprototxt, model_step_two, deployprototxt_step_two in zip(config.dataset, config.models,
config.deployprototxt,
config.models_step_two,
config.deployprototxt_step_two):
logging.info("Starting new fold")
# Lists to save scores for each volume of this fold
foldscore_lesion_crf = []
foldscore_liver_crf = []
foldscore_liver = []
foldscore_lesion = []
# Iterate volumes in fold
for volidx, volpaths in enumerate(fold):
logging.info("Loading Network for Step 1")
# load new network for this fold
try:
del net # it is a good idea to delete the net object to free up memory before instantiating another one
net = caffe.Net(deployprototxt, model, caffe.TEST)
except NameError:
net = caffe.Net(deployprototxt, model, caffe.TEST)
logging.info("Loading " + volpaths[1])
imgvol = np.load(volpaths[1])
labelvol = np.load(volpaths[2])
# the raw probabilites of step 1
probvol = np.zeros((config.slice_shape[0], config.slice_shape[1], imgvol.shape[2], 2))
# the probabilites of step 2 scaled back down into the volume
pred_step_two = np.zeros((config.slice_shape[0], config.slice_shape[1], imgvol.shape[2]))
pred_step_one = np.zeros((config.slice_shape[0], config.slice_shape[1], imgvol.shape[2]))
probvol_step_two = np.zeros((config.slice_shape[0], config.slice_shape[1], imgvol.shape[2], 2))
# rotate volumes so that the networks sees them in the same orientation like during training
imgvol = np.rot90(imgvol)
labelvol = np.rot90(labelvol)
imgvol_downscaled, labelvol_downscaled = downscale_img_label(imgvol, labelvol)
# iterate slices in volume and do prediction
logging.info("Predicting " + volpaths[1])
for i in range(imgvol_downscaled.shape[2]):
slc = imgvol_downscaled[:, :, i]
# create mirrored slc for unet
slc = np.pad(slc, ((92, 92), (92, 92)), mode='reflect')
# load slc into network and do forward pass
net.blobs['data'].data[...] = slc
net.forward()
# now save raw probabilities
probvol[:, :, i, :] = net.blobs['prob'].data.transpose((0, 2, 3, 1))[0]
pred_step_one[:, :, i] = np.argmax(probvol[:, :, i, :], axis=2)
# result shape is batch_img_idx , height, width, probability_of_class
# dump probabiliteis to .npy file for future use
# np.save('./probfiles/' + ))
##FIX THIS
logging.info("Here are the liver scores before CRF:")
# calculate scores for liver
pred_to_use = np.logical_or(probvol.argmax(3) == 1, probvol.argmax(3) == 2)
label_to_use = np.logical_or(labelvol_downscaled == 1, labelvol_downscaled == 2)
#voxelspacing = volpaths[3]
volumescore_liver = scorer(pred_to_use, label_to_use)
# Run Liver CRF
logging.info("Now running CRF on Liver")
crfparams = {'max_iterations': 10, 'dynamic_z': True, 'ignore_memory': True, 'pos_x_std': 1.5,
'pos_y_std': 1.5,
'pos_z_std': 1.5, 'pos_w': 3.0, 'bilateral_x_std': 9.0, 'bilateral_y_std': 9.0,
'bilateral_z_std': 9.0, 'bilateral_intensity_std': 20.0, 'bilateral_w': 10.0}
pro = CRFProcessor.CRF3DProcessor(**crfparams)
if config.save_probability_volumes:
np.save(os.path.join(config.output_dir, os.path.basename(volpaths[1])) + ".liver.npy", probvol)
crf_pred_liver = pro.set_data_and_run(imgvol_downscaled, probvol)
# calculate scores for liver
label_to_use = np.logical_or(labelvol_downscaled == 1, labelvol_downscaled == 2)
logging.info("Here are the liver scores after CRF:")
volumescore_liver_crf = scorer(crf_pred_liver, label_to_use)
# calculate scores for lesions
# pred_to_use = probvol.argmax(3)==2
# label_to_use = labelvol_downscaled==2
# volumescore_lesion = scorer(pred_to_use,label_to_use,voxelspacing)
# OK, we're done on the first step of the cascaded networks and have evaluated them.
# Now let's get to the second step.
del net
logging.info("Deleted network for cascade step 1")
net = caffe.Net(deployprototxt_step_two, model_step_two, caffe.TEST)
logging.info("Loaded network for cascade step 2")
# we again iterate over all slices in the volume
for i in range(imgvol_downscaled.shape[2]):
slc = imgvol_downscaled[:, :, i]
# create mirrored slc for unet
# slc = np.pad(slc,((92,92),(92,92)),mode='reflect')
# now we crop and upscale the liver
slc_crf_pred_liver = crf_pred_liver[:, :, i].astype(SEG_DTYPE)
# slc_crf_pred_liver = pred_to_use[:,:,i].astype(SEG_DTYPE)
# slc_crf_pred_liver = labelvol_downscaled[:,:,i]
if np.count_nonzero(slc_crf_pred_liver) == 0:
probvol_step_two[:, :, i, :] = 0
else:
slc, bbox = zoomliver_UNET_processor(slc, slc_crf_pred_liver)
# load slc into network and do forward pass
net.blobs['data'].data[...] = slc
net.forward()
# scale output back down and insert into the probability volume
x1, x2, y1, y2 = bbox
leftpad, rightpad = x1, 388 - x2
toppad, bottompad = y1, 388 - y2
width, height = int(x2 - x1), int(y2 - y1)
# now save probabilities
prob = net.blobs['prob'].data.transpose((0, 2, 3, 1))[0]
# probvol[:,:,i,:] = prob
slc_pred_step_two = np.argmax(prob, axis=2).astype(SEG_DTYPE)
slc_pred_step_two = to_scale(slc_pred_step_two, (height, width))
slc_pred_step_two = np.pad(slc_pred_step_two, ((toppad, bottompad), (leftpad, rightpad)),
mode='constant')
pred_step_two[:, :, i] = slc_pred_step_two
prob0 = prob[:, :, 0].astype(
IMG_DTYPE) # use IMG_DTYPE bcoz we've probabiblities, not hard labels
prob0 = to_scale(prob0, (height, width))
prob0 = np.pad(prob0, ((toppad, bottompad), (leftpad, rightpad)), mode='constant')
#
#
prob1 = prob[:, :, 1].astype(IMG_DTYPE)
prob1 = to_scale(prob1, (height, width))
prob1 = np.pad(prob1, ((toppad, bottompad), (leftpad, rightpad)), mode='constant')
probvol_step_two[:, :, i, 0] = prob0
probvol_step_two[:, :, i, 1] = prob1
# probvol_step_two[bbox[0]:bbox[0] + bbox[1], bbox[2]:bbox[2] + bbox[3], i, :] =
logging.info("Lesion scores after step 2 before CRF")
# pred_to_use = probvol_step_two.argmax(3) == 2
pred_to_use = pred_step_two.astype(SEG_DTYPE)
label_to_use = labelvol_downscaled == 2
volumescore_lesion = scorer(pred_to_use, label_to_use)
# Save lesion npy probabilities
if config.save_probability_volumes:
np.save(os.path.join(config.output_dir, os.path.basename(volpaths[1])) + ".lesion.npy",
probvol_step_two)
### SAVE PLOTS
if config.plot_every_n_slices > 0:
for i in range(0, imgvol_downscaled.shape[2], config.plot_every_n_slices):
pred_vol_bothsteps = pred_step_one
pred_vol_bothsteps[pred_step_two == 1] = 2
liverdc = metric.dc(pred_step_one[:, :, i], labelvol_downscaled[:, :, i] == 1)
lesiondc = metric.dc(pred_step_two[:, :, i], labelvol_downscaled[:, :, i] == 2)
fname = os.path.join(config.output_dir, os.path.basename(volpaths[1]))
fname += "_slc" + str(i) + "_"
fname += "liv" + str(liverdc) + "_les" + str(lesiondc) + ".png"
# logging.info("Plotting "+fname)
miccaiimshow(imgvol_downscaled[:, :, i], labelvol_downscaled[:, :, i],
[labelvol_downscaled[:, :, i], pred_vol_bothsteps[:, :, i]], fname=fname,
titles=["Ground Truth", "Prediction"], plot_separate_img=True)
logging.info("Now running LESION CRF on Liver")
crf_params = {'ignore_memory': True, 'bilateral_intensity_std': 0.16982742320252908,
'bilateral_w': 6.406401876489639,
'pos_w': 2.3422381267344132, 'bilateral_x_std': 284.5377968491542,
'pos_x_std': 23.636281254341867,
'max_iterations': 10}
pro = CRFProcessor.CRF3DProcessor(**crf_params)
crf_pred_lesion = pro.set_data_and_run(imgvol_downscaled, probvol_step_two)
volumescore_lesion_crf = scorer(crf_pred_lesion, label_to_use)
# Append to results lists so that the average scores can be calculated later
foldscore_liver.append(volumescore_liver)
foldscore_lesion.append(volumescore_lesion)
foldscore_liver_crf.append(volumescore_liver_crf)
foldscore_lesion_crf.append(volumescore_lesion_crf)
overall_score_liver_crf.append(volumescore_liver_crf)
overall_score_lesion_crf.append(volumescore_lesion_crf)
overall_score_liver.append(volumescore_liver)
overall_score_lesion.append(volumescore_lesion)
logging.info("=========================================")
logging.info("Average Liver Scores before CRF for this fold: ")
logging.info("Dice " + str(get_average_score(foldscore_liver, 'dice')))
logging.info("Jaccard " + str(get_average_score(foldscore_liver, 'jaccard')))
logging.info("VOE " + str(get_average_score(foldscore_liver, 'voe')))
logging.info("RVD " + str(get_average_score(foldscore_liver, 'rvd')))
# logging.info("ASSD " + str(get_average_score(foldscore_liver, 'assd')))
#logging.info("MSD " + str(get_average_score(foldscore_liver, 'msd')))
logging.info("=========================================")
logging.info("=========================================")
logging.info("Average Liver Scores after CRF for this fold: ")
logging.info("Dice " + str(get_average_score(foldscore_liver_crf, 'dice')))
logging.info("Jaccard " + str(get_average_score(foldscore_liver_crf, 'jaccard')))
logging.info("VOE " + str(get_average_score(foldscore_liver_crf, 'voe')))
logging.info("RVD " + str(get_average_score(foldscore_liver_crf, 'rvd')))
# logging.info("ASSD " + str(get_average_score(foldscore_liver_crf, 'assd')))
# logging.info("MSD " + str(get_average_score(foldscore_liver_crf, 'msd')))
logging.info("=========================================")
logging.info("=========================================")
logging.info("Average Lesion Scores before CRF for this fold: ")
logging.info("Dice " + str(get_average_score(foldscore_lesion, 'dice')))
logging.info("Jaccard " + str(get_average_score(foldscore_lesion, 'jaccard')))
logging.info("VOE " + str(get_average_score(foldscore_lesion, 'voe')))
logging.info("RVD " + str(get_average_score(foldscore_lesion, 'rvd')))
#logging.info("ASSD " + str(get_average_score(foldscore_lesion, 'assd')))
#logging.info("MSD " + str(get_average_score(foldscore_lesion, 'msd')))
logging.info("=========================================")
logging.info("=========================================")
logging.info("Average Lesion Scores AFTER CRF for this fold: ")
logging.info("Dice " + str(get_average_score(foldscore_lesion_crf, 'dice')))
logging.info("Jaccard " + str(get_average_score(foldscore_lesion_crf, 'jaccard')))
logging.info("VOE " + str(get_average_score(foldscore_lesion_crf, 'voe')))
logging.info("RVD " + str(get_average_score(foldscore_lesion_crf, 'rvd')))
#logging.info("ASSD " + str(get_average_score(foldscore_lesion_crf, 'assd')))
# logging.info("MSD " + str(get_average_score(foldscore_lesion_crf, 'msd')))
logging.info("=========================================")
logging.info("=========================================")
logging.info("DONE WITH PROCESSING ALL FOLDS. NOW THE OVERALL RESULTS COME")
logging.info("=========================================")
logging.info("Average Liver Scores before CRF overall: ")
logging.info("Dice " + str(get_average_score(overall_score_liver, 'dice')))
logging.info("Jaccard " + str(get_average_score(overall_score_liver, 'jaccard')))
logging.info("VOE " + str(get_average_score(overall_score_liver, 'voe')))
logging.info("RVD " + str(get_average_score(overall_score_liver, 'rvd', mode='abs')))
logging.info("ASSD " + str(get_average_score(overall_score_liver, 'assd')))
# logging.info("MSD " + str(get_average_score(overall_score_liver, 'msd')))
# logging.info("=========================================")
logging.info("=========================================")
logging.info("Average Liver Scores after CRF overall: ")
logging.info("Dice " + str(get_average_score(overall_score_liver_crf, 'dice')))
logging.info("Jaccard " + str(get_average_score(overall_score_liver_crf, 'jaccard')))
logging.info("VOE " + str(get_average_score(overall_score_liver_crf, 'voe')))
logging.info("RVD " + str(get_average_score(overall_score_liver_crf, 'rvd', mode='abs')))
#logging.info("ASSD " + str(get_average_score(overall_score_liver_crf, 'assd')))
#logging.info("MSD " + str(get_average_score(overall_score_liver_crf, 'msd')))
logging.info("=========================================")
logging.info("=========================================")
logging.info("Average Lesion Scores before step2 CRF overall: ")
logging.info("Dice " + str(get_average_score(overall_score_lesion, 'dice')))
logging.info("Jaccard " + str(get_average_score(overall_score_lesion, 'jaccard')))
logging.info("VOE " + str(get_average_score(overall_score_lesion, 'voe')))
logging.info("RVD " + str(get_average_score(overall_score_lesion, 'rvd', mode='abs')))
#logging.info("ASSD " + str(get_average_score(overall_score_lesion, 'assd')))
#logging.info("MSD " + str(get_average_score(overall_score_lesion, 'msd')))
logging.info("=========================================")
logging.info("=========================================")
logging.info("Average Lesion Scores after step2 CRF overall: ")
logging.info("Dice " + str(get_average_score(overall_score_lesion_crf, 'dice')))
logging.info("Jaccard " + str(get_average_score(overall_score_lesion_crf, 'jaccard')))
logging.info("VOE " + str(get_average_score(overall_score_lesion_crf, 'voe')))
logging.info("RVD " + str(get_average_score(overall_score_lesion_crf, 'rvd', mode='abs')))
#logging.info("ASSD " + str(get_average_score(overall_score_lesion_crf, 'assd')))
# logging.info("MSD " + str(get_average_score(overall_score_lesion_crf, 'msd')))
logging.info("=========================================")
# Creating CSV
csvarray = np.zeros((len(overall_score_liver), 13))
csvarray[:, 0] = range(1, len(overall_score_liver) + 1)
# csvarray[:,1] = [s['dice'] for s in overall_score_liver]
d = overall_score_liver.iteritems()
for i in range(6):
d.next()[1]
csvarray[:, i + 1] = d.next()[1]
d = overall_score_lesion.iteritems()
for i in range(6):
d.next()[1]
csvarray[:, i + 1 + 6] = d.next()[1]
np.savetxt("Numbers.csv", csvarray, delimiter=",")
except:
logging.exception("Exception happend...")
IPython.embed()
本文详细介绍了参与LITS挑战的肝分割过程，包括配置设置、使用CRF进行后处理优化，以及多个阶段的实验数据和验证步骤，探讨了numpy_data_layer在数据处理中的应用。
扫一扫
603
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
