--- a
+++ b/CLI/CancellousAnalysis/CancellousAnalysis.py
@@ -0,0 +1,155 @@
+#!/usr/bin/env python-real
+
+import sys
+import os
+from datetime import date
+
+import numpy as np
+
+sys.path.append(os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
+from MusculoskeletalAnalysisCLITools import (
+  crop,
+  bWshape,
+  densityMap,
+  findSpheres,
+  readImg,
+  readMask,
+  updateVertices,
+  writeReport,
+)
+
+
+OUTPUT_FIELDS = [
+    ("Date Analysis Performed", "The current date"),
+    ("Input Volume", "Name of the input volume"),
+    ("Total Volume (mm^3)", "The volume of the segmented area"),
+    ("Bone Volume (mm^3)", "The volume of cancellous bone in the segmented area, calculated using marching cubes"),
+    ("Bone Volume/Total Volume", "The fraction of the volume that is bone"),
+    ("Mean Trabecular Thickness (mm)", "The mean thickess of the bone, measured using largest sphere thickness, in milimeters"),
+    ("Trabecular Thickness Standard Deviation (mm)", "The standard deviation of the mean trabecular thickness"),
+    ("Mean Trabecular Spacing (mm)", "The mean thickness of the non area not containing bone in milimeters, measured using the same method as bone thickness."),
+    ("Trabecular Spacing Standard Deviation (mm)", "The standard deviation of the mean trabecular spacing"),
+    ("Trabecular Number", "Approximated as inverse of trabecular spacing"),
+    ("Structure Model Index", "A measurement of the trabecular shape. 0 is a plate, 3 is a rod, 4 is a sphere"),
+    ("Connectivity Density", "A measurement of the number of connections per volume, based on the Euler characteristic of the bone after removing isolated components and holes"),
+    ("Tissue Mineral Density(mgHA/cm^3)", "The mean density of the bone, measured in miligrams of hydroxyapatite per cubic centimeter"),
+    ("Voxel Dimension (mm)", "The side length of one voxel, measured in milimeters"),
+    ("Lower Threshold", "The lower threshold value for bone"),
+    ("Upper Threshold", "The upper threshold value for bone"),
+]
+
+
+# Performs analysis on cancellous bone
+# image: 3D image black and white of bone
+# mask: 3D labelmap of bone area, including cavities
+# threshold: The threshold for bone in the image. Any pixels with a higher value will be considered bone.
+# voxSize: The physical side length of the voxels, in mm
+# slope, intercept and scale: parameters for density conversion
+# output: The name of the output directory
+def main(inputImg, inputMask, lower, upper, voxSize, slope, intercept, name, output):
+    from skimage import measure
+
+    imgData = readImg(inputImg)
+    (_, maskData) = readMask(inputMask)
+    
+    (maskData, imgData) = crop(maskData, imgData)
+    if np.count_nonzero(maskData) == 0:
+         raise Exception("Segmentation mask is empty.")
+    trabecular = (imgData > lower) & (imgData <= upper)
+    # Create mesh using marching squares and calculate its volume
+    boneMesh = bWshape(imgData, lower)
+    boneVolume=boneMesh.volume * voxSize**3
+    # Find volume of entire mask by counting voxels
+    totalVolume = np.count_nonzero(maskData) * voxSize**3
+    bvtv = boneVolume/totalVolume
+    print("""<filter-progress>{}</filter-progress>""".format(.20))
+    sys.stdout.flush()
+    background = np.bitwise_and(maskData, np.invert(trabecular))
+    # Get the thickness map and calculate the thickness
+    rads = findSpheres(trabecular)
+    diams = rads * 2 * voxSize
+    thickness = np.mean(diams)
+    thicknessStd = np.std(diams)
+    print("""<filter-progress>{}</filter-progress>""".format(.40))
+    sys.stdout.flush()
+    # Get the thickness map for the background
+    rads = findSpheres(background)
+    diams = rads * 2 * voxSize
+    spacing = np.mean(diams)
+    spacingStd = np.std(diams)
+    trabecularNum = 1/spacing
+    print("""<filter-progress>{}</filter-progress>""".format(.60))
+    sys.stdout.flush()
+    # SMI
+    # Calculated by finding a 3d mesh, expanding the vertices by a small amount, and calculate a value based on the relative difference in surface area
+    # Measures the characteristics of the shape; 0 for plate-like, 3 for rod-like, 4 for sphere-like
+    # Surface area
+    bS = measure.mesh_surface_area(boneMesh.vertices, boneMesh.faces)*(voxSize**2)
+    # Arbitrary small value. Reducing size increases accuracy, but risks being rounded to zero
+    dr = 0.000001
+    # Creates a new mesh by moving each vertex a small distance in the direction of its vertex normal, then find the difference in surface area
+    newVert=np.add(boneMesh.vertices, boneMesh.vertex_normals*dr)
+    boneMesh = updateVertices(boneMesh, newVert)
+    dS = (measure.mesh_surface_area(boneMesh.vertices, boneMesh.faces)*(voxSize**2))-bS
+    # Convert to mm
+    dr = dr*voxSize
+    # Apply the SMI formula
+    SMI=(6*boneVolume*(dS/dr)/(bS**2))
+    print("""<filter-progress>{}</filter-progress>""".format(.80))
+    sys.stdout.flush()
+    # Calculate density
+    density = densityMap(imgData, slope, intercept)
+    tmd = np.mean(density[trabecular])
+
+    # Connnectivity
+    # Remove disconected holes and islands
+    labelmap=measure.label(np.invert(trabecular), connectivity=1)
+    largest=np.argmax(np.bincount(labelmap[np.nonzero(labelmap)]))
+    trabecular=labelmap!=largest
+    # Holes use face connectivity, switch to vertex connectivity for islands
+    labelmap=measure.label(trabecular, connectivity=3)
+    largest=np.argmax(np.bincount(labelmap[np.nonzero(labelmap)]))
+    trabecular=labelmap==largest
+
+    # Find the euler characteristic
+    # roughly equal to 2-2*number of holes
+    phi = measure.euler_number(trabecular, connectivity=3)
+    # connD gives an approximate measure of holes/connections per volume
+    connD = (1-phi)/totalVolume
+
+
+
+
+    fPath = os.path.join(output, "cancellous.txt")
+
+    header = [field[0] for field in OUTPUT_FIELDS]
+
+    data = [
+        date.today(),
+        name,
+        totalVolume,
+        boneVolume,
+        bvtv,
+        thickness,
+        thicknessStd,
+        spacing,
+        spacingStd,
+        trabecularNum,
+        SMI,
+        connD,
+        tmd,
+        voxSize,
+        lower,
+        upper
+    ]
+
+    writeReport(fPath, header, data)
+
+if __name__ == "__main__":
+    if len(sys.argv) < 10:
+        print(sys.argv)
+        print(len(sys.argv))
+        print("Usage: CancellousAnalysis <input> <mask> <lowerThreshold> <upperThreshold> <voxelSize> <slope> <intercept> <name> <output>")
+        sys.exit(1)
+
+    main(sys.argv[1], sys.argv[2], float(sys.argv[3]), float(sys.argv[4]), float(sys.argv[5]), float(sys.argv[6]), float(sys.argv[7]), str(sys.argv[8]), str(sys.argv[9]))