#
##
##  This file is part of pyFormex 1.0.5  (Sat Feb 16 10:40:32 CET 2019)
##  pyFormex is a tool for generating, manipulating and transforming 3D
##  geometrical models by sequences of mathematical operations.
##  Home page: http://pyformex.org
##  Project page:  http://savannah.nongnu.org/projects/pyformex/
##  Copyright 2004-2018 (C) Benedict Verhegghe (benedict.verhegghe@ugent.be)
##  Distributed under the GNU General Public License version 3 or later.
##
##  This program is free software: you can redistribute it and/or modify
##  it under the terms of the GNU General Public License as published by
##  the Free Software Foundation, either version 3 of the License, or
##  (at your option) any later version.
##
##  This program is distributed in the hope that it will be useful,
##  but WITHOUT ANY WARRANTY; without even the implied warranty of
##  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
##  GNU General Public License for more details.
##
##  You should have received a copy of the GNU General Public License
##  along with this program.  If not, see http://www.gnu.org/licenses/.
##
"""Operations on triangulated surfaces using GTS functions.

This module provides access to GTS from insisde pyFormex.
"""
from __future__ import absolute_import, division, print_function


import pyformex as pf
from pyformex.coords import *
from pyformex import utils
import os
from pyformex.multi import multitask

# No warning when in sphinx/docmodule mode
if not pf.sphinx:
    if not utils.hasExternal('gts-bin'):
        utils.warn("error_no_gts_bin")
    if not utils.hasExternal('gts-extra'):
        utils.warn("error_no_gts_extra")


#
# gts commands used:
#   in Debian package: stl2gts gts2stl gtscheck
#   not in Debian package: gtssplit gtscoarsen gtsrefine gtssmooth gtsinside
#


def read_gts_intersectioncurve(fn):
    import re
    from pyformex.formex import Formex
    RE = re.compile("^VECT 1 2 0 2 0 (?P<data>.*)$")
    r = []
    for line in open(fn, 'r'):
        m = RE.match(line)
        if m:
            r.append(m.group('data'))
    nelems = len(r)
    x = fromstring('\n'.join(r), sep=' ').reshape(-1, 2, 3)
    F = Formex(x)
    return F


def boolean(self,surf,op,check=False,verbose=False):
    """Perform a boolean operation with another surface.

    Boolean operations between surfaces are a basic operation in
    free surface modeling. Both surfaces should be closed orientable
    non-intersecting manifolds.
    Use the :meth:`check` method to find out.

    The boolean operations are set operations on the enclosed volumes:
    union('+'), difference('-') or intersection('*').

    Parameters:

    - `surf`: a closed manifold surface
    - `op`: boolean operation: one of '+', '-' or '*'.
    - `filt`: a filter command to be executed on the gtsset output
    - `ext`: extension of the result file
    - `check`: boolean: check that the surfaces are not self-intersecting;
      if one of them is, the set of self-intersecting faces is written
      (as a GtsSurface) on standard output
    - `verbose`: boolean: print statistics about the surface

    Returns: a closed manifold TriSurface

    .. note: This uses the external command 'gtsset'
    """
    return self.gtsset(surf, op, filt = '', ext='.gts', check=check, verbose=verbose)


def intersection(self,surf,check=False,verbose=False):
    """Return the intersection curve of two surfaces.

    Boolean operations between surfaces are a basic operation in
    free surface modeling. Both surfaces should be closed orientable
    non-intersecting manifolds.
    Use the :meth:`check` method to find out.

    Parameters:

    - `surf`: a closed manifold surface
    - `check`: boolean: check that the surfaces are not self-intersecting;
      if one of them is, the set of self-intersecting faces is written
      (as a GtsSurface) on standard output
    - `verbose`: boolean: print statistics about the surface

    Returns: a list of intersection curves.
    """
    return self.gtsset(surf, op='*', ext='.list', curve=True, check=check, verbose=verbose)


def gtsset(self,surf,op,filt='',ext='.tmp',curve=False,check=False,verbose=False):
    """_Perform the boolean/intersection methods.

    See the boolean/intersection methods for more info.
    Parameters not explained there:

    - curve: if True, an intersection curve is computed, else the surface.

    Returns the resulting TriSurface (curve=False), or a plex-2 Formex
    (curve=True), or None if the input surfaces do not intersect.
    """
    # import here to avoid circular import
    from pyformex.trisurface import TriSurface
    op = {'+':'union', '-':'diff', '*':'inter'}[op]
    options = ''
    if curve:
        options += '-i'
    if check:
        options += ' -s'
    if not verbose:
        options += ' -v'
    tmp = utils.tempFile(suffix='.gts').name
    tmp1 = utils.tempFile(suffix='.gts').name
    tmp2 = utils.tempFile(suffix=ext).name
    print("Writing temp file %s" % tmp)
    self.write(tmp, 'gts')
    print("Writing temp file %s" % tmp1)
    surf.write(tmp1, 'gts')
    print("Performing boolean operation")
    cmd = "gtsset %s %s %s %s %s" % (options, op, tmp, tmp1, filt)
    P = utils.system(cmd, stdout=open(tmp2, 'w'))
    os.remove(tmp)
    os.remove(tmp1)
    if P.sta or verbose:
        print(P.out)
    if P.sta:
        print(P.err)
        return None
    print("Reading result from %s" % tmp2)
    if curve:
        res = read_gts_intersectioncurve(tmp2)
    else:
        res = TriSurface.read(tmp2)
    os.remove(tmp2)
    return res

# TODO: change pts to a Coords
def gtsinside(self,pts,dir=0):
    """_Test whether points are inside the surface.

    pts is a plex-1 Formex.
    dir is the shooting direction.

    Returns a list of point numbers that are inside.
    This is not intended to be used directly. Use inside instead
    """
    import os
    S = self.rollAxes(dir)
    P = pts.rollAxes(dir)
    tmp = utils.tempFile(suffix='.gts').name
    tmp1 = utils.tempFile(suffix='.dta').name
    tmp2 = utils.tempFile(suffix='.out').name
    #print("Writing temp file %s" % tmp)
    S.write(tmp, 'gts')
    #print("Writing temp file %s" % tmp1)
    with open(tmp1, 'w') as f:
        P.coords.tofile(f, sep=' ')
        f.write('\n')

    #print("Performing inside testing")
    cmd = "gtsinside %s %s" % (tmp, tmp1)
    P = utils.command(cmd, stdout=open(tmp2, 'w'))
    os.remove(tmp)
    os.remove(tmp1)
    if P.sta:
        #print("An error occurred during the testing.\nSee file %s for more details." % tmp2)
        print(P.out)
        return None
    #print("Reading results from %s" % tmp2)
    ind = fromfile(tmp2, sep=' ', dtype=Int)
    os.remove(tmp2)
    return ind


def inside(self,pts,atol='auto',multi=True):
    """Test which of the points pts are inside the surface.

    Parameters:

    - `pts`: a (usually 1-plex) Formex or a data structure that can be used
      to initialize a Formex.

    Returns an integer array with the indices of the points that are
    inside the surface. The indices refer to the onedimensional list
    of points as obtained from pts.points().
    """
    from pyformex.formex import Formex
    if not isinstance(pts, Formex):
        pts = Formex(pts)

    if atol == 'auto':
        atol = pts.dsize()*0.001

    # determine bbox of common space of surface and points
    bb = bboxIntersection(self, pts)
    if (bb[0] > bb[1]).any():
        # No bbox intersection: no points inside
        return array([], dtype=Int)

    # Limit the points to the common part
    # Add point numbers as property, to allow return of original numbers
    pts.setProp(arange(pts.nelems()))
    pts = pts.clip(pts.testBbox(bb, atol=atol))

    ins = zeros((3, pts.nelems()), dtype=bool)
    if pf.scriptMode == 'script' or not multi:
        for i in range(3):
            dirs = roll(arange(3), -i)[1:]
            # clip the surface perpendicular to the shooting direction
            # !! gtsinside seems to fail sometimes when using clipping
            S = self#.clip(testBbox(self,bb,dirs=dirs,atol=atol),compact=False)
            # find inside points shooting in direction i
            ok = gtsinside(S, pts, dir=i)
            ins[i, ok] = True
    else:
        tasks = [(gtsinside, (self, pts, i)) for i in range(3)]
        ind = multitask(tasks, 3)
        for i in range(3):
            ins[i, ind[i]] = True

    ok = where(ins.sum(axis=0) > 1)[0]
    return pts.prop[ok]


# End
