import numpy
from scipy.spatial import ConvexHull
from radiomics.shape import RadiomicsShape
from .mybase import MyRadiomicsFeaturesBase


class MyRadiomicsShape(MyRadiomicsFeaturesBase):
    def __init__(self, inputImage, inputMask, **kwargs):
        super().__init__(inputImage, inputMask, **kwargs)

    def _initVoxelBasedCalculation(self):
        return RadiomicsShape._initVoxelBasedCalculation(self)

    def _initSegmentBasedCalculation(self):
        return RadiomicsShape._initSegmentBasedCalculation(self)

    def _calculateSurfaceArea(self):
        return RadiomicsShape._calculateSurfaceArea(self)

    def calculateDiameters(self):
        return RadiomicsShape.calculateDiameters(self)

    def _interpolate(self, grid, p1, p2):
        return RadiomicsShape._interpolate(self, grid, p1, p2)

    def getCompactness2FeatureValue(self):
        r"""
        **6. Compactness 2**

        .. math::
        \textit{compactness 2} = 36 \pi \frac{V^2}{A^3}

        Similar to Sphericity and Compactness 1, Compactness 2 is a measure of how compact the shape of the tumor is
        relative to a sphere (most compact). It is a dimensionless measure, independent of scale and orientation. The value
        range is :math:`0 < compactness\ 2 \leq 1`, where a value of 1 indicates a perfect sphere.

        By definition, :math:`compactness\ 2 = (sphericity)^3`

        .. note::
        This feature is correlated to Compactness 1, Sphericity and Spherical Disproportion.
        Therefore, this feature is marked, so it is not enabled by default (i.e. this feature will not be enabled if no
        individual features are specified (enabling 'all' features), but will be enabled when individual features are
        specified, including this feature). To include this feature in the extraction, specify it by name in the enabled features.
        """
        if self.SurfaceArea is None:
            self.SurfaceArea = self._calculateSurfaceArea()
        return (36.0 * numpy.pi) * (self.Volume ** 2.0) / (self.SurfaceArea ** 3.0)

    def getEccentricityFeatureValue(self):
        r"""
        **17. Eccentricity**

        This is a custom shape feature, implemented for comparison with Valieres et al.
        Eccentricity is defined as: 

        .. math::
        \textit{Eccentricity} = \sqrt{1-\frac{a\times b}{c^2}}

        Here, where c is the longest semi-principal axes of the ellipsoid,
        and a and b are the second and third longest semi-principal axes of the ellipsoid.
        """

        major_axis = numpy.sqrt(self.eigenValues[2]) * 4
        minor_axis = numpy.sqrt(self.eigenValues[1]) * 4
        least_axis = numpy.sqrt(self.eigenValues[0]) * 4

        return numpy.sqrt(1.0 - (minor_axis * least_axis / (major_axis ** 2)))

    def getSolidityFeatureValue(self):
        r"""
        **18. Solidity**

        This is a custom shape feature, implemented for comparison with Valieres et al.
        Eccentricity is defined as: 

        Solidity is defined as:
        
        Ratio of the number of voxels in the tumour region to the number of voxels in the 3D convex hull
        of the tumour region (smallest polyhedron containing the tumour region)

        """

        Np = len(self.labelledVoxelCoordinates[0])
        coordinates = numpy.array(self.labelledVoxelCoordinates, dtype='int').transpose(
            (1, 0)
        )  # Transpose equals zip(*a)
        physicalCoordinates = coordinates * self.pixelSpacing[None, :]
        physicalCoordinates -= numpy.mean(physicalCoordinates, axis=0)  # Centered at 0
        physicalCoordinates /= numpy.sqrt(Np)

        ch = ConvexHull(physicalCoordinates)

        simplices = numpy.column_stack(
            (numpy.repeat(ch.vertices[0], ch.nsimplex), ch.simplices)
        )
        tets = ch.points[simplices]
        a, b, c, d = tets[:, 0], tets[:, 1], tets[:, 2], tets[:, 3]
        tetrahedron_vol = (
            numpy.abs(numpy.einsum('ij,ij->i', a - d, numpy.cross(b - d, c - d))) / 6
        )

        ROIVolume = self.Volume

        return ROIVolume / numpy.sum(tetrahedron_vol)