#!/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]))