import numpy
from oct2py import Oct2Py, Oct2PyError
from scipy import ndimage

from radiomics.firstorder import RadiomicsFirstOrder

from .mybase import MyRadiomicsFeaturesBase


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

        self.pixelSpacing = inputImage.GetSpacing()
        self.voxelArrayShift = kwargs.get('voxelArrayShift', 0)

    def _initCalculation(self):
        return super()._initSegmentBasedCalculation()

        self.targetVoxelArray = self.imageArray[self.labelledVoxelCoordinates].astype(
            'float'
        )
        self.discretizedTargetVoxelArray = None  # Lazy instantiation

        self.logger.debug('First order feature class initialized')

    def _moment(self, a, moment=1, axis=0):
        return RadiomicsFirstOrder._moment(self, a, moment, axis)

    def _getDiscretizedTargetVoxelArray(self):
        return RadiomicsFirstOrder._getDiscretizedTargetVoxelArray(self)

    def getSuvMaxFeatureValue(self):
        r"""
        **20. SUV max**

        Implemented as one of the intensity features extracted by Vallieres et al.
        Maximum SUV of the tumour region. 

        .. note::
        Extracted from PET scans and not used in the CT feature set.
        """
        ROIPet = self.imageArray
        mask = self.maskArray
        ROIPet[~mask] = numpy.nan

        return numpy.max(ROIPet[~numpy.isnan(ROIPet)])

    def getSuvMeanFeatureValue(self):
        r"""
        **21. SUV mean**

        Implemented as one of the intensity features extracted by valieres et al.
        Average SUV of the tumour region. 

        .. note::
        Extracted from PET scans and not used in the CT feature set.
        """
        ROIPet = self.imageArray
        mask = self.maskArray
        ROIPet[~mask] = numpy.nan

        return numpy.mean(ROIPet[~numpy.isnan(ROIPet)])

    def getTLGFeatureValue(self):
        r"""
        **22. TLG**

        Total lesion glycolysis. 
        Defined as 
        
        .. math::
        \textit{SUVMean}\times \textit{total volume of the tumour region}
    
        .. note::
        Extracted from PET scans and not used in the CT feature set.
        """
        z, y, x = self.pixelSpacing
        Np = len(self.labelledVoxelCoordinates[0])
        volume = Np * (z * x * y)

        ROIPet = self.imageArray
        mask = self.maskArray
        ROIPet[~mask] = numpy.nan

        return numpy.mean(ROIPet[~numpy.isnan(ROIPet)]) * volume

    def getAUCCSHFeatureValue(self):
        r"""
        **23. Inactive volume **

        Area under the curve of the cumulative SUV-volume histogram describing the percentage 
        of total tumour volume above a percentage threshold of maximum SUV, as defined by van Velden et al.

        .. note::
        Extracted from PET scans and not used in the CT feature set.
        """
        nBins = 1000  # By default
        ROIPet = self.imageArray
        mask = self.maskArray

        ROIPet[~mask] = numpy.nan

        outliers = numpy.where(
            ROIPet
            > (
                numpy.mean(ROIPet[~numpy.isnan(ROIPet)])
                + 3 * numpy.std(ROIPet[~numpy.isnan(ROIPet)])
            )
        )[0]

        good_voxels = numpy.where(
            ROIPet
            <= (
                numpy.mean(ROIPet[~numpy.isnan(ROIPet)])
                + 3 * numpy.std(ROIPet[~numpy.isnan(ROIPet)])
            )
        )[0]

        ROIPet[outliers] = numpy.mean(ROIPet[good_voxels])
        ROIPet = ROIPet - numpy.min(ROIPet[~numpy.isnan(ROIPet)])
        ROIPet = ROIPet / numpy.max(ROIPet[~numpy.isnan(ROIPet)])
        volume = ROIPet[~numpy.isnan(ROIPet)]
        nVoxel = len(volume)

        bins, _ = numpy.histogram(volume, nBins)

        csh = numpy.flipud(numpy.cumsum(numpy.flipud(bins)) / nVoxel)
        return numpy.sum(csh / nBins)

    def getgETUFeatureValue(self):
        r"""
        **24. gETU**

        Generalized effective total uptake, with parameter a = 0.25 as defined by Rahim et al.

        - input: 3D array representing the PET volume in SUV format, with 
                voxels outside the ROI set to NaNs.

        .. note::
        Extracted from PET scans and not used in the CT feature set.
        """
        ROIPet = self.imageArray
        mask = self.maskArray
        ROIPet[~mask] = numpy.nan

        a = 0.25
        n_voxels = numpy.sum(ROIPet[~numpy.isnan(ROIPet)])
        ROIPet = ROIPet ** a

        return numpy.sum(ROIPet[~numpy.isnan(ROIPet)] / n_voxels) ** (1 / a)

    def getSuvPeakFeatureValue(self):
        r"""
        **25. SuvPeak**

        Average of the voxel with maximum SUV within the tumour region and its 26 connected neighbours. 
    
        - input: 3D array representing the PET volume in SUV format
    
        .. note::
        Extracted from PET scans and not used in the CT feature set.
        """
        ROIPet = self.imageArray
        mask = self.maskArray

        ROIPet[~mask] = numpy.nan

        oc = Oct2Py()

        with Oct2Py() as oc:
            try:
                oc.eval('pkg load image')
                ROIPet = oc.padarray(
                    ROIPet, oc.double(numpy.array([1, 1, 1]).tolist()), numpy.nan
                )

                SUVmax = numpy.max(ROIPet[~numpy.isnan(ROIPet)])
                indMax = numpy.where(ROIPet == SUVmax)

                indMaxX, indMaxY, indMaxZ = indMax[0][0], indMax[1][0], indMax[2][0]

                connectivity = numpy.array(
                    [
                        [-1, -1, -1],
                        [0, -1, -1],
                        [1, -1, -1],
                        [-1, 0, -1],
                        [0, 0, -1],
                        [1, 0, -1],
                        [-1, 1, -1],
                        [0, 1, -1],
                        [1, 1, -1],
                        [-1, -1, 0],
                        [0, -1, 0],
                        [1, -1, 0],
                        [-1, 0, 0],
                        [0, 0, 0],
                        [1, 0, 0],
                        [-1, 1, 0],
                        [0, 1, 0],
                        [1, 1, 0],
                        [-1, -1, 1],
                        [0, -1, 1],
                        [1, -1, 1],
                        [-1, 0, 1],
                        [0, 0, 1],
                        [1, 0, 1],
                        [-1, 1, 1],
                        [0, 1, 1],
                        [1, 1, 1],
                    ]
                )
                nPeak = len(connectivity)
                neighborsMax = numpy.zeros((1, nPeak))

                for i in range(nPeak):
                    neighborsMax[0, i] = ROIPet[
                        connectivity[i, 0] + indMaxX - 1,
                        connectivity[i, 1] + indMaxY - 1,
                        connectivity[i, 2] + indMaxZ - 1,
                    ]

            except Oct2PyError as e:
                print(e)
                oc.exit()

        return numpy.mean(neighborsMax[~numpy.isnan(neighborsMax)])

    def getInactiveVolumeFeatureValue(self):
        r"""
        **26. Inactive volume **

        Percentage of the tumour region that is inactive. Ast defined by Valieres et al
        a threshold of 0.01 × (SUVmax)^2 followed by closing and opening morphological operations were used 
        to differentiate active and inactive regions on FDG-PET scans
        
        - input: 3D array representing the PET volume in SUV format

        .. note::
        Extracted from PET scans and not used in the CT feature set.
        """

        ROIPet = self.imageArray
        mask = self.maskArray

        thresh = 0.01 * (numpy.max(ROIPet[~numpy.isnan(ROIPet)]) ** 2)

        mask_thresholded = ROIPet > thresh
        mask_inactive = numpy.logical_and(mask_thresholded, mask)

        # MORPHOLOGICAL OPERATIONS
        conn = numpy.zeros([5, 5, 5])
        conn[:, :, 0] = numpy.array(
            [
                [0, 0, 0, 0, 0],
                [0, 0, 0, 0, 0],
                [0, 0, 1, 0, 0],
                [0, 0, 0, 0, 0],
                [0, 0, 0, 0, 0],
            ]
        )
        conn[:, :, 4] = conn[:, :, 0]

        conn[:, :, 1] = numpy.array(
            [
                [0, 0, 0, 0, 0],
                [0, 0, 1, 0, 0],
                [0, 1, 1, 1, 0],
                [0, 0, 1, 0, 0],
                [0, 0, 0, 0, 0],
            ]
        )
        conn[:, :, 3] = conn[:, :, 1]

        conn[:, :, 2] = numpy.array(
            [
                [0, 0, 1, 0, 0],
                [0, 1, 1, 1, 0],
                [1, 1, 1, 1, 1],
                [0, 1, 1, 1, 0],
                [0, 0, 1, 0, 0],
            ]
        )

        with Oct2Py() as oc:
            try:
                oc.eval('pkg load image')
                perimeter = oc.bwperim(mask, 26)
                mask_inactive = oc.imclose(mask_inactive, conn)
                mask_inactive = oc.imopen(mask_inactive, conn)
                new_mask = mask_inactive + perimeter

                new_mask[mask == 0] = 10
                new_mask[new_mask == 1] = 10
                new_mask[new_mask == 2] = 10
                new_mask[new_mask == 0] = 1
                new_mask[new_mask == 10] = 0

                connObjects = oc.bwconncomp(new_mask, 26)

                if int(connObjects.NumObjects) > 1:

                    b = numpy.zeros((1, int(connObjects.NumObjects)))

                    for i in range(int(connObjects.NumObjects)):
                        if isinstance(connObjects.PixelIdxList[0][i], numpy.ndarray):
                            a = connObjects.PixelIdxList[0][i].shape[0]
                            # If the number of of pixel forming and object is lower than 15, reject it
                            if a < 15:
                                b[:, i] = 0
                            else:
                                b[:, i] = 1
                        else:
                            b[:, i] = 0

                    rows, cols = numpy.where(b > 0)
                    sumInactive = sum(
                        [
                            connObjects.PixelIdxList[rows, cols][i].shape[0]
                            for i in range(len(cols))
                        ]
                    )

                else:
                    sumInactive = 0
                sumVolume = numpy.sum(mask)

            except Oct2PyError as e:
                print(e)
                oc.exit()

        return sumInactive / sumVolume * 100