"""
Module for povray file format support.

See http://www.povray.org/ for details on the format.
"""
from collections.abc import Mapping, Sequence
from subprocess import check_call, DEVNULL
from os import unlink
from pathlib import Path

import numpy as np

from ase.io.utils import PlottingVariables
from ase.constraints import FixAtoms
from ase import Atoms


def pa(array):
    """Povray array syntax"""
    return '<' + ', '.join(f"{x:>6.2f}" for x in tuple(array)) + '>'


def pc(array):
    """Povray color syntax"""
    if isinstance(array, str):
        return 'color ' + array
    if isinstance(array, float):
        return f'rgb <{array:.2f}>*3'.format(array)
    l = len(array)
    if l > 2 and l < 6:
        return f"rgb{'' if l == 3 else 't' if l == 4 else 'ft'} <" +\
            ', '.join(f"{x:.2f}" for x in tuple(array)) + '>'


def get_bondpairs(atoms, radius=1.1):
    """Get all pairs of bonding atoms

    Return all pairs of atoms which are closer than radius times the
    sum of their respective covalent radii.  The pairs are returned as
    tuples::

      (a, b, (i1, i2, i3))

    so that atoms a bonds to atom b displaced by the vector::

        _     _     _
      i c + i c + i c ,
       1 1   2 2   3 3

    where c1, c2 and c3 are the unit cell vectors and i1, i2, i3 are
    integers."""

    from ase.data import covalent_radii
    from ase.neighborlist import NeighborList
    cutoffs = radius * covalent_radii[atoms.numbers]
    nl = NeighborList(cutoffs=cutoffs, self_interaction=False)
    nl.update(atoms)
    bondpairs = []
    for a in range(len(atoms)):
        indices, offsets = nl.get_neighbors(a)
        bondpairs.extend([(a, a2, offset)
                          for a2, offset in zip(indices, offsets)])
    return bondpairs


def set_high_bondorder_pairs(bondpairs, high_bondorder_pairs=None):
    """Set high bondorder pairs

    Modify bondpairs list (from get_bondpairs((atoms)) to include high
    bondorder pairs.

    Parameters:
    -----------
    bondpairs: List of pairs, generated from get_bondpairs(atoms)
    high_bondorder_pairs: Dictionary of pairs with high bond orders
                          using the following format:
                          { ( a1, b1 ): ( offset1, bond_order1, bond_offset1),
                            ( a2, b2 ): ( offset2, bond_order2, bond_offset2),
                            ...
                          }
                          offset, bond_order, bond_offset are optional.
                          However, if they are provided, the 1st value is
                          offset, 2nd value is bond_order,
                          3rd value is bond_offset """

    if high_bondorder_pairs is None:
        high_bondorder_pairs = dict()
    bondpairs_ = []
    for pair in bondpairs:
        (a, b) = (pair[0], pair[1])
        if (a, b) in high_bondorder_pairs.keys():
            bondpair = [a, b] + [item for item in high_bondorder_pairs[(a, b)]]
            bondpairs_.append(bondpair)
        elif (b, a) in high_bondorder_pairs.keys():
            bondpair = [a, b] + [item for item in high_bondorder_pairs[(b, a)]]
            bondpairs_.append(bondpair)
        else:
            bondpairs_.append(pair)
    return bondpairs_


class POVRAY:
    material_styles_dict = dict(
        simple='finish {phong 0.7}',
        pale=('finish {ambient 0.5 diffuse 0.85 roughness 0.001 '
              'specular 0.200 }'),
        intermediate=('finish {ambient 0.3 diffuse 0.6 specular 0.1 '
                      'roughness 0.04}'),
        vmd=('finish {ambient 0.0 diffuse 0.65 phong 0.1 phong_size 40.0 '
             'specular 0.5 }'),
        jmol=('finish {ambient 0.2 diffuse 0.6 specular 1 roughness 0.001 '
              'metallic}'),
        ase2=('finish {ambient 0.05 brilliance 3 diffuse 0.6 metallic '
              'specular 0.7 roughness 0.04 reflection 0.15}'),
        ase3=('finish {ambient 0.15 brilliance 2 diffuse 0.6 metallic '
              'specular 1.0 roughness 0.001 reflection 0.0}'),
        glass=('finish {ambient 0.05 diffuse 0.3 specular 1.0 '
               'roughness 0.001}'),
        glass2=('finish {ambient 0.01 diffuse 0.3 specular 1.0 '
                'reflection 0.25 roughness 0.001}'),
    )

    def __init__(self, cell, cell_vertices, positions, diameters, colors,
                 image_width, image_height, constraints=tuple(), isosurfaces=[],
                 display=False, pause=True, transparent=True, canvas_width=None,
                 canvas_height=None, camera_dist=50., image_plane=None,
                 camera_type='orthographic', point_lights=[],
                 area_light=[(2., 3., 40.), 'White', .7, .7, 3, 3],
                 background='White', textures=None, transmittances=None,
                 depth_cueing=False, cue_density=5e-3,
                 celllinewidth=0.05, bondlinewidth=0.10, bondatoms=[],
                 exportconstraints=False):
        """
        # x, y is the image plane, z is *out* of the screen
        cell: ase.cell
            cell object
        cell_vertices: 2-d numpy array
            contains the 8 vertices of the cell, each with three coordinates
        positions: 2-d numpy array
            number of atoms length array with three coordinates for positions
        diameters: 1-d numpy array
            diameter of atoms (in order with positions)
        colors: list of str
            colors of atoms (in order with positions)
        image_width: float
            image width in pixels
        image_height: float
            image height in pixels
        constraints: Atoms.constraints
            constraints to be visualized
        isosurfaces: list of POVRAYIsosurface
            composite object to write/render POVRAY isosurfaces
        display: bool
            display while rendering
        pause: bool
            pause when done rendering (only if display)
        transparent: bool
            make background transparent
        canvas_width: int
            width of canvas in pixels
        canvas_height: int
            height of canvas in pixels
        camera_dist: float
            distance from camera to front atom
        image_plane: float
            distance from front atom to image plane
        camera_type: str
            if 'orthographic' perspective, ultra_wide_angle
        point_lights: list of 2-element sequences
            like [[loc1, color1], [loc2, color2],...]
        area_light: 3-element sequence of location (3-tuple), color (str),
                   width (float), height (float),
                   Nlamps_x (int), Nlamps_y (int)
            example [(2., 3., 40.), 'White', .7, .7, 3, 3]
        background: str
            color specification, e.g., 'White'
        textures: list of str
            length of atoms list of texture names
        transmittances: list of floats
            length of atoms list of transmittances of the atoms
        depth_cueing: bool
            whether or not to use depth cueing a.k.a. fog
            use with care - adjust the camera_distance to be closer
        cue_density: float
            if there is depth_cueing, how dense is it (how dense is the fog)
        celllinewidth: float
            radius of the cylinders representing the cell (Ang.)
        bondlinewidth: float
            radius of the cylinders representing bonds (Ang.)
        bondatoms: list of lists (polymorphic)
            [[atom1, atom2], ... ] pairs of bonding atoms
             For bond order > 1 = [[atom1, atom2, offset,
                                    bond_order, bond_offset],
                                   ... ]
             bond_order: 1, 2, 3 for single, double,
                          and triple bond
             bond_offset: vector for shifting bonds from
                           original position. Coordinates are
                           in Angstrom unit.
        exportconstraints: bool
            honour FixAtoms and mark?"""

        # attributes from initialization
        self.area_light = area_light
        self.background = background
        self.bondatoms = bondatoms
        self.bondlinewidth = bondlinewidth
        self.camera_dist = camera_dist
        self.camera_type = camera_type
        self.celllinewidth = celllinewidth
        self.cue_density = cue_density
        self.depth_cueing = depth_cueing
        self.display = display
        self.exportconstraints = exportconstraints
        self.isosurfaces = isosurfaces
        self.pause = pause
        self.point_lights = point_lights
        self.textures = textures
        self.transmittances = transmittances
        self.transparent = transparent

        self.image_width = image_width
        self.image_height = image_height
        self.colors = colors
        self.cell = cell
        self.diameters = diameters

        # calculations based on passed inputs

        z0 = positions[:, 2].max()
        self.offset = (image_width / 2, image_height / 2, z0)
        self.positions = positions - self.offset

        if cell_vertices is not None:
            self.cell_vertices = cell_vertices - self.offset
            self.cell_vertices.shape = (2, 2, 2, 3)
        else:
            self.cell_vertices = None

        ratio = float(self.image_width) / self.image_height
        if canvas_width is None:
            if canvas_height is None:
                self.canvas_width = min(self.image_width * 15, 640)
                self.canvas_height = min(self.image_height * 15, 640)
            else:
                self.canvas_width = canvas_height * ratio
                self.canvas_height = canvas_height
        elif canvas_height is None:
            self.canvas_width = canvas_width
            self.canvas_height = self.canvas_width / ratio
        else:
            raise RuntimeError("Can't set *both* width and height!")

        # Distance to image plane from camera
        if image_plane is None:
            if self.camera_type == 'orthographic':
                self.image_plane = 1 - self.camera_dist
            else:
                self.image_plane = 0
        self.image_plane += self.camera_dist

        self.constrainatoms = []
        for c in constraints:
            if isinstance(c, FixAtoms):
                # self.constrainatoms.extend(c.index) # is this list-like?
                for n, i in enumerate(c.index):
                    self.constrainatoms += [i]

    @classmethod
    def from_PlottingVariables(cls, pvars, **kwargs):
        cell = pvars.cell
        cell_vertices = pvars.cell_vertices
        if 'colors' in kwargs.keys():
            colors = kwargs.pop('colors')
        else:
            colors = pvars.colors
        diameters = pvars.d
        image_height = pvars.h
        image_width = pvars.w
        positions = pvars.positions
        constraints = pvars.constraints
        return cls(cell=cell, cell_vertices=cell_vertices, colors=colors,
                   constraints=constraints, diameters=diameters,
                   image_height=image_height, image_width=image_width,
                   positions=positions, **kwargs)

    @classmethod
    def from_atoms(cls, atoms, **kwargs):
        return cls.from_plotting_variables(
            PlottingVariables(atoms, scale=1.0), **kwargs)

    def write_ini(self, path):
        """Write ini file."""

        ini_str = f"""\
Input_File_Name={path.with_suffix('.pov').name}
Output_to_File=True
Output_File_Type=N
Output_Alpha={'on' if self.transparent else 'off'}
; if you adjust Height, and width, you must preserve the ratio
; Width / Height = {self.canvas_width/self.canvas_height:f}
Width={self.canvas_width}
Height={self.canvas_height}
Antialias=True
Antialias_Threshold=0.1
Display={self.display}
Pause_When_Done={self.pause}
Verbose=False
"""
        with open(path, 'w') as fd:
            fd.write(ini_str)
        return path

    def write_pov(self, path):
        """Write pov file."""

        point_lights = '\n'.join(f"light_source {{{pa(loc)} {pc(rgb)}}}"
                                 for loc, rgb in self.point_lights)

        area_light = ''
        if self.area_light is not None:
            loc, color, width, height, nx, ny = self.area_light
            area_light += f"""\nlight_source {{{pa(loc)} {pc(color)}
  area_light <{width:.2f}, 0, 0>, <0, {height:.2f}, 0>, {nx:n}, {ny:n}
  adaptive 1 jitter}}"""

        fog = ''
        if self.depth_cueing and (self.cue_density >= 1e-4):
            # same way vmd does it
            if self.cue_density > 1e4:
                # larger does not make any sense
                dist = 1e-4
            else:
                dist = 1. / self.cue_density
            fog += f'fog {{fog_type 1 distance {dist:.4f} '\
                   f'color {pc(self.background)}}}'

        mat_style_keys = (f'#declare {k} = {v}'
                          for k, v in self.material_styles_dict.items())
        mat_style_keys = '\n'.join(mat_style_keys)

        # Draw unit cell
        cell_vertices = ''
        if self.cell_vertices is not None:
            for c in range(3):
                for j in ([0, 0], [1, 0], [1, 1], [0, 1]):
                    p1 = self.cell_vertices[tuple(j[:c]) + (0,) + tuple(j[c:])]
                    p2 = self.cell_vertices[tuple(j[:c]) + (1,) + tuple(j[c:])]

                    distance = np.linalg.norm(p2 - p1)
                    if distance < 1e-12:
                        continue

                    cell_vertices += f'cylinder {{{pa(p1)}, {pa(p2)}, '\
                                     f'Rcell pigment {{Black}}}}\n'
                    # all strings are f-strings for consistency
            cell_vertices = cell_vertices.strip('\n')

        # Draw atoms
        a = 0
        atoms = ''
        for loc, dia, col in zip(self.positions, self.diameters, self.colors):
            tex = 'ase3'
            trans = 0.
            if self.textures is not None:
                tex = self.textures[a]
            if self.transmittances is not None:
                trans = self.transmittances[a]
            atoms += f'atom({pa(loc)}, {dia/2.:.2f}, {pc(col)}, '\
                     f'{trans}, {tex}) // #{a:n}\n'
            a += 1
        atoms = atoms.strip('\n')

        # Draw atom bonds
        bondatoms = ''
        for pair in self.bondatoms:
            # Make sure that each pair has 4 componets: a, b, offset,
            #                                           bond_order, bond_offset
            # a, b: atom index to draw bond
            # offset: original meaning to make offset for mid-point.
            # bond_oder: if not supplied, set it to 1 (single bond).
            #            It can be  1, 2, 3, corresponding to single,
            #            double, triple bond
            # bond_offset: displacement from original bond position.
            #              Default is (bondlinewidth, bondlinewidth, 0)
            #              for bond_order > 1.
            if len(pair) == 2:
                a, b = pair
                offset = (0, 0, 0)
                bond_order = 1
                bond_offset = (0, 0, 0)
            elif len(pair) == 3:
                a, b, offset = pair
                bond_order = 1
                bond_offset = (0, 0, 0)
            elif len(pair) == 4:
                a, b, offset, bond_order = pair
                bond_offset = (self.bondlinewidth, self.bondlinewidth, 0)
            elif len(pair) > 4:
                a, b, offset, bond_order, bond_offset = pair
            else:
                raise RuntimeError('Each list in bondatom must have at least '
                                   '2 entries. Error at %s' % pair)

            if len(offset) != 3:
                raise ValueError('offset must have 3 elements. '
                                 'Error at %s' % pair)
            if len(bond_offset) != 3:
                raise ValueError('bond_offset must have 3 elements. '
                                 'Error at %s' % pair)
            if bond_order not in [0, 1, 2, 3]:
                raise ValueError('bond_order must be either 0, 1, 2, or 3. '
                                 'Error at %s' % pair)

            # Up to here, we should have all a, b, offset, bond_order,
            # bond_offset for all bonds.

            # Rotate bond_offset so that its direction is 90 deg. off the bond
            # Utilize Atoms object to rotate
            if bond_order > 1 and np.linalg.norm(bond_offset) > 1.e-9:
                tmp_atoms = Atoms('H3')
                tmp_atoms.set_cell(self.cell)
                tmp_atoms.set_positions([
                    self.positions[a],
                    self.positions[b],
                    self.positions[b] + np.array(bond_offset),
                ])
                tmp_atoms.center()
                tmp_atoms.set_angle(0, 1, 2, 90)
                bond_offset = tmp_atoms[2].position - tmp_atoms[1].position

            R = np.dot(offset, self.cell)
            mida = 0.5 * (self.positions[a] + self.positions[b] + R)
            midb = 0.5 * (self.positions[a] + self.positions[b] - R)
            if self.textures is not None:
                texa = self.textures[a]
                texb = self.textures[b]
            else:
                texa = texb = 'ase3'

            if self.transmittances is not None:
                transa = self.transmittances[a]
                transb = self.transmittances[b]
            else:
                transa = transb = 0.

            # draw bond, according to its bond_order.
            # bond_order == 0: No bond is plotted
            # bond_order == 1: use original code
            # bond_order == 2: draw two bonds, one is shifted by bond_offset/2,
            #                  and another is shifted by -bond_offset/2.
            # bond_order == 3: draw two bonds, one is shifted by bond_offset,
            #                  and one is shifted by -bond_offset, and the
            #                  other has no shift.
            # To shift the bond, add the shift to the first two coordinate in
            # write statement.

            posa = self.positions[a]
            posb = self.positions[b]
            cola = self.colors[a]
            colb = self.colors[b]

            if bond_order == 1:
                draw_tuples = (posa, mida, cola, transa, texa),\
                              (posb, midb, colb, transb, texb)

            elif bond_order == 2:
                bs = [x / 2 for x in bond_offset]
                draw_tuples = (posa - bs, mida - bs, cola, transa, texa),\
                              (posb - bs, midb - bs, colb, transb, texb),\
                              (posa + bs, mida + bs, cola, transa, texa),\
                              (posb + bs, midb + bs, colb, transb, texb)

            elif bond_order == 3:
                bs = bond_offset
                draw_tuples = (posa, mida, cola, transa, texa),\
                              (posb, midb, colb, transb, texb),\
                              (posa + bs, mida + bs, cola, transa, texa),\
                              (posb + bs, midb + bs, colb, transb, texb),\
                              (posa - bs, mida - bs, cola, transa, texa),\
                              (posb - bs, midb - bs, colb, transb, texb)

            bondatoms += ''.join(f'cylinder {{{pa(p)}, '
                                 f'{pa(m)}, Rbond texture{{pigment '
                                 f'{{color {pc(c)} '
                                 f'transmit {tr}}} finish{{{tx}}}}}}}\n'
                                 for p, m, c, tr, tx in
                                 draw_tuples)

        bondatoms = bondatoms.strip('\n')

        # Draw constraints if requested
        constraints = ''
        if self.exportconstraints:
            for a in self.constrainatoms:
                dia = self.diameters[a]
                loc = self.positions[a]
                trans = 0.0
                if self.transmittances is not None:
                    trans = self.transmittances[a]
                constraints += f'constrain({pa(loc)}, {dia/2.:.2f}, Black, '\
                    f'{trans}, {tex}) // #{a:n} \n'
        constraints = constraints.strip('\n')

        pov = f"""#include "colors.inc"
#include "finish.inc"

global_settings {{assumed_gamma 1 max_trace_level 6}}
background {{{pc(self.background)}{' transmit 1.0' if self.transparent else ''}}}
camera {{{self.camera_type}
  right -{self.image_width:.2f}*x up {self.image_height:.2f}*y
  direction {self.image_plane:.2f}*z
  location <0,0,{self.camera_dist:.2f}> look_at <0,0,0>}}
{point_lights}
{area_light if area_light != '' else '// no area light'}
{fog if fog != '' else '// no fog'}
{mat_style_keys}
#declare Rcell = {self.celllinewidth:.3f};
#declare Rbond = {self.bondlinewidth:.3f};

#macro atom(LOC, R, COL, TRANS, FIN)
  sphere{{LOC, R texture{{pigment{{color COL transmit TRANS}} finish{{FIN}}}}}}
#end
#macro constrain(LOC, R, COL, TRANS FIN)
union{{torus{{R, Rcell rotate 45*z texture{{pigment{{color COL transmit TRANS}} finish{{FIN}}}}}}
     torus{{R, Rcell rotate -45*z texture{{pigment{{color COL transmit TRANS}} finish{{FIN}}}}}}
     translate LOC}}
#end

{cell_vertices if cell_vertices != '' else '// no cell vertices'}
{atoms}
{bondatoms}
{constraints if constraints != '' else '// no constraints'}
"""  # noqa: E501

        with open(path, 'w') as fd:
            fd.write(pov)

        return path

    def write(self, pov_path):
        pov_path = require_pov(pov_path)
        ini_path = pov_path.with_suffix('.ini')
        self.write_ini(ini_path)
        self.write_pov(pov_path)
        if self.isosurfaces is not None:
            with open(pov_path, 'a') as fd:
                for iso in self.isosurfaces:
                    fd.write(iso.format_mesh())
        return POVRAYInputs(ini_path)


def require_pov(path):
    path = Path(path)
    if path.suffix != '.pov':
        raise ValueError(f'Expected .pov path, got {path}')
    return path


class POVRAYInputs:
    def __init__(self, path):
        self.path = path

    def render(self, povray_executable='povray', stderr=DEVNULL,
               clean_up=False):
        cmd = [povray_executable, str(self.path)]

        check_call(cmd, stderr=stderr)
        png_path = self.path.with_suffix('.png').absolute()
        if not png_path.is_file():
            raise RuntimeError(f'Povray left no output PNG file "{png_path}"')

        if clean_up:
            unlink(self.path)
            unlink(self.path.with_suffix('.pov'))

        return png_path


class POVRAYIsosurface:
    def __init__(self, density_grid, cut_off, cell, cell_origin,
                 closed_edges=False, gradient_ascending=False,
                 color=(0.85, 0.80, 0.25, 0.2), material='ase3'):
        """
        density_grid: 3D float ndarray
            A regular grid on that spans the cell. The first dimension
            corresponds to the first cell vector and so on.
        cut_off: float
            The density value of the isosurface.
        cell: 2D float ndarray or ASE cell object
            The 3 vectors which give the cell's repetition
        cell_origin: 4 float tuple
            The cell origin as used by POVRAY object
        closed_edges: bool
            Setting this will fill in isosurface edges at the cell boundaries.
            Filling in the edges can help with visualizing
            highly porous structures.
        gradient_ascending: bool
            Lets you pick the area you want to enclose, i.e., should the denser
            or less dense area be filled in.
        color: povray color string, float, or float tuple
            1 float is interpreted as grey scale, a 3 float tuple is rgb,
            4 float tuple is rgbt, and 5 float tuple is rgbft, where
            t is transmission fraction and f is filter fraction.
            Named Povray colors are set in colors.inc
            (http://wiki.povray.org/content/Reference:Colors.inc)
        material: string
            Can be a finish macro defined by POVRAY.material_styles
            or a full Povray material {...} specification. Using a
            full material specification willoverride the color parameter.
        """

        self.gradient_direction = 'ascent' if gradient_ascending else 'descent'
        self.color = color
        self.material = material
        self.closed_edges = closed_edges
        self._cut_off = cut_off

        if self.gradient_direction == 'ascent':
            cv = 2 * cut_off
        else:
            cv = 0

        if closed_edges:
            shape_old = density_grid.shape
            # since well be padding, we need to keep the data at origin
            cell_origin += -(1.0 / np.array(shape_old)) @ cell
            density_grid = np.pad(
                density_grid, pad_width=(
                    1,), mode='constant', constant_values=cv)
            shape_new = density_grid.shape
            s = np.array(shape_new) / np.array(shape_old)
            cell = cell @ np.diag(s)

        self.cell = cell
        self.cell_origin = cell_origin
        self.density_grid = density_grid
        self.spacing = tuple(1.0 / np.array(self.density_grid.shape))

        scaled_verts, faces, normals, values = self.compute_mesh(
            self.density_grid,
            self.cut_off,
            self.spacing,
            self.gradient_direction)

        # The verts are scaled by default, this is the super easy way of
        # distributing them in real space but it's easier to do affine
        # transformations/rotations on a unit cube so I leave it like that
        # verts = scaled_verts.dot(self.cell)
        self.verts = scaled_verts
        self.faces = faces

    @property
    def cut_off(self):
        return self._cut_off

    @cut_off.setter
    def cut_off(self, value):
        raise Exception("Use the set_cut_off method")

    def set_cut_off(self, value):
        self._cut_off = value

        if self.gradient_direction == 'ascent':
            cv = 2 * self.cut_off
        else:
            cv = 0

        if self.closed_edges:
            shape_old = self.density_grid.shape
            # since well be padding, we need to keep the data at origin
            self.cell_origin += -(1.0 / np.array(shape_old)) @ self.cell
            self.density_grid = np.pad(
                self.density_grid, pad_width=(
                    1,), mode='constant', constant_values=cv)
            shape_new = self.density_grid.shape
            s = np.array(shape_new) / np.array(shape_old)
            self.cell = self.cell @ np.diag(s)

        self.spacing = tuple(1.0 / np.array(self.density_grid.shape))

        scaled_verts, faces, _, _ = self.compute_mesh(
            self.density_grid,
            self.cut_off,
            self.spacing,
            self.gradient_direction)

        self.verts = scaled_verts
        self.faces = faces

    @classmethod
    def from_POVRAY(cls, povray, density_grid, cut_off, **kwargs):
        return cls(cell=povray.cell,
                   cell_origin=povray.cell_vertices[0, 0, 0],
                   density_grid=density_grid,
                   cut_off=cut_off, **kwargs)

    @staticmethod
    def wrapped_triples_section(triple_list,
                                triple_format="<{:f}, {:f}, {:f}>".format,
                                triples_per_line=4):

        triples = [triple_format(*x) for x in triple_list]
        n = len(triples)
        s = ''
        tpl = triples_per_line
        c = 0

        while c < n - tpl:
            c += tpl
            s += '\n     '
            s += ', '.join(triples[c - tpl:c])
        s += '\n    '
        s += ', '.join(triples[c:])
        return s

    @staticmethod
    def compute_mesh(density_grid, cut_off, spacing, gradient_direction):
        """

        Import statement is in this method and not file header
        since few users will use isosurface rendering.

        Returns scaled_verts, faces, normals, values. See skimage docs.

        """

        from skimage import measure
        return measure.marching_cubes_lewiner(
            density_grid,
            level=cut_off,
            spacing=spacing,
            gradient_direction=gradient_direction,
            allow_degenerate=False)

    def format_mesh(self):
        """Returns a formatted data output for POVRAY files

        Example:
        material = '''
          material { // This material looks like pink jelly
            texture {
              pigment { rgbt <0.8, 0.25, 0.25, 0.5> }
              finish{ diffuse 0.85 ambient 0.99 brilliance 3 specular 0.5 roughness 0.001
                reflection { 0.05, 0.98 fresnel on exponent 1.5 }
                conserve_energy
              }
            }
            interior { ior 1.3 }
          }
          photons {
              target
              refraction on
              reflection on
              collect on
          }'''
        """  # noqa: E501

        if self.material in POVRAY.material_styles_dict:
            material = f"""material {{
        texture {{
          pigment {{ {pc(self.color)} }}
          finish {{ {self.material} }}
        }}
      }}"""
        else:
            material = self.material

        # Start writing the mesh2
        vertex_vectors = self.wrapped_triples_section(
            triple_list=self.verts,
            triple_format="<{:f}, {:f}, {:f}>".format,
            triples_per_line=4)

        face_indices = self.wrapped_triples_section(
            triple_list=self.faces,
            triple_format="<{:n}, {:n}, {:n}>".format,
            triples_per_line=5)

        cell = self.cell
        cell_or = self.cell_origin
        mesh2 = f"""\n\nmesh2 {{
    vertex_vectors {{  {len(self.verts):n},
    {vertex_vectors}
    }}
    face_indices {{ {len(self.faces):n},
    {face_indices}
    }}
{material if material != '' else '// no material'}
  matrix < {cell[0][0]:f}, {cell[0][1]:f}, {cell[0][2]:f},
           {cell[1][0]:f}, {cell[1][1]:f}, {cell[1][2]:f},
           {cell[2][0]:f}, {cell[2][1]:f}, {cell[2][2]:f},
           {cell_or[0]:f}, {cell_or[1]:f}, {cell_or[2]:f}>
    }}
    """
        return mesh2


def pop_deprecated(dct, name):
    import warnings
    if name in dct:
        del dct[name]
        warnings.warn(f'The "{name}" keyword of write_pov() is deprecated '
                      'and has no effect; this will raise an error in the '
                      'future.', FutureWarning)


def write_pov(filename, atoms, *,
              povray_settings=None, isosurface_data=None,
              **generic_projection_settings):

    for name in ['run_povray', 'povray_path', 'stderr', 'extras']:
        pop_deprecated(generic_projection_settings, name)

    if povray_settings is None:
        povray_settings = {}

    pvars = PlottingVariables(atoms, scale=1.0, **generic_projection_settings)
    pov_obj = POVRAY.from_PlottingVariables(pvars, **povray_settings)

    if isinstance(isosurface_data, Mapping):
        pov_obj.isosurfaces = [POVRAYIsosurface.from_POVRAY(
            pov_obj, **isosurface_data)]
    elif isinstance(isosurface_data, Sequence):
        pov_obj.isosurfaces = [POVRAYIsosurface.from_POVRAY(
            pov_obj, **isodata) for isodata in isosurface_data]

    return pov_obj.write(filename)