tritools.py

"""
Tools for triangular grids.
"""

import numpy as np

from matplotlib import cbook
from matplotlib.tri import Triangulation


class TriAnalyzer:
    """
    Define basic tools for triangular mesh analysis and improvement.

    A TriAnalyzer encapsulates a :class:`~matplotlib.tri.Triangulation`
    object and provides basic tools for mesh analysis and mesh improvement.

    Parameters
    ----------
    triangulation : :class:`~matplotlib.tri.Triangulation` object
        The encapsulated triangulation to analyze.

    Attributes
    ----------
    `scale_factors`

    """
    def __init__(self, triangulation):
        cbook._check_isinstance(Triangulation, triangulation=triangulation)
        self._triangulation = triangulation

    @property
    def scale_factors(self):
        """
        Factors to rescale the triangulation into a unit square.

        Returns *k*, tuple of 2 scale factors.

        Returns
        -------
        k : tuple of 2 floats (kx, ky)
            Tuple of floats that would rescale the triangulation :
            ``[triangulation.x * kx, triangulation.y * ky]``
            fits exactly inside a unit square.

        """
        compressed_triangles = self._triangulation.get_masked_triangles()
        node_used = (np.bincount(np.ravel(compressed_triangles),
                                 minlength=self._triangulation.x.size) != 0)
        return (1 / np.ptp(self._triangulation.x[node_used]),
                1 / np.ptp(self._triangulation.y[node_used]))

    def circle_ratios(self, rescale=True):
        """
        Returns a measure of the triangulation triangles flatness.

        The ratio of the incircle radius over the circumcircle radius is a
        widely used indicator of a triangle flatness.
        It is always ``<= 0.5`` and ``== 0.5`` only for equilateral
        triangles. Circle ratios below 0.01 denote very flat triangles.

        To avoid unduly low values due to a difference of scale between the 2
        axis, the triangular mesh can first be rescaled to fit inside a unit
        square with :attr:`scale_factors` (Only if *rescale* is True, which is
        its default value).

        Parameters
        ----------
        rescale : boolean, optional
            If True, a rescaling will be internally performed (based on
            :attr:`scale_factors`, so that the (unmasked) triangles fit
            exactly inside a unit square mesh. Default is True.

        Returns
        -------
        circle_ratios : masked array
            Ratio of the incircle radius over the
            circumcircle radius, for each 'rescaled' triangle of the
            encapsulated triangulation.
            Values corresponding to masked triangles are masked out.

        """
        # Coords rescaling
        if rescale:
            (kx, ky) = self.scale_factors
        else:
            (kx, ky) = (1.0, 1.0)
        pts = np.vstack([self._triangulation.x*kx,
                         self._triangulation.y*ky]).T
        tri_pts = pts[self._triangulation.triangles]
        # Computes the 3 side lengths
        a = tri_pts[:, 1, :] - tri_pts[:, 0, :]
        b = tri_pts[:, 2, :] - tri_pts[:, 1, :]
        c = tri_pts[:, 0, :] - tri_pts[:, 2, :]
        a = np.hypot(a[:, 0], a[:, 1])
        b = np.hypot(b[:, 0], b[:, 1])
        c = np.hypot(c[:, 0], c[:, 1])
        # circumcircle and incircle radii
        s = (a+b+c)*0.5
        prod = s*(a+b-s)*(a+c-s)*(b+c-s)
        # We have to deal with flat triangles with infinite circum_radius
        bool_flat = (prod == 0.)
        if np.any(bool_flat):
            # Pathologic flow
            ntri = tri_pts.shape[0]
            circum_radius = np.empty(ntri, dtype=np.float64)
            circum_radius[bool_flat] = np.inf
            abc = a*b*c
            circum_radius[~bool_flat] = abc[~bool_flat] / (
                4.0*np.sqrt(prod[~bool_flat]))
        else:
            # Normal optimized flow
            circum_radius = (a*b*c) / (4.0*np.sqrt(prod))
        in_radius = (a*b*c) / (4.0*circum_radius*s)
        circle_ratio = in_radius/circum_radius
        mask = self._triangulation.mask
        if mask is None:
            return circle_ratio
        else:
            return np.ma.array(circle_ratio, mask=mask)

    def get_flat_tri_mask(self, min_circle_ratio=0.01, rescale=True):
        """
        Eliminates excessively flat border triangles from the triangulation.

        Returns a mask *new_mask* which allows to clean the encapsulated
        triangulation from its border-located flat triangles
        (according to their :meth:`circle_ratios`).
        This mask is meant to be subsequently applied to the triangulation
        using :func:`matplotlib.tri.Triangulation.set_mask`.
        *new_mask* is an extension of the initial triangulation mask
        in the sense that an initially masked triangle will remain masked.

        The *new_mask* array is computed recursively; at each step flat
        triangles are removed only if they share a side with the current mesh
        border. Thus no new holes in the triangulated domain will be created.

        Parameters
        ----------
        min_circle_ratio : float, optional
            Border triangles with incircle/circumcircle radii ratio r/R will
            be removed if r/R < *min_circle_ratio*. Default value: 0.01
        rescale : boolean, optional
            If True, a rescaling will first be internally performed (based on
            :attr:`scale_factors` ), so that the (unmasked) triangles fit
            exactly inside a unit square mesh. This rescaling accounts for the
            difference of scale which might exist between the 2 axis. Default
            (and recommended) value is True.

        Returns
        -------
        new_mask : array-like of booleans
            Mask to apply to encapsulated triangulation.
            All the initially masked triangles remain masked in the
            *new_mask*.

        Notes
        -----
        The rationale behind this function is that a Delaunay
        triangulation - of an unstructured set of points - sometimes contains
        almost flat triangles at its border, leading to artifacts in plots
        (especially for high-resolution contouring).
        Masked with computed *new_mask*, the encapsulated
        triangulation would contain no more unmasked border triangles
        with a circle ratio below *min_circle_ratio*, thus improving the
        mesh quality for subsequent plots or interpolation.
        """
        # Recursively computes the mask_current_borders, true if a triangle is
        # at the border of the mesh OR touching the border through a chain of
        # invalid aspect ratio masked_triangles.
        ntri = self._triangulation.triangles.shape[0]
        mask_bad_ratio = self.circle_ratios(rescale) < min_circle_ratio

        current_mask = self._triangulation.mask
        if current_mask is None:
            current_mask = np.zeros(ntri, dtype=bool)
        valid_neighbors = np.copy(self._triangulation.neighbors)
        renum_neighbors = np.arange(ntri, dtype=np.int32)
        nadd = -1
        while nadd != 0:
            # The active wavefront is the triangles from the border (unmasked
            # but with a least 1 neighbor equal to -1
            wavefront = (np.min(valid_neighbors, axis=1) == -1) & ~current_mask
            # The element from the active wavefront will be masked if their
            # circle ratio is bad.
            added_mask = wavefront & mask_bad_ratio
            current_mask = added_mask | current_mask
            nadd = np.sum(added_mask)

            # now we have to update the tables valid_neighbors
            valid_neighbors[added_mask, :] = -1
            renum_neighbors[added_mask] = -1
            valid_neighbors = np.where(valid_neighbors == -1, -1,
                                       renum_neighbors[valid_neighbors])

        return np.ma.filled(current_mask, True)

    def _get_compressed_triangulation(self, return_tri_renum=False,
                                      return_node_renum=False):
        """
        Compress (if masked) the encapsulated triangulation.

        Returns minimal-length triangles array (*compressed_triangles*) and
        coordinates arrays (*compressed_x*, *compressed_y*) that can still
        describe the unmasked triangles of the encapsulated triangulation.

        Parameters
        ----------
        return_tri_renum : boolean, optional
            Indicates whether a renumbering table to translate the triangle
            numbers from the encapsulated triangulation numbering into the
            new (compressed) renumbering will be returned.
        return_node_renum : boolean, optional
            Indicates whether a renumbering table to translate the nodes
            numbers from the encapsulated triangulation numbering into the
            new (compressed) renumbering will be returned.

        Returns
        -------
        compressed_triangles : array-like
            the returned compressed triangulation triangles
        compressed_x : array-like
            the returned compressed triangulation 1st coordinate
        compressed_y : array-like
            the returned compressed triangulation 2nd coordinate
        tri_renum : array-like of integers
            renumbering table to translate the triangle numbers from the
            encapsulated triangulation into the new (compressed) renumbering.
            -1 for masked triangles (deleted from *compressed_triangles*).
            Returned only if *return_tri_renum* is True.
        node_renum : array-like of integers
            renumbering table to translate the point numbers from the
            encapsulated triangulation into the new (compressed) renumbering.
            -1 for unused points (i.e. those deleted from *compressed_x* and
            *compressed_y*). Returned only if *return_node_renum* is True.

        """
        # Valid triangles and renumbering
        tri_mask = self._triangulation.mask
        compressed_triangles = self._triangulation.get_masked_triangles()
        ntri = self._triangulation.triangles.shape[0]
        tri_renum = self._total_to_compress_renum(tri_mask, ntri)

        # Valid nodes and renumbering
        node_mask = (np.bincount(np.ravel(compressed_triangles),
                                 minlength=self._triangulation.x.size) == 0)
        compressed_x = self._triangulation.x[~node_mask]
        compressed_y = self._triangulation.y[~node_mask]
        node_renum = self._total_to_compress_renum(node_mask)

        # Now renumbering the valid triangles nodes
        compressed_triangles = node_renum[compressed_triangles]

        # 4 cases possible for return
        if not return_tri_renum:
            if not return_node_renum:
                return compressed_triangles, compressed_x, compressed_y
            else:
                return (compressed_triangles, compressed_x, compressed_y,
                        node_renum)
        else:
            if not return_node_renum:
                return (compressed_triangles, compressed_x, compressed_y,
                        tri_renum)
            else:
                return (compressed_triangles, compressed_x, compressed_y,
                        tri_renum, node_renum)

    @staticmethod
    def _total_to_compress_renum(mask, n=None):
        """
        Parameters
        ----------
        mask : 1d boolean array or None
            mask
        n : integer
            length of the mask. Useful only id mask can be None

        Returns
        -------
        renum : integer array
            array so that (`valid_array` being a compressed array
            based on a `masked_array` with mask *mask*) :

                  - For all i such as mask[i] = False:
                    valid_array[renum[i]] = masked_array[i]
                  - For all i such as mask[i] = True:
                    renum[i] = -1 (invalid value)

        """
        if n is None:
            n = np.size(mask)
        if mask is not None:
            renum = np.full(n, -1, dtype=np.int32)  # Default num is -1
            valid = np.arange(n, dtype=np.int32)[~mask]
            renum[valid] = np.arange(np.size(valid, 0), dtype=np.int32)
            return renum
        else:
            return np.arange(n, dtype=np.int32)