from __future__ import annotations

import numpy as np

from radiomics import base, cMatrices


class RadiomicsGLSZM(base.RadiomicsFeaturesBase):
    r"""
  A Gray Level Size Zone (GLSZM) quantifies gray level zones in an image. A gray level zone is defined as a the number
  of connected voxels that share the same gray level intensity. A voxel is considered connected if the distance is 1
  according to the infinity norm (26-connected region in a 3D, 8-connected region in 2D).
  In a gray level size zone matrix :math:`P(i,j)` the :math:`(i,j)^{\text{th}}` element equals the number of zones
  with gray level :math:`i` and size :math:`j` appear in image. Contrary to GLCM and GLRLM, the GLSZM is rotation
  independent, with only one matrix calculated for all directions in the ROI.

  As a two dimensional example, consider the following 5x5 image, with 5 discrete gray levels:

  .. math::
    \textbf{I} = \begin{bmatrix}
    5 & 2 & 5 & 4 & 4\\
    3 & 3 & 3 & 1 & 3\\
    2 & 1 & 1 & 1 & 3\\
    4 & 2 & 2 & 2 & 3\\
    3 & 5 & 3 & 3 & 2 \end{bmatrix}

  The GLSZM then becomes:

  .. math::
    \textbf{P} = \begin{bmatrix}
    0 & 0 & 0 & 1 & 0\\
    1 & 0 & 0 & 0 & 1\\
    1 & 0 & 1 & 0 & 1\\
    1 & 1 & 0 & 0 & 0\\
    3 & 0 & 0 & 0 & 0 \end{bmatrix}

  Let:

  - :math:`N_g` be the number of discrete intensity values in the image
  - :math:`N_s` be the number of discrete zone sizes in the image
  - :math:`N_p` be the number of voxels in the image
  - :math:`N_z` be the number of zones in the ROI, which is equal to :math:`\sum^{N_g}_{i=1}\sum^{N_s}_{j=1}
    {\textbf{P}(i,j)}` and :math:`1 \leq N_z \leq N_p`
  - :math:`\textbf{P}(i,j)` be the size zone matrix
  - :math:`p(i,j)` be the normalized size zone matrix, defined as :math:`p(i,j) = \frac{\textbf{P}(i,j)}{N_z}`

  .. note::
    The mathematical formulas that define the GLSZM features correspond to the definitions of features extracted from
    the GLRLM.

  References

  - Guillaume Thibault; Bernard Fertil; Claire Navarro; Sandrine Pereira; Pierre Cau; Nicolas Levy; Jean Sequeira;
    Jean-Luc Mari (2009). "Texture Indexes and Gray Level Size Zone Matrix. Application to Cell Nuclei Classification".
    Pattern Recognition and Information Processing (PRIP): 140-145.
  - `<https://en.wikipedia.org/wiki/Gray_level_size_zone_matrix>`_
  """

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

        self.P_glszm = None
        self.imageArray = self._applyBinning(self.imageArray)

    def _initCalculation(self, voxelCoordinates=None):
        self.P_glszm = self._calculateMatrix(voxelCoordinates)

        self._calculateCoefficients()

        self.logger.debug(
            "GLSZM feature class initialized, calculated GLSZM with shape %s",
            self.P_glszm.shape,
        )

    def _calculateMatrix(self, voxelCoordinates=None):
        """
        Number of times a region with a
        gray level and voxel count occurs in an image. P_glszm[level, voxel_count] = # occurrences

        For 3D-images this concerns a 26-connected region, for 2D an 8-connected region
        """
        self.logger.debug("Calculating GLSZM matrix in C")
        Ng = self.coefficients["Ng"]
        Ns = np.sum(self.maskArray)

        matrix_args = [
            self.imageArray,
            self.maskArray,
            Ng,
            Ns,
            self.settings.get("force2D", False),
            self.settings.get("force2Ddimension", 0),
        ]
        if self.voxelBased:
            matrix_args += [self.settings.get("kernelRadius", 1), voxelCoordinates]

        P_glszm = cMatrices.calculate_glszm(*matrix_args)  # shape (Nvox, Ng, Ns)

        # Delete rows that specify gray levels not present in the ROI
        NgVector = range(1, Ng + 1)  # All possible gray values
        GrayLevels = self.coefficients["grayLevels"]  # Gray values present in ROI
        emptyGrayLevels = np.array(
            list(set(NgVector) - set(GrayLevels)), dtype=int
        )  # Gray values NOT present in ROI

        return np.delete(P_glszm, emptyGrayLevels - 1, 1)

    def _calculateCoefficients(self):
        self.logger.debug("Calculating GLSZM coefficients")

        ps = np.sum(self.P_glszm, 1)  # shape (Nvox, Ns)
        pg = np.sum(self.P_glszm, 2)  # shape (Nvox, Ng)

        ivector = self.coefficients["grayLevels"].astype(float)  # shape (Ng,)
        jvector = np.arange(
            1, self.P_glszm.shape[2] + 1, dtype=np.float64
        )  # shape (Ns,)

        # Get the number of zones in this GLSZM
        Nz = np.sum(self.P_glszm, (1, 2))  # shape (Nvox,)
        Nz[Nz == 0] = 1  # set sum to numpy.spacing(1) if sum is 0?

        # Get the number of voxels represented by this GLSZM: Multiply the zones by their size and sum them
        Np = np.sum(ps * jvector[None, :], 1)  # shape (Nvox, )
        Np[Np == 0] = 1

        # Delete columns that specify zone sizes not present in the ROI
        emptyZoneSizes = np.where(np.sum(ps, 0) == 0)
        self.P_glszm = np.delete(self.P_glszm, emptyZoneSizes, 2)
        jvector = np.delete(jvector, emptyZoneSizes)
        ps = np.delete(ps, emptyZoneSizes, 1)

        self.coefficients["Np"] = Np
        self.coefficients["Nz"] = Nz
        self.coefficients["ps"] = ps
        self.coefficients["pg"] = pg
        self.coefficients["ivector"] = ivector
        self.coefficients["jvector"] = jvector

    def getSmallAreaEmphasisFeatureValue(self):
        r"""
        **1. Small Area Emphasis (SAE)**

        .. math::
          \textit{SAE} = \frac{\sum^{N_g}_{i=1}\sum^{N_s}_{j=1}{\frac{\textbf{P}(i,j)}{j^2}}}{N_z}

        SAE is a measure of the distribution of small size zones, with a greater value indicative of more smaller size zones
        and more fine textures.
        """
        ps = self.coefficients["ps"]
        jvector = self.coefficients["jvector"]
        Nz = self.coefficients["Nz"]

        return np.sum(ps / (jvector[None, :] ** 2), 1) / Nz

    def getLargeAreaEmphasisFeatureValue(self):
        r"""
        **2. Large Area Emphasis (LAE)**

        .. math::
          \textit{LAE} = \frac{\sum^{N_g}_{i=1}\sum^{N_s}_{j=1}{\textbf{P}(i,j)j^2}}{N_z}

        LAE is a measure of the distribution of large area size zones, with a greater value indicative of more larger size
        zones and more coarse textures.
        """
        ps = self.coefficients["ps"]
        jvector = self.coefficients["jvector"]
        Nz = self.coefficients["Nz"]

        return np.sum(ps * (jvector[None, :] ** 2), 1) / Nz

    def getGrayLevelNonUniformityFeatureValue(self):
        r"""
        **3. Gray Level Non-Uniformity (GLN)**

        .. math::
          \textit{GLN} = \frac{\sum^{N_g}_{i=1}\left(\sum^{N_s}_{j=1}{\textbf{P}(i,j)}\right)^2}{N_z}

        GLN measures the variability of gray-level intensity values in the image, with a lower value indicating more
        homogeneity in intensity values.
        """
        pg = self.coefficients["pg"]
        Nz = self.coefficients["Nz"]

        return np.sum(pg**2, 1) / Nz

    def getGrayLevelNonUniformityNormalizedFeatureValue(self):
        r"""
        **4. Gray Level Non-Uniformity Normalized (GLNN)**

        .. math::
          \textit{GLNN} = \frac{\sum^{N_g}_{i=1}\left(\sum^{N_s}_{j=1}{\textbf{P}(i,j)}\right)^2}{N_z^2}

        GLNN measures the variability of gray-level intensity values in the image, with a lower value indicating a greater
        similarity in intensity values. This is the normalized version of the GLN formula.
        """
        pg = self.coefficients["pg"]
        Nz = self.coefficients["Nz"]

        return np.sum(pg**2, 1) / Nz**2

    def getSizeZoneNonUniformityFeatureValue(self):
        r"""
        **5. Size-Zone Non-Uniformity (SZN)**

        .. math::
          \textit{SZN} = \frac{\sum^{N_s}_{j=1}\left(\sum^{N_g}_{i=1}{\textbf{P}(i,j)}\right)^2}{N_z}

        SZN measures the variability of size zone volumes in the image, with a lower value indicating more homogeneity in
        size zone volumes.
        """
        ps = self.coefficients["ps"]
        Nz = self.coefficients["Nz"]

        return np.sum(ps**2, 1) / Nz

    def getSizeZoneNonUniformityNormalizedFeatureValue(self):
        r"""
        **6. Size-Zone Non-Uniformity Normalized (SZNN)**

        .. math::
          \textit{SZNN} = \frac{\sum^{N_s}_{j=1}\left(\sum^{N_g}_{i=1}{\textbf{P}(i,j)}\right)^2}{N_z^2}

        SZNN measures the variability of size zone volumes throughout the image, with a lower value indicating more
        homogeneity among zone size volumes in the image. This is the normalized version of the SZN formula.
        """
        ps = self.coefficients["ps"]
        Nz = self.coefficients["Nz"]

        return np.sum(ps**2, 1) / Nz**2

    def getZonePercentageFeatureValue(self):
        r"""
        **7. Zone Percentage (ZP)**

        .. math::
          \textit{ZP} = \frac{N_z}{N_p}

        ZP measures the coarseness of the texture by taking the ratio of number of zones and number of voxels in the ROI.

        Values are in range :math:`\frac{1}{N_p} \leq ZP \leq 1`, with higher values indicating a larger portion of the ROI
        consists of small zones (indicates a more fine texture).
        """
        Nz = self.coefficients["Nz"]
        Np = self.coefficients["Np"]

        return Nz / Np

    def getGrayLevelVarianceFeatureValue(self):
        r"""
        **8. Gray Level Variance (GLV)**

        .. math::
          \textit{GLV} = \displaystyle\sum^{N_g}_{i=1}\displaystyle\sum^{N_s}_{j=1}{p(i,j)(i - \mu)^2}

        Here, :math:`\mu = \displaystyle\sum^{N_g}_{i=1}\displaystyle\sum^{N_s}_{j=1}{p(i,j)i}`

        GLV measures the variance in gray level intensities for the zones.
        """
        ivector = self.coefficients["ivector"]
        Nz = self.coefficients["Nz"]
        pg = (
            self.coefficients["pg"] / Nz[:, None]
        )  # divide by Nz to get the normalized matrix

        u_i = np.sum(pg * ivector[None, :], 1, keepdims=True)
        return np.sum(pg * (ivector[None, :] - u_i) ** 2, 1)

    def getZoneVarianceFeatureValue(self):
        r"""
        **9. Zone Variance (ZV)**

        .. math::
          \textit{ZV} = \displaystyle\sum^{N_g}_{i=1}\displaystyle\sum^{N_s}_{j=1}{p(i,j)(j - \mu)^2}

        Here, :math:`\mu = \displaystyle\sum^{N_g}_{i=1}\displaystyle\sum^{N_s}_{j=1}{p(i,j)j}`

        ZV measures the variance in zone size volumes for the zones.
        """
        jvector = self.coefficients["jvector"]
        Nz = self.coefficients["Nz"]
        ps = (
            self.coefficients["ps"] / Nz[:, None]
        )  # divide by Nz to get the normalized matrix

        u_j = np.sum(ps * jvector[None, :], 1, keepdims=True)
        return np.sum(ps * (jvector[None, :] - u_j) ** 2, 1)

    def getZoneEntropyFeatureValue(self):
        r"""
        **10. Zone Entropy (ZE)**

        .. math::
          \textit{ZE} = -\displaystyle\sum^{N_g}_{i=1}\displaystyle\sum^{N_s}_{j=1}{p(i,j)\log_{2}(p(i,j)+\epsilon)}

        Here, :math:`\epsilon` is an arbitrarily small positive number (:math:`\approx 2.2\times10^{-16}`).

        ZE measures the uncertainty/randomness in the distribution of zone sizes and gray levels. A higher value indicates
        more heterogeneity in the texture patterns.
        """
        eps = np.spacing(1)
        Nz = self.coefficients["Nz"]
        p_glszm = (
            self.P_glszm / Nz[:, None, None]
        )  # divide by Nz to get the normalized matrix

        return -np.sum(p_glszm * np.log2(p_glszm + eps), (1, 2))

    def getLowGrayLevelZoneEmphasisFeatureValue(self):
        r"""
        **11. Low Gray Level Zone Emphasis (LGLZE)**

        .. math::
          \textit{LGLZE} = \frac{\sum^{N_g}_{i=1}\sum^{N_s}_{j=1}{\frac{\textbf{P}(i,j)}{i^2}}}{N_z}

        LGLZE measures the distribution of lower gray-level size zones, with a higher value indicating a greater proportion
        of lower gray-level values and size zones in the image.
        """
        pg = self.coefficients["pg"]
        ivector = self.coefficients["ivector"]
        Nz = self.coefficients["Nz"]

        return np.sum(pg / (ivector[None, :] ** 2), 1) / Nz

    def getHighGrayLevelZoneEmphasisFeatureValue(self):
        r"""
        **12. High Gray Level Zone Emphasis (HGLZE)**

        .. math::
          \textit{HGLZE} = \frac{\sum^{N_g}_{i=1}\sum^{N_s}_{j=1}{\textbf{P}(i,j)i^2}}{N_z}

        HGLZE measures the distribution of the higher gray-level values, with a higher value indicating a greater proportion
        of higher gray-level values and size zones in the image.
        """
        pg = self.coefficients["pg"]
        ivector = self.coefficients["ivector"]
        Nz = self.coefficients["Nz"]

        return np.sum(pg * (ivector[None, :] ** 2), 1) / Nz

    def getSmallAreaLowGrayLevelEmphasisFeatureValue(self):
        r"""
        **13. Small Area Low Gray Level Emphasis (SALGLE)**

        .. math::
          \textit{SALGLE} = \frac{\sum^{N_g}_{i=1}\sum^{N_s}_{j=1}{\frac{\textbf{P}(i,j)}{i^2j^2}}}{N_z}

        SALGLE measures the proportion in the image of the joint distribution of smaller size zones with lower gray-level
        values.
        """
        ivector = self.coefficients["ivector"]
        jvector = self.coefficients["jvector"]
        Nz = self.coefficients["Nz"]

        return (
            np.sum(
                self.P_glszm
                / ((ivector[None, :, None] ** 2) * (jvector[None, None, :] ** 2)),
                (1, 2),
            )
            / Nz
        )

    def getSmallAreaHighGrayLevelEmphasisFeatureValue(self):
        r"""
        **14. Small Area High Gray Level Emphasis (SAHGLE)**

        .. math::
          \textit{SAHGLE} = \frac{\sum^{N_g}_{i=1}\sum^{N_s}_{j=1}{\frac{\textbf{P}(i,j)i^2}{j^2}}}{N_z}

        SAHGLE measures the proportion in the image of the joint distribution of smaller size zones with higher gray-level
        values.
        """
        ivector = self.coefficients["ivector"]
        jvector = self.coefficients["jvector"]
        Nz = self.coefficients["Nz"]

        return (
            np.sum(
                self.P_glszm
                * (ivector[None, :, None] ** 2)
                / (jvector[None, None, :] ** 2),
                (1, 2),
            )
            / Nz
        )

    def getLargeAreaLowGrayLevelEmphasisFeatureValue(self):
        r"""
        **15. Large Area Low Gray Level Emphasis (LALGLE)**

        .. math::
          \textit{LALGLE} = \frac{\sum^{N_g}_{i=1}\sum^{N_s}_{j=1}{\frac{\textbf{P}(i,j)j^2}{i^2}}}{N_z}

        LALGLE measures the proportion in the image of the joint distribution of larger size zones with lower gray-level
        values.
        """
        ivector = self.coefficients["ivector"]
        jvector = self.coefficients["jvector"]
        Nz = self.coefficients["Nz"]

        return (
            np.sum(
                self.P_glszm
                * (jvector[None, None, :] ** 2)
                / (ivector[None, :, None] ** 2),
                (1, 2),
            )
            / Nz
        )

    def getLargeAreaHighGrayLevelEmphasisFeatureValue(self):
        r"""
        **16. Large Area High Gray Level Emphasis (LAHGLE)**

        .. math::
          \textit{LAHGLE} = \frac{\sum^{N_g}_{i=1}\sum^{N_s}_{j=1}{\textbf{P}(i,j)i^2j^2}}{N_z}

        LAHGLE measures the proportion in the image of the joint distribution of larger size zones with higher gray-level
        values.
        """
        ivector = self.coefficients["ivector"]
        jvector = self.coefficients["jvector"]
        Nz = self.coefficients["Nz"]

        return (
            np.sum(
                self.P_glszm
                * (ivector[None, :, None] ** 2)
                * (jvector[None, None, :] ** 2),
                (1, 2),
            )
            / Nz
        )
