diff --git a/docs/Makefile b/docs/Makefile
new file mode 100644
index 0000000..4ce41e6
--- /dev/null
+++ b/docs/Makefile
@@ -0,0 +1,20 @@
+# Minimal makefile for Sphinx documentation
+#
+
+# You can set these variables from the command line.
+SPHINXOPTS =
+SPHINXBUILD = sphinx-build
+SPHINXPROJ = pylidc
+SOURCEDIR = source
+BUILDDIR = build
+
+# Put it first so that "make" without argument is like "make help".
+help:
+ @$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
+
+.PHONY: help Makefile
+
+# Catch-all target: route all unknown targets to Sphinx using the new
+# "make mode" option. $(O) is meant as a shortcut for $(SPHINXOPTS).
+%: Makefile
+ @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
\ No newline at end of file
diff --git a/docs/index.html b/docs/index.html
new file mode 100644
index 0000000..aa2b848
--- /dev/null
+++ b/docs/index.html
@@ -0,0 +1,6 @@
+
+
+
+
+
+
diff --git a/docs/source/_static/copybutton.js b/docs/source/_static/copybutton.js
new file mode 100644
index 0000000..a2b21d1
--- /dev/null
+++ b/docs/source/_static/copybutton.js
@@ -0,0 +1,63 @@
+$(document).ready(function() {
+ /* Add a [>>>] button on the top-right corner of code samples to hide
+ * the >>> and ... prompts and the output and thus make the code
+ * copyable. */
+ var div = $('.highlight-python .highlight,' +
+ '.highlight-python3 .highlight')
+ var div = $('.highlight-default > .highlight')
+ var pre = div.find('pre');
+
+ // get the styles from the current theme
+ pre.parent().parent().css('position', 'relative');
+ var hide_text = 'Hide the prompts and output';
+ var show_text = 'Show the prompts and output';
+ var border_width = pre.css('border-top-width');
+ var border_style = pre.css('border-top-style');
+ var border_color = pre.css('border-top-color');
+ var button_styles = {
+ 'cursor':'pointer', 'position': 'absolute', 'top': '0', 'right': '0',
+ 'border-color': border_color, 'border-style': border_style,
+ 'border-width': border_width, 'color': border_color, 'text-size': '75%',
+ 'font-family': 'monospace', 'padding-left': '0.2em', 'padding-right': '0.2em',
+ 'border-radius': '0 3px 0 0'
+ }
+
+ // create and add the button to all the code blocks that contain >>>
+ div.each(function(index) {
+ console.log('blah')
+ var jthis = $(this);
+ if (jthis.find('.gp').length > 0) {
+ var button = $('>>>');
+ button.css(button_styles)
+ button.attr('title', hide_text);
+ button.data('hidden', 'false');
+ jthis.prepend(button);
+ }
+ // tracebacks (.gt) contain bare text elements that need to be
+ // wrapped in a span to work with .nextUntil() (see later)
+ jthis.find('pre:has(.gt)').contents().filter(function() {
+ return ((this.nodeType == 3) && (this.data.trim().length > 0));
+ }).wrap('');
+ });
+
+ // define the behavior of the button when it's clicked
+ $('.copybutton').click(function(e){
+ e.preventDefault();
+ var button = $(this);
+ if (button.data('hidden') === 'false') {
+ // hide the code output
+ button.parent().find('.go, .gp, .gt').hide();
+ button.next('pre').find('.gt').nextUntil('.gp, .go').css('visibility', 'hidden');
+ button.css('text-decoration', 'line-through');
+ button.attr('title', show_text);
+ button.data('hidden', 'true');
+ } else {
+ // show the code output
+ button.parent().find('.go, .gp, .gt').show();
+ button.next('pre').find('.gt').nextUntil('.gp, .go').css('visibility', 'visible');
+ button.css('text-decoration', 'none');
+ button.attr('title', hide_text);
+ button.data('hidden', 'false');
+ }
+ });
+});
diff --git a/docs/source/annotation.rst b/docs/source/annotation.rst
new file mode 100644
index 0000000..33cb4c8
--- /dev/null
+++ b/docs/source/annotation.rst
@@ -0,0 +1,5 @@
+The Annotation Model
+====================
+
+.. autoclass:: pylidc.Annotation
+ :members:
diff --git a/docs/source/conf.py b/docs/source/conf.py
new file mode 100644
index 0000000..11f6049
--- /dev/null
+++ b/docs/source/conf.py
@@ -0,0 +1,180 @@
+# -*- coding: utf-8 -*-
+#
+# pylidc documentation build configuration file, created by
+# sphinx-quickstart on Mon Mar 26 15:44:42 2018.
+#
+# This file is execfile()d with the current directory set to its
+# containing dir.
+#
+# Note that not all possible configuration values are present in this
+# autogenerated file.
+#
+# All configuration values have a default; values that are commented out
+# serve to show the default.
+
+# If extensions (or modules to document with autodoc) are in another directory,
+# add these directories to sys.path here. If the directory is relative to the
+# documentation root, use os.path.abspath to make it absolute, like shown here.
+
+import os
+import sys
+sys.path.insert(0, os.path.abspath('../pylidc'))
+
+
+# -- General configuration ------------------------------------------------
+
+# If your documentation needs a minimal Sphinx version, state it here.
+#
+# needs_sphinx = '1.0'
+
+# Add any Sphinx extension module names here, as strings. They can be
+# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom
+# ones.
+extensions = ['sphinx.ext.autodoc',
+ 'sphinx.ext.doctest',
+ 'sphinx.ext.todo',
+ 'sphinx.ext.coverage',
+ 'sphinx.ext.ifconfig',
+ 'sphinx.ext.viewcode',
+ 'sphinx.ext.napoleon',
+ 'sphinx.ext.githubpages']
+
+# Add any paths that contain templates here, relative to this directory.
+templates_path = ['_templates']
+
+# The suffix(es) of source filenames.
+# You can specify multiple suffix as a list of string:
+#
+# source_suffix = ['.rst', '.md']
+source_suffix = '.rst'
+
+# The master toctree document.
+master_doc = 'index'
+
+# General information about the project.
+project = u'pylidc'
+copyright = u'2018, Matt Hancock'
+author = u'Matt Hancock'
+
+# The version info for the project you're documenting, acts as replacement for
+# |version| and |release|, also used in various other places throughout the
+# built documents.
+#
+# The short X.Y version.
+version = u''
+# The full version, including alpha/beta/rc tags.
+release = u''
+
+# The language for content autogenerated by Sphinx. Refer to documentation
+# for a list of supported languages.
+#
+# This is also used if you do content translation via gettext catalogs.
+# Usually you set "language" from the command line for these cases.
+language = None
+
+# List of patterns, relative to source directory, that match files and
+# directories to ignore when looking for source files.
+# This patterns also effect to html_static_path and html_extra_path
+exclude_patterns = []
+
+# The name of the Pygments (syntax highlighting) style to use.
+pygments_style = 'sphinx'
+
+# If true, `todo` and `todoList` produce output, else they produce nothing.
+todo_include_todos = True
+
+
+# -- Options for HTML output ----------------------------------------------
+
+# The theme to use for HTML and HTML Help pages. See the documentation for
+# a list of builtin themes.
+#
+#html_theme = 'haiku'
+#html_theme = 'alabaster'
+html_theme = 'bizstyle'
+#html_theme = 'pyramid'
+#html_theme = 'nature'
+#html_theme = 'traditional'
+#html_theme = 'scrolls'
+#html_theme = 'sphinxdoc'
+
+# Theme options are theme-specific and customize the look and feel of a theme
+# further. For a list of options available for each theme, see the
+# documentation.
+#
+html_theme_options = {
+}
+
+# Add any paths that contain custom static files (such as style sheets) here,
+# relative to this directory. They are copied after the builtin static files,
+# so a file named "default.css" will overwrite the builtin "default.css".
+html_static_path = ['_static']
+
+# Custom sidebar templates, must be a dictionary that maps document names
+# to template names.
+#
+# This is required for the alabaster theme
+# refs: http://alabaster.readthedocs.io/en/latest/installation.html#sidebars
+html_sidebars = {
+ '**': [
+ 'relations.html', # needs 'show_related': True theme option to display
+ 'searchbox.html',
+ ]
+}
+
+
+# -- Options for HTMLHelp output ------------------------------------------
+
+# Output file base name for HTML help builder.
+htmlhelp_basename = 'pylidcdoc'
+
+
+# -- Options for LaTeX output ---------------------------------------------
+
+latex_elements = {
+ # The paper size ('letterpaper' or 'a4paper').
+ #
+ # 'papersize': 'letterpaper',
+
+ # The font size ('10pt', '11pt' or '12pt').
+ #
+ # 'pointsize': '10pt',
+
+ # Additional stuff for the LaTeX preamble.
+ #
+ # 'preamble': '',
+
+ # Latex figure (float) alignment
+ #
+ # 'figure_align': 'htbp',
+}
+
+# Grouping the document tree into LaTeX files. List of tuples
+# (source start file, target name, title,
+# author, documentclass [howto, manual, or own class]).
+latex_documents = [
+ (master_doc, 'pylidc.tex', u'pylidc Documentation',
+ u'Matt Hancock', 'manual'),
+]
+
+
+# -- Options for manual page output ---------------------------------------
+
+# One entry per manual page. List of tuples
+# (source start file, name, description, authors, manual section).
+man_pages = [
+ (master_doc, 'pylidc', u'pylidc Documentation',
+ [author], 1)
+]
+
+
+# -- Options for Texinfo output -------------------------------------------
+
+# Grouping the document tree into Texinfo files. List of tuples
+# (source start file, target name, title, author,
+# dir menu entry, description, category)
+texinfo_documents = [
+ (master_doc, 'pylidc', u'pylidc Documentation',
+ author, 'pylidc', 'One line description of project.',
+ 'Miscellaneous'),
+]
diff --git a/docs/source/contour.rst b/docs/source/contour.rst
new file mode 100644
index 0000000..2920703
--- /dev/null
+++ b/docs/source/contour.rst
@@ -0,0 +1,5 @@
+The Contour Model
+====================
+
+.. autoclass:: pylidc.Contour
+ :members:
diff --git a/docs/source/index.rst b/docs/source/index.rst
new file mode 100644
index 0000000..50d8a78
--- /dev/null
+++ b/docs/source/index.rst
@@ -0,0 +1,18 @@
+pylidc
+======
+
+API
+===
+.. toctree::
+ :maxdepth: 7
+
+ scan
+ annotation
+ contour
+
+
+Indices and tables
+==================
+
+* :ref:`genindex`
+* :ref:`search`
diff --git a/docs/source/scan.rst b/docs/source/scan.rst
new file mode 100644
index 0000000..08b7705
--- /dev/null
+++ b/docs/source/scan.rst
@@ -0,0 +1,5 @@
+The Scan Model
+==============
+
+.. autoclass:: pylidc.Scan
+ :members:
diff --git a/pylidc/Annotation.py b/pylidc/Annotation.py
index af83afd..9d5217b 100644
--- a/pylidc/Annotation.py
+++ b/pylidc/Annotation.py
@@ -4,7 +4,6 @@
from ._Base import Base
from .Scan import Scan
-import dicom
import numpy as np
import matplotlib.pyplot as plt
@@ -50,28 +49,59 @@ class Annotation(Base):
has many contours, each of which refers to the contour drawn for
nodule in each scan slice.
- Example:
- >>> import pylidc as pl
- >>> # Get the first annotation with spiculation value greater than 3.
- >>> ann = pl.query(pl.Annotation).filter(pl.Annotation.spiculation > 3).first()
- >>>
- >>> print(ann.spiculation)
- >>> # => 4
- >>>
- >>> # Each nodule feature has a corresponding property
- >>> # to print the semantic value.
- >>> print(ann.Spiculation)
- >>> # => Medium-High Spiculation
- >>>
- >>> q = pl.query(pl.Annotation).join(pl.Scan)
- >>> q = q.filter(pl.Scan.slice_thickness <= 1,
- >>> pl.Annotation.malignancy == 5)
- >>> print(q.count())
- >>> # => 58
- >>>
- >>> ann = q.first()
- >>> print("%.2f, %.2f, %.2f" % (ann.diameter, ann.surface_area, ann.volume))
- >>> # => 17.98, 1221.40, 1033.70
+ Parameters
+ ----------
+ subtlety: int
+ blah
+
+ internalStructure: int
+ blah
+
+ calcification: int
+ blah
+
+ sphericity: int
+ blah
+
+ margin: int
+ blah
+
+ lobulation: int
+ blah
+
+ spiculation: int
+ blah
+
+ texture: int
+ blah
+
+ malignancy: int
+ blah
+
+ Example
+ ---------
+ ::
+ import pylidc as pl
+ # Get the first annotation with spiculation value greater than 3.
+ ann = pl.query(pl.Annotation).filter(pl.Annotation.spiculation > 3).first()
+
+ print(ann.spiculation)
+ # => 4
+
+ # Each nodule feature has a corresponding property
+ # to print the semantic value.
+ print(ann.Spiculation)
+ # => Medium-High Spiculation
+
+ q = pl.query(pl.Annotation).join(pl.Scan)
+ q = q.filter(pl.Scan.slice_thickness <= 1,
+ pl.Annotation.malignancy == 5)
+ print(q.count())
+ # => 58
+
+ ann = q.first()
+ print("%.2f, %.2f, %.2f" % (ann.diameter, ann.surface_area, ann.volume))
+ # => 17.98, 1221.40, 1033.70
"""
__tablename__ = 'annotations'
id = sq.Column('id', sq.Integer, primary_key=True)
@@ -238,45 +268,200 @@ def print_formatted_feature_table(self):
print('%-18s | %-24s | %-2d'%(fnames[i].title(),
fstrings[i], fvals[i]))
-
- def bbox(self, image_coords=True):
+ def bbox(self, pad=None):
"""
- Return a 3 by 2 matrix, corresponding to the bounding box of the
- annotation within the scan. If `scan_slice` is a numpy array
- containing aslice of the scan, each slice of the annotation is
- contained within the box:
-
- bbox[0,0]:bbox[0,1]+1, bbox[1,0]:bbox[1,1]+1
-
- If `image_coords` is `False` then each annotation slice is
- instead contained within:
+ Returns a tuple of Python `slice` objects that can be used to index
+ into the image volume corresponding to the extent of the
+ (padded) bounding box.
+
+ Parameters
+ ----------
+
+ pad: int, list of ints, or float, default=None
+ * If None (default), then no padding is used.
+ * If an integer is provided, then the slices are padded
+ uniformly by this integer amount.
+ * If a list of integers is provided, then it is of the form
+ `[(i1,i2), (j1,j2), (k1,k2)]`
+ and indicates the pad amounts along each coordinate axis.
+ * If a float is provided, then the slices are padded such
+ that the bounding box occupies at least `pad` physical units
+ (using the corresponding scan `pixel_spacing` and `slice_spacing`
+ parameters). This means the returned Slice indices will
+ yield a bounding box that is at least `pad` millimeters along
+ each coordinate axis direction.
+
+ Note
+ ----
+ In the various `pad` cases above, borders are handled so that if a
+ pad beyond the image borders is requested, then it is set
+ to the maximum (or minimum, depending on the direction)
+ possible index.
+
+ Examples
+ --------
+
+ The example below illustrates the various `pad` argument types::
+
+ import pylidc as pl
+
+ ann = pl.query(pl.Annotation).first()
+ vol = ann.scan.to_volume()
+
+ print ann.bbox()
+ # => (slice(151, 185, None), slice(349, 376, None), slice(44, 50, None))
+
+ print(vol[ann.bbox()].shape)
+ # => (34, 27, 6)
- bbox[1,0]:bbox[1,1]+1, bbox[0,0]:bbox[0,1]+1
+ print(vol[ann.bbox(pad=2)].shape)
+ # => (38, 31, 10)
+
+ print(vol[ann.bbox(pad=[(1,2), (3,0), (2,4)])].shape)
+ # => (37, 30, 12)
+
+ print(max(ann.bbox_dims()))
+ # => 21.45
+
+ print(vol[ann.bbox(pad=30.0)].shape)
+ # => (48, 49, 12)
+
+ print(ann.bbox_dims(pad=30.0))
+ # => [30.55, 31.200000000000003, 33.0]
+ """
+ # Error checking ...
+ if pad is not None:
+ if not isinstance(pad, (int, list, float)):
+ raise TypeError("`pad` is incorrect type.")
+ if isinstance(pad, list):
+ if len(pad) != 3:
+ raise ValueError("`pad` list length should be 3.")
+ for p in pad:
+ msg = "`pad` list elements should be (int, int)"
+ if len(p) != 2:
+ raise ValueError(msg)
+ if not isinstance(p[0], int) or not isinstance(p[1], int):
+ raise TypeError(msg)
+
+ # The index limits for the scan.
+ limits = [(0,511), (0,511), (0,self.scan.slice_zvals.shape[0]-1)]
+
+ cmatrix = self.contours_matrix
+ imin,jmin,kmin = cmatrix.min(axis=0)
+ imax,jmax,kmax = cmatrix.max(axis=0)
+
+ # Adding the padding for each respective case, handling the
+ # borders as needed.
+ if isinstance(pad, int):
+ imin = max(imin-pad, limits[0][0])
+ imax = min(imax+pad, limits[0][1])
+ jmin = max(jmin-pad, limits[1][0])
+ jmax = min(jmax+pad, limits[1][1])
+ kmin = max(kmin-pad, limits[2][0])
+ kmax = min(kmax+pad, limits[2][1])
+ elif isinstance(pad, list):
+ imin = max(imin-pad[0][0], limits[0][0])
+ imax = min(imax+pad[0][1], limits[0][1])
+ jmin = max(jmin-pad[1][0], limits[1][0])
+ jmax = min(jmax+pad[1][1], limits[1][1])
+ kmin = max(kmin-pad[2][0], limits[2][0])
+ kmax = min(kmax+pad[2][1], limits[2][1])
+ elif isinstance(pad, float):
+ # In this instance, we compute the extend the limits
+ # until the required physical size is met (or until we can
+ # no long extend the index).
+ rij = self.scan.pixel_spacing
+ rk = self.scan.slice_spacing
+
+ # Check if the desired bbox size is not smaller than is possible.
+ if isinstance(pad, float):
+ minsize = max(self.bbox_dims(pad=None))
+ if pad < minsize:
+ raise ValueError(("Requested `bbox` size (%.4f mm) is "
+ "less than minimal possible size "
+ "(%.4f mm).") % (pad, minsize))
+ while (imax-imin)*rij < pad:
+ imin -= 1 if imin > limits[0][0] else 0
+ imax += 1 if imax < limits[0][1] else 0
+ if imin == limits[0][0] and imax == limits[0][1]:
+ break
+ while (jmax-jmin)*rij < pad:
+ jmin -= 1 if jmin > limits[1][0] else 0
+ jmax += 1 if jmax < limits[1][1] else 0
+ if jmin == limits[1][0] and jmax == limits[1][1]:
+ break
+ while (kmax-kmin)*rk < pad:
+ kmin -= 1 if kmin > limits[2][0] else 0
+ kmax += 1 if kmax < limits[2][1] else 0
+ if kmin == limits[2][0] and kmax == limits[2][1]:
+ break
+
+ return (slice(imin,imax+1),
+ slice(jmin,jmax+1),
+ slice(kmin,kmax+1))
+
+
+ def bbox_dims(self, pad=None):
+ """
+ Return the physical dimensions of the nodule bounding box in mm.
- The last row of `bbox` give the inclusive lower and upper
- bounds of the `image_z_position`.
+ pad: int, list, or float, default=None
+ See :meth:`pylidc.Annotation.bbox` for a
+ description of this argument.
"""
- matrix = self.contours_to_matrix()
- bbox = np.c_[matrix.min(axis=0), matrix.max(axis=0)]
- return bbox if not image_coords else bbox[[1,0,2]]
+ res = [self.scan.pixel_spacing,]*2 + [self.scan.slice_spacing]
+ return [(b.stop-1-b.start)*r for r,b in zip(res, self.bbox(pad=pad))]
+
- def bbox_dimensions(self, image_coords=True):
+ def bbox_matrix(self, pad=None):
"""
- Return the dimensions of the nodule bounding box in mm.
+ The `bbox` function returns a tuple of slices to be used to index
+ into an image volume. On the other hand, `bbox_array` returns
+ a 3x2 matrix where each row is the (start, stop) indices of the
+ i, j, and k axes.
+
+ pad: int, list, or float
+ See `Annotation.bbox` for description of this argument.
+
+ NOTE: The indices return by `bbox_array` are *inclusive*, whereas
+ the indices of the slice objects in the tuple return by `bbox`
+ are offset by +1 in the "stop" index.
+
+ Example:
+ >>> import pylidc as pl
+ >>> ann = pl.query(pl.Annotation).first()
+ >>>
+ >>> bb = ann.bbox()
+ >>> bm = ann.bbox_matrix()
+ >>>
+ >>> print(all([bm[i,0] == bb[i].start for i in range(3)]))
+ >>> # => True
+ >>>
+ >>> print(all([bm[i,1]+1 == bb[i].stop for i in range(3)]))
+ >>> # => True
"""
- bb = self.bbox(image_coords)
- df = np.diff(bb)[:,0]
- df[:2] = df[:2]*self.scan.pixel_spacing
- return df
+ return np.array([[sl.start, sl.stop-1] for sl in self.bbox(pad=pad)])
@property
def centroid(self):
"""
Return the center of mass of the nodule as determined by its
- radiologist-drawn contours. Note that coordinates are image
- domain coordinates (i.e., (i,j,k), not (x,y,z)).
+ radiologist-drawn contours.
+
+ Example:
+ >>> import pylidc as pl
+ >>> import matplotlib.pyplot as plt
+ >>>
+ >>> ann = pl.query(pl.Annotation).first()
+ >>> i,j,k = ann.centroid.round().astype(np.int)
+ >>> vol = ann.scan.to_volume()
+ >>>
+ >>> plt.imshow(vol[:,:,k], cmap=plt.cm.gray)
+ >>> plt.plot(j, i, '.r', label="Nodule centroid")
+ >>> plt.legend()
+ >>> plt.show()
"""
- return self.contours_to_matrix().mean(axis=0)
+ return self.contours_matrix.mean(axis=0)
@property
def diameter(self):
@@ -318,13 +503,13 @@ def surface_area(self):
Estimate the surface area by summing the areas of a trianglation
of the nodules surface in 3d. Returned units are mm^2.
"""
- mask = self.get_boolean_mask()
+ mask = self.boolean_mask()
mask = np.pad(mask, [(1,1), (1,1), (1,1)], 'constant') # Cap the ends.
- dist = dtrans(mask) - dtrans(~mask)
+ mask = mask.astype(np.float)
- rxy = self.scan.pixel_spacing
- rz = self.scan.slice_thickness
- verts, faces, _, _ = marching_cubes(dist, 0, spacing=(rxy, rxy, rz))
+ rij = self.scan.pixel_spacing
+ rk = self.scan.slice_thickness
+ verts, faces, _, _ = marching_cubes(mask, 0.5, spacing=(rij, rij, rk))
return mesh_surface_area(verts, faces)
@property
@@ -414,7 +599,7 @@ def visualize_in_3d(self, edgecolor='0.2', cmap='viridis',
if cmap not in plt.cm.cmap_d.keys():
raise ValueError("Invalid `cmap`. See `plt.cm.cmap_d.keys()`.")
- mask = self.get_boolean_mask()
+ mask = self.boolean_mask()
mask = np.pad(mask, [(1,1), (1,1), (1,1)], 'constant') # Cap the ends.
rxy = self.scan.pixel_spacing
@@ -449,9 +634,9 @@ def visualize_in_3d(self, edgecolor='0.2', cmap='viridis',
elif backend == 'mayavi':
try:
from mayavi import mlab
- sf = mlab.pipeline.scalar_field(dist)
+ sf = mlab.pipeline.scalar_field(mask.astype(np.float))
sf.spacing = [rxy, rxy, rz]
- mlab.pipeline.iso_surface(sf, contours=[0.0])
+ mlab.pipeline.iso_surface(sf, contours=[0.5])
mlab.show()
except ImportError:
print("Mayavi could not be imported. Is it installed?")
@@ -488,7 +673,7 @@ def visualize_in_scan(self, verbose=True):
# plots every time we update the image.
for i,c in enumerate(contours):
arr = c.to_matrix()
- cc, = ax_image.plot(arr[:,0], arr[:,1], '-r')
+ cc, = ax_image.plot(arr[:,1], arr[:,0], '-r')
cc.set_visible(i==0) # Set the first contour visible.
contour_lines.append( cc )
ax_image.set_xlim(-0.5,511.5); ax_image.set_ylim(511.5,-0.5)
@@ -589,67 +774,106 @@ def update_contours(_):
plt.show()
- def contours_to_matrix(self, image_coords=True):
+ @property
+ def contour_slice_zvals(self):
+ """Returns an array of unique z-coordinates for the contours."""
+ return np.unique([c.image_z_position for c in self.contours])
+
+ @property
+ def contour_slice_indices(self):
"""
- Return all the contours in a 3D numpy array.
+ Returns a list of indices into the scan where each contour
+ belongs. An example should clarify:
- image_coords: bool, default True
- If True, the first two coordinates of each point are given
- in image coordinates (i.e., 2d array index). If False, these
- values are scaled by the respective `pixel_spacing` attribute
- of the Annotation's Scan.
+ >>> import pylidc as pl
+ >>>
+ >>> ann = pl.query(pl.Annotation)
+ >>>
+ >>> zvals = ann.contour_slice_zvals
+ >>> kvals = ann.contour_slice_indices
+ >>> scan_zvals = ann.scan.slice_zvals
+ >>>
+ >>> for k,z in zip(kvals, zvals):
+ >>> # the two z values should the same (up to machine precision)
+ >>> print(k, z, scan_zvals[k])
"""
- pts = np.vstack([c.to_matrix() for c in self.contours])
- if not image_coords:
- pts[:,:2] *= self.scan.pixel_spacing
- return pts
-
+ slc = self.scan.slice_zvals
+ return [np.abs(slc-z).argmin() for z in self.contour_slice_zvals]
- def get_boolean_mask(self, return_bbox=False):
+ @property
+ def contours_matrix(self):
"""
- Return a boolean volume which corresponds to the bounding box
- containing the nodule annotation. The slices of the volume are
- ordered by increasing `image_z_position` of the contour
- annotations.
-
- Note that this method doesn't account for a case where the nodule
- contour annotations "skip a slice".
-
- returns: mask, bounding_box
- `mask` is the boolean volume. In the original
- 512 x 512 x num_slices dicom volume, `mask` is a boolean
- mask over the region, `bbox[i,0]:bbox[i,1]+1`, i=0,1,2
+ Return all the contours in a 3D numpy array.
"""
- bbox = self.bbox(image_coords=False)
- zs = np.unique([c.image_z_position for c in self.contours])
- z_to_index = dict(zip(zs,range(len(zs))))
+ return np.vstack([c.to_matrix(include_k=True)
+ for c in sorted(self.contours, key=lambda c: c.image_z_position)])
+
+ def boolean_mask(self, pad=None, bbox=None):
+ """
+ A boolean volume where 1 indicates nodule and 0 indicates
+ non-nodule. The `mask` volume covers the extent of the voxels
+ in the image volume given by `annotation.bbox`, i.e., the `mask`
+ volume would be placed in the full image volume according to
+ the `bbox` attribute.
+
+ pad: int, list, or float, default=None
+ See `Annotation.bbox` for description of this argument.
+
+ bbox: 3x2 NumPy array, default=None
+ If `bbox` is provided, then `pad` is ignored. This argument allows
+ for more fine-tuned control of placement of the mask in a volume,
+ or for pre-computation of bbox when working with multiple
+ Annotation object.
+
+ Example:
+ >>> import pylidc as pl
+ >>> import matplotlib.pyplot as plt
+ >>>
+ >>> ann = pl.query(pl.Annotation).first()
+ >>> vol = ann.scan.to_volume()
+ >>>
+ >>> mask = ann.boolean_mask()
+ >>> bbox = ann.bbox()
+ >>>
+ >>> print("Avg HU inside nodule: %.1f" % vol[bbox][mask].mean())
+ >>> # => Avg HU inside nodule: -280.0
+ >>> print("Avg HU outside nodule: %.1f" % vol[bbox][~mask].mean())
+ >>> # => Avg HU outside nodule: -732.2
+ """
+ bb = self.bbox_matrix(pad=pad) if bbox is None else bbox
+
+ czs = self.contour_slice_zvals
+ cks = self.contour_slice_indices
+
+ zs = self.scan.slice_zvals
+ zs = zs[cks[0]:cks[-1]+1]
+
+ # Lambda to map a z-value to its appropriate index in the volume.
+ z_to_index = lambda z: dict(zip(czs,cks))[z] - bb[2,0]#cks[0]
# Get dimensions, initialize mask.
- nx,ny = np.diff(bbox[:2], axis=1).astype(int) + 1
- nx = int(nx); ny = int(ny)
- nz = int(zs.shape[0])
- mask = np.zeros((nx,ny,nz), dtype=np.bool)
+ ni,nj,nk = np.diff(bb, axis=1).astype(int)[:,0] + 1
+ mask = np.zeros((ni,nj,nk), dtype=np.bool)
# We check if these points are enclosed within each contour
# for a given slice. `test_points` is a list of image coordinate
# points, offset by the bounding box.
- test_points = bbox[:2,0] + np.c_[ np.where(~mask[:,:,0]) ]
+ ii,jj = np.indices(mask.shape[:2])
+ test_points = bb[:2,0] + np.c_[ii.flatten(), jj.flatten()]
# First we "turn on" pixels enclosed by inclusion contours.
for contour in self.contours:
if contour.inclusion:
- zi = z_to_index[contour.image_z_position]
- contour_matrix = contour.to_matrix()[:,:2]
+ zi = z_to_index(contour.image_z_position)
+ C = contour.to_matrix(include_k=False)
# Turn the contour closed if it's not.
- if (contour_matrix[0] != contour_matrix[-1]).any():
- contour_matrix = np.append(contour_matrix,
- contour_matrix[0].reshape(1,2),
- axis=0)
+ if (C[0] != C[-1]).any():
+ C = np.append(C, c[0].reshape(1,2), axis=0)
# Create path object and test all pixels
# within the contour's bounding box.
- path = mplpath.Path(contour_matrix, closed=True)
+ path = mplpath.Path(C, closed=True)
contains_pts = path.contains_points(test_points)
contains_pts = contains_pts.reshape(mask.shape[:2])
# The logical or here prevents the cases where a single
@@ -659,26 +883,19 @@ def get_boolean_mask(self, return_bbox=False):
# Second, we "turn off" pixels enclosed by exclusion contours.
for contour in self.contours:
if not contour.inclusion:
- zi = z_to_index[contour.image_z_position]
- contour_matrix = contour.to_matrix()[:,:2]
+ zi = z_to_index(contour.image_z_position)
+ C = contour.to_matrix(include_k=False)
# Turn the contour closed if it's not.
- if (contour_matrix[0] != contour_matrix[-1]).any():
- contour_matrix = np.append(contour_matrix,
- contour_matrix[0].reshape(1,2),
- axis=0)
+ if (C[0] != C[-1]).any():
+ C = np.append(C, C[0].reshape(1,2), axis=0)
- path = mplpath.Path(contour_matrix, closed=True)
+ path = mplpath.Path(C, closed=True)
not_contains_pts = ~path.contains_points(test_points)
not_contains_pts = not_contains_pts.reshape(mask.shape[:2])
mask[:,:,zi] = np.logical_and(mask[:,:,zi], not_contains_pts)
- # The first and second axes have to
- # be swapped because of the reshape.
- if return_bbox:
- return mask.swapaxes(0,1), bbox[[1,0,2]]
- else:
- return mask.swapaxes(0,1)
+ return mask
def _as_set(self):
"""
@@ -773,8 +990,8 @@ def uniform_cubic_resample(self, side_length=None, resample_vol=True,
returned, depending on which flags are set (see above).
Example:
- >>> self = pl.query(pl.Annotation).first()
- >>> ct_volume, mask = self.uniform_cubic_resample(side_length=70)
+ >>> ann = pl.query(pl.Annotation).first()
+ >>> ct_volume, mask = ann.uniform_cubic_resample(side_length=70)
>>> print(ct_volume.shape, mask.shape)
>>> # => (71, 71, 71), (71, 71, 71)
>>> # (Nodule is centered at (35,35,35).)
@@ -783,8 +1000,8 @@ def uniform_cubic_resample(self, side_length=None, resample_vol=True,
>>> plt.imshow( ct_volume[:,:,35] * (0.2 + 0.8*mask[:,:,35]) )
>>> plt.show()
"""
- bbox = self.bbox(image_coords=True)
- bboxd = self.bbox_dimensions(image_coords=True)
+ bbox = self.bbox
+ bboxd = self.bbox_dimensions
rxy = self.scan.pixel_spacing
imin,imax = bbox[0].astype(int)
@@ -799,7 +1016,7 @@ def uniform_cubic_resample(self, side_length=None, resample_vol=True,
side_length = np.ceil(bboxd.max())
else:
if not isinstance(side_length, int):
- raise ValueError('`side_length` must be an integer.')
+ raise TypeError('`side_length` must be an integer.')
if side_length < bboxd.max():
raise ValueError('`side_length` must be greater\
than any bounding box dimension.')
@@ -820,7 +1037,7 @@ def uniform_cubic_resample(self, side_length=None, resample_vol=True,
kmax = np.where(zmax == img_zs)[0][0]
# Initialize the boolean mask.
- mask = self.get_boolean_mask()
+ mask = self.boolean_mask()
########################################################
# { Begin mask corrections.
diff --git a/pylidc/Contour.py b/pylidc/Contour.py
index 2633260..4711b57 100644
--- a/pylidc/Contour.py
+++ b/pylidc/Contour.py
@@ -10,8 +10,8 @@
class Contour(Base):
"""
- The Contour model class holds a contour for a single scan slice
- of a Nodule object, as drawn by the annotating radiologist.
+ The Contour class holds the coordinate data for the vertices of
+ a :class:`pylidc.Annotation` object for a single slice in the CT volume.
inclusion: bool
If True, the area inside the contour is included as part of
@@ -35,12 +35,16 @@ class Contour(Base):
Example:
>>> import pylidc as pl
- >>> scan = pl.query(pl.Scan).first()
- >>> ann = scan.annotations[0]
-
- >>> for c in ann.contours:
- >>> print(c.image_z_position, c.to_matrix().mean(axis=0))
-
+ >>> import matplotlib.pyplot as plt
+ >>>
+ >>> ann = pl.query(pl.Annotation).first()
+ >>> vol = ann.scan.to_volume()
+ >>> ii,jj = ann.contours[0].to_matrix(include_k=False).T
+ >>>
+ >>> plt.imshow(vol[:,:,46], cmap=plt.cm.gray)
+ >>> plt.plot(jj, ii, '-r', lw=1, label="Nodule Boundary")
+ >>> plt.legend()
+ >>> plt.show()
"""
__tablename__ = 'contours'
id = sq.Column('id', sq.Integer, primary_key=True)
@@ -62,14 +66,31 @@ def __setattr__(self, name, value):
else:
super(Contour,self).__setattr__(name,value)
- def to_matrix(self):
+ def to_matrix(self, include_k=True):
"""
Return the contour-annotation coordinates as a matrix where
- each row contains an (x,y,z) coordinate.
+ each row contains an (i,j,k) *index* coordinate into the image volume.
+
+ include_k: bool, default=True
+ Set `include_k=False` to omit the `k` axis coordinate. In this
+ case, each row of the returned matrix is an (i,j) *index*
+ coordinate into the `k`th image. The index `k` can be
+ recovered via:
+ z = contour.image_z_position
+ k = contour.annotation.contour_slice_indexes[z]
"""
- xy = np.array([list(map(int,c.split(','))) for c in self.coords.split('\n')])
- z = np.ones(xy.shape[0])*self.image_z_position
- return np.c_[xy,z]
+ # The reversal [::-1] is because the coordinates from the LIDC XML
+ # are stored as (x,y), not (i,j).
+ ij = np.array([[int(cc) for cc in c.split(',')][::-1]
+ for c in self.coords.split('\n')])
+ if not include_k:
+ return ij
+ else:
+ k = np.ones(ij.shape[0])
+ zs = self.annotation.contour_slice_zvals
+ kk = np.abs(zs-self.image_z_position).argmin()
+ k *= self.annotation.contour_slice_indices[kk]
+ return np.c_[ij, k].astype(np.int)
Annotation.contours = relationship('Contour',
order_by=Contour.id,
diff --git a/pylidc/Scan.py b/pylidc/Scan.py
index 0b779df..e5dd57d 100644
--- a/pylidc/Scan.py
+++ b/pylidc/Scan.py
@@ -5,7 +5,7 @@
from ._Base import Base
import os, warnings
-import dicom
+import pydicom as dicom
import numpy as np
import matplotlib.pyplot as plt
@@ -15,9 +15,10 @@
from scipy.sparse.csgraph import connected_components
from .annotation_distance_metrics import metrics
+from scipy.stats import mode
-# Load the configuration file and get the dicom file path.
+# Load the configuration file and get the dicom file path.
try:
import configparser
except ImportError:
@@ -60,40 +61,56 @@
class Scan(Base):
"""
The Scan model class refers to the top-level XML file from the LIDC.
- A scan has many annotations, each of which is the `unblindedReadNodule`
- for the scan.
-
- Attributes:
- study_instance_uid: string
- series_instance_uid: string
+ A scan has many :class:`pylidc.Annotation` objects, which correspond
+ to the `unblindedReadNodule` XML attributes for the scan.
- patient_id: string
- Identifier of the from `LIDC-IDRI-#`
+ Attributes
+ ==========
- slice_thickness: float
- Dicom attribute, `(0018,0050)`. Note that this may not be
- equal to `spacing_betwee_slices`.
+ study_instance_uid: string
+ series_instance_uid: string
- pixel_spacing: float
- Dicom attribute, `(0028,0030)`. This is normally two
- values. All scans in the LIDC have equal resolutions
- in the transverse plane, so only one value is taken here.
+ patient_id: string
+ Identifier of the from `LIDC-IDRI-####`
- contrast_used: bool
- If the dicom file for the scan had any Contrast tag,
- this is marked as `True`.
+ slice_thickness: float
+ DICOM attribute, `(0018,0050)`. Note that this may not be
+ equal to the `slice_spacing` attribute (see below).
- is_from_initial: bool
- Indicates whether or not this PatientID was tagged as
- part of the initial 399 release.
+ slice_zvals: ndarray
+ The "z-values" for the slices of the scan (i.e.,
+ the last coordinate of the ImagePositionPatient DICOM attribute)
+ as a NumPy array sorted in increasing order.
- sorted_dicom_file_names: string
- This is a string containing a comma-separated list
- like `[number].dcm`. It is a list of the dicom file
- names for this scan in order of increasing z-coordinate
- of dicom attribute, `(0020,0032)`. In rare cases where
- a scan includes multiple files with the same z-coordinate,
- the one with the lesser `InstanceNumber` is used.
+ slice_spacing: float
+ This computes the median of the difference
+ between the slice coordinates, i.e., `scan.slice_zvals`.
+ NOTE: that this attribute is typically (but not always!) the
+ same as the `slice_thickness` attribute. Furthermore,
+ the `slice_spacing` does NOT necessarily imply that all the
+ slices are spaced with spacing (although they often are).
+
+ pixel_spacing: float
+ Dicom attribute, `(0028,0030)`. This is normally two
+ values. All scans in the LIDC have equal resolutions
+ in the transverse plane, so only one value is taken here.
+
+ contrast_used: bool
+ If the DICOM file for the scan had any Contrast tag,
+ this is marked as `True`.
+
+ is_from_initial: bool
+ Indicates whether or not this PatientID was tagged as
+ part of the initial 399 release.
+
+ sorted_dicom_file_names: string
+ This attribute is no longer used, but left for historical
+ reasons. This is a string containing a comma-separated list
+ like `[number].dcm`. It is a list of the dicom file
+ names for this scan in order of increasing z-coordinate
+ of dicom attribute, `(0020,0032)`. In rare cases where
+ a scan includes multiple files with the same z-coordinate,
+ the one with the lesser `InstanceNumber` is used.
Example:
>>> import pylidc as pl
@@ -130,7 +147,7 @@ def __setattr__(self, name, value):
def get_path_to_dicom_files(self):
"""
- Get the path to where the dicom files are stored for this scan,
+ Get the path to where the DICOM files are stored for this scan,
relative to the root path set in the pylidc configuration file (e.g.,
`~/.pylidc` in MAC and Linux).
@@ -192,8 +209,7 @@ def get_path_to_dicom_files(self):
# all have the same series/study ids.
dicom_file = dicom_file[0]
- with open(os.path.join(dpath, dicom_file)) as f:
- dimage = dicom.read_file(f)
+ dimage = dicom.dcmread(os.path.join(dpath, dicom_file))
seid = str(dimage.SeriesInstanceUID).strip()
stid = str(dimage.StudyInstanceUID).strip()
@@ -253,8 +269,6 @@ def cluster_annotations(self, metric='min', tol=None, factor=0.9,
the number of unique physical nodules present in this Scan as
determined by this overlap method and the tolerance used.
- ---
-
More on the `tol` parameter and distance measures:
The "distance" matrix, `d[i,j]`, between all Annotations for
@@ -266,6 +280,22 @@ def cluster_annotations(self, metric='min', tol=None, factor=0.9,
adjacent when `d[i,j] <= tol`. Groups are determined by finding
the connected components of the graph associated with
this adjacency matrix.
+
+ Example::
+
+ import pylidc as pl
+
+ scan = pl.query(pl.Scan).first()
+ nodules = scan.cluster_annotations()
+ print("This can has %d nodules." % len(nodules))
+ # => This can has 4 nodules.
+
+ for i,nod in enumerate(nodules):
+ print("Nodule %d has %d annotations." % (i+1,len(nod)))
+ # => Nodule 1 has 4 annotations.
+ # => Nodule 2 has 4 annotations.
+ # => Nodule 3 has 1 annotations.
+ # => Nodule 4 has 4 annotations.
"""
assert 0 < factor < 1, "`factor` must be in the interval (0,1)."
@@ -334,15 +364,19 @@ def load_all_dicom_images(self, verbose=True):
"""
Load all the DICOM images assocated with this scan and return as list.
- Example:
- >>> scan = pl.query(pl.Scan).first()
- >>> images = scan.load_all_dicom_images()
- >>> zs = [float(img.ImagePositionPatient[2]) for img in images]
- >>> print(zs[1] - zs[0], img.SliceThickness, scan.slice_thickness)
- >>>
- >>> import matplotlib.pyplot as plt
- >>> plt.imshow( images[0].pixel_array, cmap=plt.cm.gray )
- >>> plt.show()
+ Example::
+
+ import pylidc as pl
+ import matplotlib.pyplot as plt
+
+ scan = pl.query(pl.Scan).first()
+
+ images = scan.load_all_dicom_images()
+ zs = [float(img.ImagePositionPatient[2]) for img in images]
+ print(zs[1] - zs[0], img.SliceThickness, scan.slice_thickness)
+
+ plt.imshow( images[0].pixel_array, cmap=plt.cm.gray )
+ plt.show()
"""
if verbose: print("Loading dicom files ... This may take a moment.")
@@ -351,9 +385,8 @@ def load_all_dicom_images(self, verbose=True):
if fname.endswith('.dcm')]
images = []
for fname in fnames:
- with open(os.path.join(path, fname), 'rb') as f:
- image = dicom.read_file(f)
- images.append(image)
+ image = dicom.dcmread(os.path.join(path,fname))
+ images.append(image)
# ##############################################
# Clean multiple z scans.
@@ -392,6 +425,15 @@ def visualize(self, annotation_groups=None):
annotation_groups: list of lists of Annotation objects
This argument should be supplied by the returned object from
the `cluster_annotations` method.
+
+ Example::
+
+ import pylidc as pl
+
+ scan = pl.query(pl.Scan).first()
+ nodules = scan.cluster_annotations()
+
+ scan.visualize(annotation_groups=nodules)
"""
images = self.load_all_dicom_images()
@@ -419,8 +461,8 @@ def visualize(self, annotation_groups=None):
c = centroids[i]
s = '%d Annotations'%len(annotation_groups[i])
a = ax_image.annotate(s,
- xy=(c[0]-r, c[1]-r),
- xytext=(c[0]-50, c[1]-50),
+ xy=(c[1]-r, c[0]-r),
+ xytext=(c[1]-50, c[0]-50),
bbox=dict(fc='w', ec='r'),
arrowprops=dict(arrowstyle='->',
edgecolor='r'))
@@ -454,7 +496,8 @@ def visualize(self, annotation_groups=None):
c = centroids[i]
g = annotation_groups[i]
txt.append(['Nodule %d:'%(i+1),
- '%d annotations, near z=%.2f'%(len(g),c[2])])
+ '%d annotations, near slice %d' \
+ % (len(g), int(c[2].round()))])
ann_grps_table = ax_ann_grps.table(cellText=txt, loc='center',
cellLoc='left')
# Remove cell borders.
@@ -485,9 +528,9 @@ def update(_):
# Show annotation labels if possible.
if annotation_groups is not None:
for i in range(len(annotation_groups)):
- dist = abs(z-centroids[i][2])
- arrows[i].set_visible(dist <= 3*self.slice_thickness)
-
+ centroid_z = self.slice_zvals[int(centroids[i][2].round())]
+ dist = abs(z - centroid_z)
+ arrows[i].set_visible(dist <= 3*self.slice_spacing)
fig.canvas.draw_idle()
sslice.on_changed(update)
@@ -495,7 +538,7 @@ def update(_):
plt.show()
@property
- def slice_coords(self):
+ def slice_zvals(self):
"""
The "z-values" for the slices of the scan (i.e.,
the last coordinate of the ImagePositionPatient DICOM attribute)
@@ -503,6 +546,19 @@ def slice_coords(self):
"""
return np.sort([z.val for z in self.zvals])
+ @property
+ def slice_spacing(self):
+ """
+ This computes the median of the difference
+ between the slice coordinates, i.e., `scan.slice_zvals`.
+
+ NOTE: that this attribute is typically (but not always!) the
+ same as the `slice_thickness` attribute. Furthermore,
+ the `slice_spacing` does NOT necessarily imply that all the
+ slices are spaced with spacing (although they often are).
+ """
+ return np.median(np.diff(self.slice_zvals))
+
def to_volume(self, verbose=True):
"""
Return the scan as a 3D numpy array volume.
diff --git a/pylidc/__init__.py b/pylidc/__init__.py
index 0a4a760..064f3f2 100644
--- a/pylidc/__init__.py
+++ b/pylidc/__init__.py
@@ -1,8 +1,8 @@
"""
---------------------------------------------------------
-Author: Matt Hancock, not.matt.hancock@gmail.com
+pylidc
--------------------------------------------------------
+Author: Matt Hancock, not.matt.hancock@gmail.com
This python module implements an (ORM) object relational mapping
to an sqlite database containing the annotation information from
the XML files provided by the LIDC dataset. The purpose of this
@@ -13,7 +13,7 @@
The ORM is implemented using sqlalchemy. There are three data models:
-Scan, Annotation, and Contour
+:class:`Scan`, :class:`Annotation`, and :class:`Contour`
The relationships are "one to many" for each model going left
to right, i.e., scans have many annotations and annotations
@@ -43,24 +43,29 @@
from .Zval import Zval
from .Annotation import feature_names as annotation_feature_names
-from utils import consensus
+from .utils import consensus
def query(*args):
"""
- Wraps the sqlalchemy session object. Some example usage:
+ Wraps the sqlalchemy session object. Some example usage::
- >>> import pylidc as pl
- >>> qu = pl.query(pl.Scan).filter(pl.Scan.slice_thickness <= 1.)
- >>> print qu.count()
- >>> # => 97
- >>> scan = qu.first()
- >>> print len(scan.annotations)
- >>> # => 11
- >>> qu = pl.query(pl.Annotation).filter((pl.Annotation.malignancy > 3), (pl.Annotation.spiculation < 3))
- >>> print qu.count()
- >>> # => 1083
- >>> annotation = qu.first()
- >>> print annotation.volume
- >>> # => 5230.33874999
+ import pylidc as pl
+
+ qu = pl.query(pl.Scan).filter(pl.Scan.slice_thickness <= 1.)
+ print qu.count()
+ # => 97
+
+ scan = qu.first()
+ print(len(scan.annotations))
+ # => 11
+
+ qu = pl.query(pl.Annotation).filter(pl.Annotation.malignancy > 3,
+ pl.Annotation.spiculation < 3)
+ print(qu.count())
+ # => 1083
+
+ annotation = qu.first()
+ print(annotation.volume)
+ # => 5230.33874999
"""
return _session.query(*args)
diff --git a/pylidc/annotation_distance_metrics.py b/pylidc/annotation_distance_metrics.py
index 0b456c6..fbcabe2 100644
--- a/pylidc/annotation_distance_metrics.py
+++ b/pylidc/annotation_distance_metrics.py
@@ -12,8 +12,8 @@ def pairdist(ann1, ann2, which):
which: str
One of 'min', 'max', or 'avg'.
"""
- dists = cdist(ann1.contours_to_matrix(0),
- ann2.contours_to_matrix(0))
+ dists = cdist(ann1.contours_matrix,
+ ann2.contours_matrix)
if which == 'min':
return dists.min()
@@ -46,8 +46,8 @@ def centroid_xy(ann1, ann2, which):
which: str
One of 'min', 'max', or 'avg'.
"""
- P1 = ann1.contours_to_matrix(0)
- P2 = ann2.contours_to_matrix(0)
+ P1 = ann1.contours_matrix
+ P2 = ann2.contours_matrix
zvals1 = set(np.unique(P1[:,2]).tolist())
zvals2 = set(np.unique(P2[:,2]).tolist())
@@ -82,8 +82,8 @@ def hausdorff(ann1, ann2):
[1]: https://en.wikipedia.org/wiki/Hausdorff_distance
"""
- C = cdist(ann1.contours_to_matrix(0),
- ann2.contours_to_matrix(0))
+ C = cdist(ann1.contours_matrix,
+ ann2.contours_matrix)
return max(C.min(0).max(), C.min(1).max())
metrics['hausdorff'] = hausdorff
diff --git a/pylidc/utils.py b/pylidc/utils.py
index 937712e..b6a6aa1 100644
--- a/pylidc/utils.py
+++ b/pylidc/utils.py
@@ -1,7 +1,7 @@
import numpy as np
from .Annotation import Annotation
-def consensus(anns, clevel=0.5, extent=None, verbose=True):
+def consensus(anns, clevel=0.5, pad=None, ret_masks=True, verbose=True):
"""
Return the boolean-valued consensus volume amongst the
provided annotations (`anns`) at a particular consensus level
@@ -18,59 +18,41 @@ def consensus(anns, clevel=0.5, extent=None, verbose=True):
when >= 50% of the segmentations include that voxel, and 0
otherwise.
- extent: float, default=None
- If a float, the volumes will be padded so that they have
- physical dimensions >= `extent` units in each coordinate
- axis direction. The default (None) doesn't perform this
- operation.
+ pad: int, list, or float, default=None
+ See `Annotation.bbox` for description for this argument.
+
+ ret_masks: bool, default=True
+ If True, a list of masks is also returned corresponding to
+ all the annotations. Note that this slightly different than calling
+ `boolean_mask` on each respective Annotation object because these
+ volumes will be the same shape and in a common reference frame.
verbose: bool, default=True
Turns the DICOM image loading message on/off.
- """
- if not all([isinstance(a, Annotation) for a in anns]):
- raise TypeError("`anns` should be list of `pylidc.Annotation`s.")
-
- clevel = float(clevel)
-
- if not (0.0 <= clevel <= 1.0):
- raise ValueError("`clevel` should be between 0 and 1.")
-
- scan = anns[0].scan
- rij = scan.pixel_spacing
- rk = scan.slice_thickness
-
- # Load the images. Get the z positions.
- vol = scan.to_volume()
- img_zs = scan.slice_coords
-
- bboxs = np.array([a.bbox(image_coords=1) for a in anns])
- # Last dimension is z dist, not index.
- for i in range(bboxs.shape[0]):
- # Float comparison ... probably not the best way.
- bboxs[i,2,0] = img_zs.searchsorted(bboxs[i,2,0])
- bboxs[i,2,1] = img_zs.searchsorted(bboxs[i,2,1])
-
- # Cast to int now that everyone is an index.
- bboxs = bboxs.astype(np.int)
-
- imin,jmin,kmin = np.array([b[:,0] for b in bboxs]).min(0)
- imax,jmax,kmax = np.array([b[:,1] for b in bboxs]).max(0)
-
- while (imax-imin)*rij < extent:
- imin -= 1 if imin > 0 else 0
- imax += 1 if imax < 511 else 0
-
- while (jmax-jmin)*rij < extent:
- jmin -= 1 if jmin > 0 else 0
- jmax += 1 if jmax < 511 else 0
-
- while (kmax-kmin)*rk < extent:
- kmin -= 1 if kmin > 0 else 0
- kmax += 1 if kmax < (img_zs.shape[0]-1) else 0
-
- return_bbox = np.array([[imin, imax],
- [jmin, jmax],
- [kmin, kmax]])
-
- raise NotImplementedError()
+ returns: consensus_mask, consensus_bbox[, masks]
+ `consensus_mask` is the boolean-valued volume of the annotation
+ masks at `clevel` consensus. `consensus_bbox` is a 3-tuple of
+ slices that can be used to index into the image volume at the
+ corresponding location of `consensus_mask`. `masks` is a list of
+ boolean-valued mask volumes corresponding to each Annotation object.
+ Each mask in the `masks` list has the same shape and is sampled in
+ the common reference frame provided by `consensus_bbox`.
+ """
+ bmats = np.array([a.bbox_matrix(pad=pad) for a in anns])
+ imin,jmin,kmin = bmats[:,:,0].min(axis=0)
+ imax,jmax,kmax = bmats[:,:,1].max(axis=0)
+
+ # consensus_bbox
+ cbbox = np.array([[imin,imax],
+ [jmin,jmax],
+ [kmin,kmax]])
+
+ masks = [a.boolean_mask(bbox=cbbox) for a in anns]
+ cmask = np.mean(masks, axis=0) >= clevel
+ cbbox = tuple(slice(cb[0], cb[1]+1, None) for cb in cbbox)
+
+ if ret_masks:
+ return cmask, cbbox, masks
+ else:
+ return cmask, cbbox
diff --git a/setup.py b/setup.py
index df8c62c..b954584 100644
--- a/setup.py
+++ b/setup.py
@@ -32,7 +32,7 @@
packages=find_packages(exclude=['contrib', 'doc', 'tests*']),
install_requires=[
'sqlalchemy>=1.1.5', 'numpy>=1.12.0', 'scipy>=0.18.1',
- 'matplotlib>=2.0.0', 'pydicom>=0.9.9', 'scikit-image>=0.13'
+ 'matplotlib>=2.0.0', 'pydicom>=1.0.0', 'scikit-image>=0.13'
],
package_data={
'pylidc': ['pylidc.sqlite'],