Module Classes

The autodocumentation for the pyCellAnalyst module is provided below.

class pyCellAnalyst.Volume(vol_dir, output_dir=None, regions=None, pixel_dim=[0.411, 0.411, 0.6835], stain='Foreground', segmentation='Geodesic', smoothing_method='Curvature Diffusion', smoothing_parameters={}, two_dim=False, bright=False, enhance_edge=False, depth_adjust=False, display=True, handle_overlap=True, debug=False, opening=True, fillholes=True)

This class will segment objects from 3-D images using user-specified routines. The intended purpose is for laser scanning fluorescence microscopy of chondrocytes and/or their surrounding matrices. Nevertheless, this can be generalized to any 3-D object using any imaging modality; however, it is likely the segmentation parameters will need to be adjusted. Therefore, in this case, the user should set segmentation=’User’ during Class instantiation, and call the segmentaion method with appropriate parameters.

Parameters:
  • vol_dir (str) – This is required. Currently it is the path to a directory containing a stack of TIFF images or a single NifTi (.nii) file. Other formats may be supported in the future.
  • output_dir (str, optional) – The directory in which to save results. If not specified, a directory vol_dir + ‘_results’ will be created and used.
  • regions ([,[int, int, int, int, int, int], ..], optional) – Cropped regions bounding a single object to segment. In terms of voxel indices the order for each region is: [top left corner x, y, z, box edge length Lx, Ly, Lz]. If not specified, the entire image is considered as the region.
  • pixel_dim ([float=0.411, float=0.411, float=0.6835], optional) – The physical dimensions of the voxel ordered x, y, and z. If there is a need to correct a dimesion such as for the depth distortion in laser scanning microscopy, it should be incorporated here. Defaults to [0.411, 0.411, 0.6835]
  • stain (str='Foreground', optional) –
    • ‘Foreground’ indicates the objects of interest appear bright in the image.
    • ’Background’ indicates the objects of interest appear dark.
  • segmentation (str='Geodesic', optional) –
    • ‘Threshold’ – indicates to threshold the image at 0.4\times intensity_{max}.
    • ’Geodesic’ – (default) peform a geodesic active contour segmentation with default settings.
    • ’Edge-Free’ – perform an edge-free active contour segmentation with default dettins.
    • ’User’ – The user will invoke calls to segmentation function with custom settings.
  • smoothing_method (str='Curvature Diffusion', optional) –

    Smoothing method to use on regions of interest.

    • ’None’ – No smoothing will be performed.
    • ’Gaussian’ – Perform Gaussian smoothing.
    • ’Median’ – Apply a median filter.
    • ’Curvature Diffusion’ – Perform curvature-based anisotropic diffusion smoothing.
    • ’Gradient Diffusion’ – Peform classical anisotropic diffusion smoothing.
    • ’Bilateral’ – Apply a bilateral filter.
    • ’Patch-based’ – Perform patch-based denoising.
  • smoothing_parameters (dict, optional) –

    Depends on smoothing_method. Field keys are documented in methods smoothRegion().

    • ’Gaussian’ – fields are:
      • ’sigma’: float=0.5
    • ’Median’ – fields are:
      • ’radius’: (int=1, int=1, int=1)
    • ’Curvature Diffusion’ – fields are:
      • ’iterations’: int=10
      • ’conductance’: float=9.0
    • ’Gradient Diffusion’ – fields are:
      • ’iterations’: int=10
      • ’conductance’: float=9.0
      • ’time step’: float=0.01
    • ’Bilateral’ – fields are:
      • ’domainSigma’: float=1.5
      • ’rangeSigma’: float=40.0
      • ’samples’: int=100
    • ’Patch-based’ – fields are:
      • ’radius’: int=3
      • ’iterations’: int=10
      • ’patches’: int=20
      • ’noise model: str=’poisson’
        • options: (‘none’, ‘gaussian’, ‘poisson’, ‘rician’)
  • two_dim (bool=False, optional) – If True, will consider each 2-D slice in stack independently in smoothing and segmentation. This is not recommended except in special cases.
  • bright (bool=False, optional) – If True, will perform bright spot removal replacing voxels with intensities \ge 98^{th} percentile with median filtered (radius=6) value.
  • enhance_edge (bool=False, optional) – If True, will enhance edges after smoothing using Laplacian sharpening.
  • depth_adjust (bool=False, optional) – If True, will perform a linear correction for intensity degradation with depth.
  • opening (bool=True, optional) – If True, will perform a morphological binary opening following thresholding to remove spurious connections and islands. If object of interest is thin, this may cause problems in which case, setting this to False may help.
  • fillholes (bool=False, optional) – If True, all holes completely internal to the segmented object will be considered as part of the object.
  • display (bool=True, optional) – If True, will spawn a 3-D interactive window rendering of segmented object surfaces.
  • handle_overlap (bool=True, optional) – If True, overlapping segmented objects will be reclassified using a Support Vector Machine.
  • debug (bool=False, optional) – If True, will write additional images to disk in NifTi format for debugging purposes.
cells

SimpleITK Image – An image containing the segmented objects as integer labels. Has the same properties as the input image stack.

thresholds

[, int, …] – The threshold level for each cell

volumes

[, float, …] – List of the physical volumes of the segmented objects.

centroids

[,[float, float, float], …]1 – List of centroids of segmented objects in physical space.

surfaces

[,vtkPolyData, …] – List containing VTK STL objects.

dimensions

[,[float, float, float],…] – List containing the ellipsoid axis lengths of segmented objects. These are determined from the segmented binary images. It is recommended to use the values calculated from a 3-D mesh in the CellMech class.

_classifyShared(i, cells, previous)

If segmented objects overlap and handle_overlap is True, this will attempt to reclassify the shared voxels using the thresholded seed to train a support vector machine. Of course, this relies on the seed to not overlap. The user strategy to get good results from this would be to use an active contour method, with an aggressive thresholding method to produce the seed.

Returns:
  • A modified version of cells attribute with the overlapping objects
  • reclassified.
_flattenBorder(img)

To help ensure the segmentation does not touch the cropped region border, the voxel intensities of the 6 border slices are replaced by the intensity of their 1st percentile.

_getLabelShape(img)
_getMinMax(img)
_parseStack()
_replaceSeed(seed)
adjustForDepth()

Iterates over slices in 3-D image stack and appends each slice’s maximum pixel intensity to a list. Then performs a weighted linear curve fit with slices below the 2^{nd} percentile and above the 98^{th} percentile assigned zero weights with all others equally weighted. Ratios are then calculated from this fit for all z-depths as:

\frac{a_0}{a_1 z + a_0}

and each slice of the image is multiplied by its corresponding weight and reassembled into a 3-D image.

Returns:
  • Replaces image read from disk with and image that is corrected for
  • intensity change with depth.
edgeFreeSegmentation(upsampling=2, seed_method='Percentage', adaptive=True, ratio=0.4, lambda1=1.0, lambda2=1.1, curvature=0.0, iterations=20)

Performs a segmentation using the SimpleITK implementation of the Active Contours Without Edges method described in (Chan and Vese. 2001.) Please also consult SimpleITK’s documentation of ScalarChanAndVeseDenseLevelSetImageFilter.

Parameters:
  • upsampling (int=2, optional) – Resample image splitting original voxels this many times. Resampling will always be performed to make voxels isotropic, because anisotropic voxels can degrade the performance of this algorithm.
  • seed_method (str='Percentage') – Thresholding method used to determine seed image; same as thresholdSegmentation() method parameter. Please consult its documentation.
  • adaptive (bool=True) – If true will adaptively adjust threshold the threshold value until resulting segmentation no longer touches the region of interest bounds.
  • ratio (float=0.7) – The ratio to use with ‘Percentage’ threshold method. This plays no role with other seed methods.
  • lambda1 (float=1.0) – Weight for internal levelset term contribution to the total energy.
  • lambda2 (float=1.1) – Weight for external levelset term contribution to the total energy.
  • curvature (float=0.0) – Weight for curvature. Higher results in smoother levelsets, but less ability to capture fine features.
  • iterations (int=20) – The number of iterations the active contour method will conduct.
geodesic2D(seed, simg, cannyLower, cannyUpper, canny_variance, upsampling, active_iterations, rms, propagation, curvature, advection)

A 2-D implementation of geodesicSegmentation() that operates on each slice in the 3-D stack independently.

geodesicSegmentation(upsampling=2, seed_method='Percentage', adaptive=True, ratio=0.7, canny_variance=(0.05, 0.05, 0.05), cannyUpper=0.0, cannyLower=0.0, propagation=0.15, curvature=0.2, advection=1.0, rms=0.01, active_iterations=200)

Performs a segmentation using the SimpleITK implementation of the Geodesic Active Contour Levelset Segmentation method described in (Caselles et al. 1997.) Please also consult SimpleITK’s documentation of GeodesicActiveContourLevelSetImageFilter. This method will establish the initial levelset function by calling the thresholdSegmentation() method, and calculating a distance map from the resulting binary image.

Parameters:
  • propagation (float=0.15) – Weight for propagation term in active contour functional. Higher values result in faster expansion.
  • curvature (float=0.2) – Weight for curvature term in active contour functional. Higher values result in smoother segmentation.
  • advection (float=1.0) – Weight for advective term in active contour functional. Higher values causes the levelset evolution to be drawn and stick to edges.
  • rms (float=0.01) – The change in root-mean-square difference at which iterations will terminate. This value is divided by the upsampling value to account the effect of voxel size.
  • active_iterations (int=200) – The maximum number of iterations the active contour will conduct.
  • upsampling (int=2, optional) – Resample image splitting original voxels this many times. Resampling will always be performed to make voxels isotropic, because anisotropic voxels can degrade the performance of this algorithm.
  • seed_method (str='Percentage') – Thresholding method used to determine seed image; same as thresholdSegmentation() method parameter. Please consult its documentation.
  • adaptive (bool=True) – If true will adaptively adjust threshold the threshold value until resulting segmentation no longer touches the region of interest bounds.
  • ratio (float=0.7) – The ratio to use with ‘Percentage’ threshold method. This plays no role with other seed methods.
  • canny_variance ([float=0.05, float=0.05, float=0.05]) – The Gaussian variance for canny edge detection used to generate the edge map for this method. Gaussian smoothing is performed during edge detection, but if another smoothing method was already performed this can be set low. High values results in smoother edges, but risk losing edges when other objects are close.
  • cannyUpper (float=0.0) – Ensures voxels in the image gradient with a value higher than this will always be considered edges, and never discarded.
  • cannyLower (float=0.0) – Ensures voxels in the image gradient with a value lower than this will be discarded.
getDimensions()
scale2D(img, thresh)
smooth2D(img)
smoothRegion(img)
threshold2D(img, thres, ratio)
thresholdSegmentation(method='Percentage', adaptive=True, ratio=0.4)

Segments object of interest from image using user-specified method.

Parameters:
  • method (str='Percentage') –

    The thresholding method to use. Options are:

    • ’Percentage’ – Threshold at percentage of the maximum voxel intensity.
    • ’Otsu’ – Threshold using Otsu’s method
    • ’Huang’ –
    • ‘IsoData’
    • ’Li’
    • ’MaxEntropy’ – Sets the threshold value such that the sum of information entropy (Shannon) in the foreground and background is maximized.
    • ’KittlerIllingworth’
    • ’Moments’
    • ’Yen’
    • ’RenyiEntropy’ – The same as ‘MaxEntropy’, but uses the Renyi entropy function.
    • ’Shanbhag’ – Extends upon the entropy methods with fuzzy set theory.
  • ratio (float=0.4) – Ratio of maximum voxel intensity in the region of interest to threshold at. Only used if ‘Percentage’ method is given.
  • adaptive (bool=True) – If True will adaptively adjust the determined threshold value until the segmented object does not the region boundaries.
writeLabels()
writeSurfaces()
class pyCellAnalyst.CellMech(ref_dir=None, def_dir=None, rigidInitial=True, deformable=False, saveFEA=False, deformableSettings={'Maximum RMS': 0.01, 'Iterations': 200, 'Displacement Smoothing': 3.0, 'Precision': 0.01}, display=False)

Quantifies deformation between objects in reference and deformed states.

Object geometries are obtained by importing polygonal sufaces saved in the STL format. The user indicates a single reference directory containing the files and all corresponding deformed directories. The STL files must be named the same in each directory, so they are matched appropriately. In most cases, these will have been generated by pyCellAnalyst.Volume(), and by default will be named [image_dirname]_results.

Parameters:
  • ref_dir (str) – The directory containing the STL files corresponding to the reference (undeformed) state.
  • def_dir (str) – The directory containing the STL files corresponding to the deformed state.
  • rigidInitial (bool, optional) – If True do an initial rigid body transformation to align objects.
  • deformable (bool, optional) – If True deformable image registration will be performed. This will call deformableRegistration(), which will calculate a displacement map between images reconstructed from reference and deformed surfaces.
  • saveFEA (bool, optional) – If True will save nodes, elements, surface nodes, and displacement boundary conditions in a dictionary to cell{:02d}.pkl. This information can then be used to run finite element analysis in whatever software the user desires.
  • deformableSettings (dict, optional) –
    Settings for deformable image registration with fields:
    • Iterations: (int, 200) The maximum number of iterations for algorithm.
    • Maximum RMS: (float, 0.01) Will terminate iterations if root-mean-square error is less than.
    • Displacement Smoothing: (float, 3.0) Variance for Gaussian smoothing of displacement field result.
    • Precision: (float, 0.01) The fraction of the object bounding box in each dimension spanned by 1 voxel.
  • display (bool, optional) – If True will display 3-D interactive rendering of displacement fields.
rsurfs

[,vtkPolyData,…] – Polygonal surfaces of objects in reference state.

dsurfs

[,vtkPolyData,…] – Polygonal surfaces of objects in deformed state.

rmeshes

[,vtkUnstructuredGrid,…] – Tetrahedral meshes of undeformed geometries.

rcentroids

[,ndarray(3, float),…] – Volumetric centroids of objects in reference state.

dcentroids

[,ndarray(3, float),…] – Volumetric centroids of objects in deformed state.

cell_strains

[,ndarray((3,3), float) – Green-Lagrange strain tensors for object assuming uniform deformation.

vstrains

[, float, …] – Volumetric strains of the analyzed objects.

ecm_strain

ndarray((3,3), float) – Green-Lagrange strain tensor for extracellular (extra-object) matrix assuming uniform deformation.

rvols

[,float,…] – Volumes of objects in reference state.

dvols

[,float,…] – Volumes of objects in deformed state.

raxes

[,ndarray(3, float),…] – Lengths of axes for ellipsoid with equivalent principal moments of inertia to object in reference state.

daxes

[,ndarray(3, float),…] – Lengths of axes for ellisoid with equivalent principal moments of inertia to object in deformed state.

cell_fields

[,vtkUnstructuredGrid,…] – Displacement vectors determined by deformable image registration interpolated to the vertices of object mesh in reference state.

rigid_transforms

[,vtkTransform,…] – Initial rigid body transforms applied to object in reference state to align with deformed state. This is only populated if rigidInitial is True.

_deform()

Calculates the affine transform that best maps the reference polygonal surface to its corresponding deformed surface. This transform is calculated through an interactive closest point optimization, that seeks to minimize the sum of distances between the reference surface vertices and the current affine transformed surface.

Assuming a uniform deformation, the non-translational elements of this affine transform compose the deformation gradient \mathbf{F}. The Green-Lagrange strain tensor is then defined as

\mathbf{E} = \frac{1}{2}(\mathbf{F}^T.\mathbf{F} - \mathbf{1}),

where \mathbf{1} is the identity.

Returns:
Return type:cell_strains
_getECMstrain()

Generates tetrahedrons from object centroids in the reference and deformed states. The highest quality tetrahedron (edge ratio closest to 1) is used to construct a linear system of equations,

\|\mathbf{w}\|^2 - \|\mathbf{W}\|^2 = \mathbf{W}.\mathbf{E}.\mathbf{W},

where, \mathbf{W} are the reference tetrahedron edges (as vectors) and \mathbf{w} are the deformed tetrahedron edges, to solve for Green-Lagrange strain, \mathbf{E}.

Returns:
Return type:ecm_strain
_make3Dmesh(filename, frame, vConst)

Generates a 3-D tetrahedral mesh using tetmesh module build on CGAL. These meshes are then used to determine the object’s volume, centroid, and the axes of the ellipsoid that has equivalent principal moments of inertia.

Parameters:
  • filename (str) – The path and filename of the STL surface currently being analyzed. This is necessary since MeshPy has to read the STL from disk in its native format.
  • frame (str) – Indicates the curent state, ‘MATERIAL’ or ‘SPATIAL’.
Returns:

  • **if frame==’MATERIAL’** –
    • rmeshes
    • rcentroids
    • rvols
    • raxes
  • **else** –
    • dcentroids
    • dvols
    • daxes
_poly2img(ind)

Helper function called by deformableRegistration() that generates images from polygonal surfaces in reference/deformed pairs. The voxel dimension of these images is determined by the value for Precision in deformableSettings.

Parameters:ind (int) – The list index for the current object pair being analyzed.
Returns:
Return type:(Reference Image, Deformed Image, Tranformed Reference Surface)
_readstls()

Reads in all STL files contained in directories indicated by ref_dir and def_dir. Also calls _make3Dmesh() to create 3-D tetrahedral meshes.

Returns:
Return type:rsurfs, dsurfs
animate(pd, ind)

Helper function called by deformableRegistration if animate is True. Spawns a window with an interactive 3-D rendering of the current analyzed object in its reference state. The displacements calculated from the deformable image registration can be applied to this object to animate the deformation by pressing the RIGHT-ARROW. Pressing the UP-ARROW will animate and also save the frames to disk.

Parameters:
  • pd (vtkPolyData) – The current analyzed object’s reference geometry.
  • ind (int) – The index of the current polydata in rsurfs. Necessary for naming directory created if animation frames are saved.
deformableRegistration()

Performs deformable image registration on images reconstructed from polygonal surfaces at a user-specified precision. If animate parameter is True, this also spawns a VTK interative rendering that can animate the deformation.

Returns:
Return type:cell_fields