#
##
##  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/.
##
"""A structured collection of 3D coordinates.

The :mod:`coords` module defines the :class:`Coords` class, which is the basic
data structure in pyFormex to store the coordinates of points in a 3D space.

This module implements a data class for storing large sets of 3D coordinates
and provides an extensive set of methods for transforming these coordinates.
Most of pyFormex's classes which represent geometry (e.g.
:class:`~formex.Formex`, :class:`~mesh.Mesh`, :class:`~trisurface.TriSurface`,
:class:`~plugins.curve.Curve`) use a
:class:`Coords` object to store their coordinates, and thus inherit all the
transformation methods of this class.

While the user will mostly use the higher level classes, he might occasionally
find good reason to use the :class:`Coords` class directly as well.
"""
from __future__ import absolute_import, division, print_function

import numpy as np

from pyformex import utils

from pyformex.arraytools import *


###########################################################################
##
##   class Coords
##
#########################
#

class Coords(ndarray):
    #
    # :DEV
    # Because we have a __new__ constructor and no __init__,
    # we have to put the signature of the object creation explicitely
    # in the first line of the docstring.
    #
    """Coords(data=None,dtyp=Float,copy=False)

    A structured collection of points in a 3D cartesian space.

    The :class:`Coords` class is the basic data structure used throughout
    pyFormex to store the coordinates of points in a 3D space.
    It is used by other classes, such as :class:`~formex.Formex`,
    :class:`~mesh.Mesh`, :class:`~trisurface.TriSurface`,
    :class:`~plugins.curve.Curve`, which thus inherit the same transformation
    capabilities. Applications will mostly use the higher level
    classes, which have more elaborated consistency checking
    and error handling.

    :class:`Coords` is implemented as a subclass of :class:`numpy.ndarray`,
    and thus inherits all its methods and atttributes.
    The last axis of the :class:`Coords` however always has a length equal
    to 3. Each set of 3 values along the last axis are the coordinates (in a
    global 3D cartesian coordinate system) of a single point in space.
    The full Coords array thus is a collection of points. It the array is
    2-dimensional, the Coords is a flat list of points. But if the array
    has more dimensions, the collection of points itself becomes structured.

    The float datatype is only checked at creation time. It is the
    responsibility of the user to keep this consistent
    throughout the lifetime of the object.

    Note
    ----
    Methods that transform a Coords object, like
    :meth:`scale`, :meth:`translate`, :meth:`rotate`, ...
    do not change the original Coords object, but return a new object.
    Some methods however have an `inplace` option that allows the user
    to force coordinates to be changed in place. This option is seldom
    used however: rather we conveniently use statements like::

        X = X.some_transform()

    and Python can immediately free and recollect the memory used for the old
    object X.


    Parameters
    ----------
    data: float array-like, or string
        Data to initialize the Coords. The last axis should have a length
        of 1, 2 or 3, but will be expanded to 3 if it is less, filling the
        missing coordinates with zeros. Thus, if you only specify two
        coordinates, all points are lying in the z=0 plane. Specifying only
        one coordinate creates points along the x-axis.

        As a convenience, data may also be entered as a string, which will be
        passed to the :func:`pattern` function to create the actual
        coordinates of the points.

        If no data are provided, an empty Coords with shape (0,3) is created.
    dtyp: float datatype, optional
        It not provided, the datatype of ``data`` is used, or the default
        :py:attr:`~arraytools.Float` (which is equivalent to
        :data:`numpy.float32`).
    copy: bool
        If True, the data are copied. The default setting will try to use
        the original data if possible, e.g. if `data` is a correctly shaped and
        typed :class:`numpy.ndarray`.

    Returns
    -------
    : Coords
        An instance of the Coords class, which is basically an ndarray
        of floats, with the last axis having a length of 3.


    The Coords instance has a number of attributes that provide views on
    (part of) the data. They are a notational convenience over using indexing.
    These attributes can be used to set all or some of the coordinates
    by direct assignment. The assigned data should however be broadcast
    compatible with the assigned shape: the shape of the Coords can
    not be changed.

    Attributes
    ----------
    xyz: float array
        The full coordinate array as an ndarray.
    x: float array
        The X coordinates of the points as an ndarray with shape
        :meth:`pshape()`.
    y: float array
        The Y coordinates of the points as an ndarray with shape
        :meth:`pshape()`.
    z: float array
        The Z coordinates of the points as an ndarray with shape
        :meth:`pshape()`.
    xy: float array
        The X and Y coordinates of the points as an ndarray with shape
        :meth:`pshape()` + (2,).
    xz: float array
        The X and Z coordinates of the points as an ndarray with shape
        :meth:`pshape()` + (2,).
    yz: float array
        The Y and Z coordinates of the points as an ndarray with shape
        :meth:`pshape()` + (2,).


    Examples
    --------
    >>> Coords([1.,2.])
    Coords([ 1., 2., 0.])
    >>> X = Coords(arange(6).reshape(2,3))
    >>> print(X)
    [[ 0. 1. 2.]
     [ 3. 4. 5.]]
    >>> print(X.y)
    [ 1. 4.]
    >>> X.z[1] = 9.
    >>> print(X)
    [[ 0. 1. 2.]
     [ 3. 4. 9.]]
    >>> print(X.xz)
    [[ 0. 2.]
     [ 3. 9.]]
    >>> X.x = 0.
    >>> print(X)
    [[ 0. 1. 2.]
     [ 0. 4. 9.]]

    >>> Y = Coords(X)             # Y shares its data with X
    >>> Z = Coords(X, copy=True)  # Z is independent
    >>> Y.y = 5
    >>> Z.z = 6
    >>> print(X)
    [[ 0.  5.  2.]
     [ 0.  5.  9.]]
    >>> print(Y)
    [[ 0.  5.  2.]
     [ 0.  5.  9.]]
    >>> print(Z)
    [[ 0.  1.  6.]
     [ 0.  4.  6.]]
    >>> X.coords is X
    True
    >>> Z.xyz = [1,2,3]
    >>> print(Z)
    [[ 1.  2.  3.]
     [ 1.  2.  3.]]

    >>> print(Coords('0123'))  # initialize with string
    [[ 0. 0. 0.]
     [ 1. 0. 0.]
     [ 1. 1. 0.]
     [ 0. 1. 0.]]

    """

    _exclude_members_ = [ 'actor' ]

    #
    # TODO
    # Because the Coords class is sticky, results that are not conforming
    # to the requirements of being a Coords array, should be converted to
    # the general array class:   e.g.   return asarray(T)
    #
    # This could probably NOT be done in an __array_finalize__ method
    # Maybe by oveloading __getitem__, but should we bother???
    #

    def __new__(clas, data=None, dtyp=Float, copy=False):
        """Create a new instance of :class:`Coords`."""
        if data is None:
            # create an empty array : we need at least a 2D array
            # because we want the last axis to have length 3 and
            # we also need an axis with length 0 to have size 0
            ar = ndarray((0, 3), dtype=dtyp)
        else:
            # turn the data into an array, and copy if requested
            # DO NOT ADD ndmin=1 HERE ! (see below)
            if isinstance(data,str):
                data = pattern(data,aslist=True)
            ar = array(data, dtype=dtyp, copy=copy)

        #
        # The Coords object needs to be at least 1-D array, no a scalar
        # We could force 'ar' above to be at least 1-D, but that would
        # turn every scalar into a 1-D vector, which would circumvent
        # detection of input errors (e.g. with translation, where input
        # can be either a vector or an axis number)
        #
        if ar.ndim == 0:
            raise ValueError("Expected array data, not a scalar")

        if ar.shape[-1] == 3:
            pass
        elif ar.shape[-1] in [1, 2]:
            # make last axis length 3, adding 0 values
            ar = growAxis(ar, 3-ar.shape[-1], -1)
        elif ar.shape[-1] == 0:
            # allow empty coords objects
            ar = ar.reshape(0, 3)
        else:
            raise ValueError("Expected a length 1,2 or 3 for last array axis")

        # Make sure dtype is a float type
        if ar.dtype.kind != 'f':
            ar = ar.astype(Float)

        # Transform 'subarr' from an ndarray to our new subclass.
        ar = ar.view(clas)

        return ar


    # for now, very simple detector of not-full coordinates
    def __getitem__(self,i):
        res = np.ndarray.__getitem__(self,i)
        if res.ndim < 1 or res.shape[-1] != 3:
            res = res.view(np.ndarray)
        return res


    # check that an array is a valid Coords
    @staticmethod
    def __is_valid__(obj):
        return obj.shape[-1] == 3


    # overload some some ndarray methods that could easily return
    # invalid Coords
    def swapaxes(self,*args,**kargs):
        res = np.ndarray.swapaxes(self,*args,**kargs)
        if not Coords.__is_valid__(res):
            res = res.view(np.ndarray)
        return res

    ################ property methods ##########

    @property
    def xyz(self):
        """Returns the coordinates of the points as an ndarray.

        Returns an ndarray with shape `self.shape` except last axis
        is reduced to 2, providing a view on all the coordinates
        of all the points.

        """
        return self.view(type=ndarray)

    @property
    def x(self):
        """Returns the X-coordinates of all points.

        Returns an ndarray with shape `self.pshape()`, providing a view
        on the X-coordinates of all the points.

        """
        return self.xyz[..., 0]

    @property
    def y(self):
        """Returns the Y-coordinates of all points.

        Returns an ndarray with shape `self.pshape()`, providing a view
        on the Y-coordinates of all the points.

        """
        return self.xyz[..., 1]

    @property
    def z(self):
        """Returns the Z-coordinates of all points.

        Returns an ndarray with shape `self.pshape()`, providing a view
        on the Z-coordinates of all the points.

        """
        return self[..., 2].view(type=ndarray)

    @property
    def xy(self):
        """Returns the X- and Y-coordinates of all points.

        Returns an ndarray with shape `self.shape` except last axis
        is reduced to 2, providing a view on the X- and Y-coordinates
        of all the points.

        """
        return self.xyz[..., :2]

    @property
    def xz(self):
        """Returns the X- and Y-coordinates of all points.

        Returns an ndarray with shape `self.shape` except last axis
        is reduced to 2, providing a view on the X- and Z-coordinates
        of all the points.

        """
        return self.xyz[..., (0,2)]

    @property
    def yz(self):
        """Returns the X- and Y-coordinates of all points.

        Returns an ndarray with shape `self.shape` except last axis
        is reduced to 2, providing a view on the Y- and Z-coordinates
        of all the points.

        """
        return self.xyz[..., 1:]

    ################ property setters ##############

    @xyz.setter
    def xyz(self, value):
        """Set the XYZ coordinates of the points"""
        self[...] = value

    @x.setter
    def x(self, value):
        """Set the X coordinates of the points"""
        self[...,0] = value

    @y.setter
    def y(self, value):
        """Set the Y coordinates of the points"""
        self[...,1] = value

    @z.setter
    def z(self, value):
        """Set the Z coordinates of the points"""
        self[...,2] = value

    @xy.setter
    def xy(self, value):
        """Set the XY coordinates of the points"""
        self[...,:2] = value

    @xz.setter
    def xz(self, value):
        """Set the XZ coordinates of the points"""
        self[...,(0,2)] = value

    @yz.setter
    def yz(self, value):
        """Set the YZ coordinates of the points"""
        self[...,1:] = value

    @property
    def coords(self):
        """Returns the `Coords` object .

        This exists only for consistency with other classes.
        """
        return self

    ################ end property methods ##########


    def __repr__(self):
        """String representation of a Coords

        Examples
        --------
        >>> Coords([1.0,2.0,3.0])
        Coords([ 1.,  2.,  3.])

        """
        res = np.ndarray.__repr__(self)
        if self.dtype == Float:
            res = res.replace(', dtype=float32','')
        return res


    ## THIS IS A CANDIDATE FOR THE LIBRARY
    ## (possibly in a more general arrayprint form)
    ## (could be common with calpy)

    def fprint(self,fmt="%10.3e %10.3e %10.3e"):
        """Formatted printing of the points of a :class:`Coords` object.

        Parameters
        ----------
        fmt: string
            Format to be used to print a single point.
            The supplied format should contain exactly 3 formatting
            sequences, ome for each of the three coordinates.

        Examples
        --------
        >>> x = Coords([[[0.,0.],[1.,0.]],[[0.,1.],[0.,2.]]])
        >>> x.fprint()
        0.000e+00  0.000e+00  0.000e+00
        1.000e+00  0.000e+00  0.000e+00
        0.000e+00  1.000e+00  0.000e+00
        0.000e+00  2.000e+00  0.000e+00
        >>> x.fprint("%5.2f"*3)
        0.00 0.00 0.00
        1.00 0.00 0.00
        0.00 1.00 0.00
        0.00 2.00 0.00

        """
        for p in self.points():
            print(fmt % tuple(p))


###########################################################################
    #
    #   Methods that return information about a Coords object or other
    #   views on the object data, without changing the object itself.

    # General

    def pshape(self):
        """Return the points shape of the :class:`Coords` object.

        This is the shape of the :class:`numpy.ndarray` with the last axis
        removed.

        Note
        ----
        The full shape of the Coords array can be obtained from
        the inherited (NumPy) shape attribute.

        Examples
        --------
        >>> X = Coords(arange(12).reshape(2,1,2,3))
        >>> X.shape
        (2, 1, 2, 3)
        >>> X.pshape()
        (2, 1, 2)

        """
        return self.shape[:-1]


    def points(self):
        """Return the :class:`Coords` object as a flat set of points.

        Returns
        -------
        : Coords
            The Coords reshaped to a 2-dimensional array, flattening
            the structure of the points.

        Examples
        --------
        >>> X = Coords(arange(12).reshape(2,1,2,3))
        >>> X.shape
        (2, 1, 2, 3)
        >>> X.points().shape
        (4, 3)
        """
        return self.reshape((-1, 3))

    def npoints(self):
        """Return the total number of points in the Coords.

        Notes
        -----
        `npoints` and `ncoords` are equivalent. The latter exists to provide
        a common interface with other geometry classes.

        Examples
        --------
        >>> Coords(arange(12).reshape(2,1,2,3)).npoints()
        4
        """
        return array(self.shape[:-1]).prod()

    ncoords = npoints

    # Size, Bounds

    def bbox(self):
        """Return the bounding box of a set of points.

        The bounding box is the smallest rectangular volume in the global
        coordinates, such that no points of the :class:`Coords` are outside
        that volume.

        Returns
        -------
        : Coords (2,3)
            Coords array with two points: the first point contains the
            minimal coordinates, the second has the maximal ones.

        See Also
        --------
        center: return the center of the bounding box
        bboxPoint: return a corner or middle point of the bounding box
        bboxPoints: return all corners of the bounding box

        Examples
        --------
        >>> X = Coords([[0.,0.,0.],[3.,0.,0.],[0.,3.,0.]])
        >>> print(X.bbox())
        [[ 0. 0. 0.]
         [ 3. 3. 0.]]
        """
        if self.size > 0:
            x = self.points()
            bb = row_stack([ x.min(axis=0), x.max(axis=0) ])
        else:
            o = origin()
            bb = [o, o]
        return Coords(bb)


    def center(self):
        """Return the center of the :class:`Coords`.

        The center of a Coords is the center of its bbox().
        The return value is a (3,) shaped :class:`Coords` object.

        See also
        --------
        bbox: return the bounding box of the Coords
        centroid: return the average coordinates of the points

        Examples
        --------
        >>> X = Coords([[0.,0.,0.],[3.,0.,0.],[0.,3.,0.]])
        >>> print(X.center())
        [ 1.5 1.5 0. ]
        """
        X0, X1 = self.bbox()
        return 0.5 * (X0+X1)


    def bboxPoint(self, position):
        """Return a bounding box point of a Coords.

        Bounding box points are points whose coordinates are either
        the minimal value, the maximal value or the middle value
        for the Coords.
        Combining the three values in three dimensions results in
        3**3 = 27 alignment points. The corner points of the
        bounding box are a subset of these.

        Parameters
        ----------
        position: str
            String of three characters, one for each direction 0, 1, 2.
            Each character should be one of the following

                - '-': use the minimal value for that coordinate,
                - '+': use the minimal value for that coordinate,
                - '0': use the middle value for that coordinate.

            Any other character will set the corresponding coordinate to zero.

        Notes
        -----
        A string '000' is equivalent with center(). The values '---' and
        '+++' give the points of the bounding box.

        See Also
        --------
        Coords.align: translate Coords by bboxPoint

        Examples
        --------
        >>> X = Coords([[0.,0.,0.],[1.,1.,1.]])
        >>> print(X.bboxPoint('-0+'))
        [ 0.  0.5 1. ]

        """
        bb = self.bbox()
        al = { '-': bb[0], '+': bb[1], '0': 0.5*(bb[0]+bb[1]) }
        pt = zeros(3)
        for i, c in enumerate(position):
            if c in al:
                pt[i] = al[c][i]
        return Coords(pt)


    def bboxPoints(self):
        """Return all the corners of the bounding box point of a Coords.

        Returns
        -------
        : Coords (8,3)
            A Coords with the eight corners of the bounding box, in the order
            of a :py:attr:`elements.Hex8`.

        See also
        --------
        bbox: return only two points, with the minimum and maximum coordinates

        Examples
        --------
        >>> X = Coords([[0.,0.,0.],[3.,0.,0.],[0.,3.,0.]])
        >>> print(X.bboxPoints())
        [[ 0.  0.  0.]
         [ 3.  0.  0.]
         [ 3.  3.  0.]
         [ 0.  3.  0.]
         [ 0.  0.  0.]
         [ 3.  0.  0.]
         [ 3.  3.  0.]
         [ 0.  3.  0.]]
        """
        from pyformex.simple import cuboid
        return cuboid(*self.bbox()).coords.reshape(-1,3)


    def average(self,wts=None,axis=None):
        """Returns a (weighted) average of the :class:`Coords`.

        The average of a Coords is a Coords that is obtained by averaging
        the points along some or all axes.
        Weights can be specified to get a weighted average.

        Parameters
        ----------
        wts: float array-like, optional
            Weight to be attributed to the points. If provided, and `axis` is
            an int, `wts` should be 1-dim with the same length as the
            specified axis. Else, it has a shape equal to self.shape or
            self.shape[:-1].
        axis: int or tuple of ints, optional
            If provided, the average is computed along the specified
            axis/axes only. Else, the average is taken over all the points,
            thus over all the axes of the array except the last.

        Notes
        -----
        Averaging over the -1 axis does not make much sense.

        Examples
        --------
        >>> X = Coords([[[0.,0.,0.],[1.,0.,0.],[2.,0.,0.]], \
                [[4.,0.,0.],[5.,0.,0.],[6.,0.,0.]]])
        >>> X = Coords(arange(6).reshape(3,2,1))
        >>> X
        Coords([[[ 0.,  0.,  0.],
                 [ 1.,  0.,  0.]],
        <BLANKLINE>
                [[ 2.,  0.,  0.],
                 [ 3.,  0.,  0.]],
        <BLANKLINE>
                [[ 4.,  0.,  0.],
                 [ 5.,  0.,  0.]]])
        >>> print(X.average())
        [ 2.5  0.   0. ]
        >>> print(X.average(axis=0))
        [[ 2. 0. 0.]
         [ 3. 0. 0.]]
        >>> print(X.average(axis=1))
        [[ 0.5  0.   0. ]
         [ 2.5  0.   0. ]
         [ 4.5  0.   0. ]]
        >>> print(X.average(wts=[0.5,0.25,0.25],axis=0))
        [[ 1.5  0.   0. ]
         [ 2.5  0.   0. ]]
        >>> print(X.average(wts=[3,1],axis=1))
        [[ 0.25  0.    0.  ]
         [ 2.25  0.    0.  ]
         [ 4.25  0.    0.  ]]
        >>> print(X.average(wts=multiplex([3,1],3,0)))
        [ 2.25  0.    0.  ]
        """
        if axis == None:
            axis = tuple(range(self.ndim-1))
            if wts is not None:
                wts = np.asarray(wts)
                if wts.shape == self.shape:
                    pass
                elif wts.shape == self.shape[:-1]:
                    wts = multiplex(wts,3,-1)
                else:
                    raise ValueError("Shape of wts should match self.shape or self.pshape()")

        return np.average(self, weights=wts, axis=axis)


    def centroid(self):
        """Return the centroid of the :class:`Coords`.

        The centroid of Coords is the point whose coordinates
        are the mean values of all points.

        Returns
        -------
        : Coords (3,)
            A single point that is the centroid of the Coords.

        See also
        --------
        center: return the center of the bounding box.

        Examples
        --------
        >>> print(Coords([[0.,0.,0.],[3.,0.,0.],[0.,3.,0.]]).centroid())
        [ 1. 1. 0.]
        """
        return self.points().mean(axis=0)


    def centroids(self):
        """Return the Coords itself.

        Notes
        -----
        This method exists only to have a common interface with
        other geometry classes.
        """
        return self


    def sizes(self):
        """Return the bounding box sizes of the :class:`Coords`.

        Returns
        -------
        : array (3,)
            The length of the bounding box along the three global axes.

        See Also
        --------
        dsize: The diagonal size of the bounding box.
        principalSizes: the sizes of the bounding box along the principal axes

        Examples
        --------
        >>> print(Coords([[0.,0.,0.],[3.,0.,0.],[0.,3.,0.]]).sizes())
        [ 3. 3. 0.]

        """
        X0, X1 = self.bbox()
        return X1-X0


    def maxsize(self):
        """Return the maximum size of a Coords in any coordinate direction.

        Returns
        -------
        : float
            The maximum length of any edge of the bounding box.

        Notes
        -----
        This is a convenient shorthand for `self.sizes().max()`.

        See Also
        --------
        sizes: return the length of the bounding box along global axes
        bbox: return the bounding box

        Examples
        --------
        >>> print(Coords([[0.,0.,0.],[3.,0.,0.],[0.,3.,0.]]).maxsize())
        3.0

        """
        return self.sizes().max()


    def dsize(self):
        """Return the diagonal size of the bounding box of the :class:`Coords`.

        Returns
        -------
        : float
            The length of the diagonal of the bounding box.

        Notes
        -----
        All the points of the Coords are inside a sphere with the
        :meth:`center` as center and the :meth:`dsize` as
        length of the diameter (though it is not necessarily the
        smallest bouding sphere).
        :meth:`dsize` is in general a good estimate for the maximum size
        of the cross section to be expected when the object can be rotated
        freely around its center. It is conveniently used to zoom the
        camera on an object, while guaranteeing that the full object
        remains visible during rotations.

        See Also
        --------
        bsphere: return radius of smallest sphere encompassing all points
        sizes: return the length of the bounding box along global axes
        bbox: return the bounding box

        Examples
        --------
        >>> print(Coords([[0.,0.,0.],[3.,0.,0.],[0.,3.,0.]]).dsize())
        4.24264

        """
        X0, X1 = self.bbox()
        return length(X1-X0)


    def bsphere(self):
        """Return the radius of the bounding sphere of the :class:`Coords`.

        The bounding sphere used here is the smallest sphere with center
        in the center() of the :class:`Coords`, and such that no points of the
        Coords are lying outside the sphere.

        Returns
        -------
        : float
            The maximum distance of any point to the `Coords.center`.

        Notes
        -----
        This is not necessarily the absolute smallest bounding sphere,
        because we use the center from lookin only in the global axes
        directions.

        Examples
        --------
        >>> X = Coords([[0.,0.,0.],[3.,0.,0.],[0.,3.,0.]])
        >>> print(X.dsize(), X.bsphere())
        4.24264 2.12132
        >>> X = Coords([[0.5,0.],[1.,0.5],[0.5,1.0],[0.0,0.5]])
        >>> print(X.dsize(), X.bsphere())
        1.41421 0.5

        """
        return self.distanceFromPoint(self.center()).max()


    def bboxes(self):
        """Return the bboxes of all subsets of points in the Coords.

        Subsets of points are 2-dim subarrays of the Coords, taken
        along the two last axes. If the Coords has ndim==2, there
        is only one subset: the full Coords.

        Returns
        -------
        float array
            Array with shape (...,2,3). The elements along the penultimate
            axis are the minimal and maximal values of the Coords
            along that axis.

        Examples
        --------
        >>> X = Coords(arange(18).reshape(2,3,3))
        >>> print(X)
         [[[  0.   1.   2.]
          [  3.   4.   5.]
          [  6.   7.   8.]]
        <BLANKLINE>
         [[  9.  10.  11.]
          [ 12.  13.  14.]
          [ 15.  16.  17.]]]
        >>> print(X.bboxes())
        [[[  0.   1.   2.]
          [  6.   7.   8.]]
        <BLANKLINE>
         [[  9.  10.  11.]
          [ 15.  16.  17.]]]
        """
        return minmax(self, axis=1)


    # Inertia

    @utils.warning("warn_inertia_changed")
    def inertia(self,mass=None):
        """Return inertia related quantities of the :class:`Coords`.

        Parameters
        ----------
        mass: float array, optional
            If provided, it is a 1-dim array with :meth:`npoints` weight
            values for the points, in the order of the :meth:`points`.
            The default is to attribute a weight 1.0 to each point.

        Returns
        -------
        :class:`~inertia.Inertia`
            The Inertia object has the following attributes:

            - ``mass``: the total mass (float)
            - ``ctr``: the center of mass: float (3,)
            - ``tensor``: the inertia tensor in the central axes: shape (3,3)

        See Also
        --------
        principalCS: Return the principal axes of the inertia tensor

        Examples
        --------
        >>> from pyformex.elements import Tet4
        >>> I = Tet4.vertices.inertia()
        >>> print(I.tensor)
        [[ 1.5  0.25 0.25]
         [ 0.25 1.5  0.25]
         [ 0.25 0.25 1.5 ]]
        >>> print(I.ctr)
        [ 0.25 0.25 0.25]
        >>> print(I.mass)
        4.0

        """
        from pyformex import inertia
        M,C,I = inertia.point_inertia(self.points(), mass)
        I = inertia.Tensor(I)
        return inertia.Inertia(I,ctr=C,mass=M)


    def principalCS(self,mass=None):
        """Return a CoordSys formed by the principal axes of inertia.

        Parameters
        ----------
        mass: 1-dim float array (:meth:`points`,), optional
            The mass to be attributed to each of the points, in the order
            of :meth:`npoints`. If not provided, a mass 1.0 will be
            attributed to each point.

        Returns
        -------
        :class:`~coordsys.CoordSys` object.
            Coordinate system aligned along the principal axes of the inertia,
            for the specified point masses. The origin of the CoordSys is
            the center of mass of the Coords.

        See Also
        --------
        centralCS: CoordSys at the center of mass, but axes along global
            directions

        Examples
        --------
        >>> from pyformex.elements import Tet4
        >>> print(Tet4.vertices.principalCS())
        CoordSys: trl=[ 0.25  0.25  0.25]; rot=[[ 0.58  0.58  0.58]
                                                [ 0.34 -0.81  0.47]
                                                [ 0.82 -0.41 -0.41]]

        """
        from pyformex.coordsys import CoordSys
        I = self.inertia(mass)
        prin, axes = I.principal()
        return CoordSys(rot=axes,trl=I.ctr)


    def principalSizes(self):
        """Return the sizes in the principal directions of the :class:`Coords`.

        Returns
        -------
        float array (3,)
            Array with the size of the bounding box along the 3
            principal axes.

        Notes
        -----
        This is a convenient shorthand for:
        ``self.toCS(self.principalCS()).sizes()``

        Examples
        --------
        >>> print(Coords([[[0.,0.,0.],[3.,0.,0.]]]).rotate(30,2).principalSizes())
        [ 0. 0. 3.]

        """
        return self.toCS(self.principalCS()).sizes()


    def centralCS(self,mass=None):
        """Returns the central coordinate system of the Coords.

        Parameters
        ----------
        mass: 1-dim float array (:meth:`points`,), optional
            The mass to be attributed to each of the points, in the order
            of :meth:`npoints`. If not provided, a mass 1.0 will be
            attributed to each point.

        Returns
        -------
        :class:`~coordsys.CoordSys` object.
            Coordinate system with origin at the center of mass of the
            Coords and axes parallel to the global axes.

        See Also
        --------
        principalCS: CoordSys aligned with principa axes of inertia tensor

        Examples
        --------
        >>> from pyformex.elements import Tet4
        >>> print(Tet4.vertices.centralCS())
        CoordSys: trl=[ 0.25  0.25  0.25]; rot=[[ 1.  0.  0.]
                                                [ 0.  1.  0.]
                                                [ 0.  0.  1.]]

        """
        from pyformex.coordsys import CoordSys
        C = self.reshape(-1,3).average(wts=mass,axis=0)
        return CoordSys(trl=C)


    #  Distance


    def distanceFromPoint(self, p):
        """Returns the distance of all points from the point p.

        Parameters
        ----------
        p: float array-like with shape (3,) or (1,3)
            Coordinates of a single point in space

        Returns
        -------
        float array
            Array with shape :meth:`pshape` holding the distance of
            each point to point p. All values are positive or zero.

        See Also
        --------
        closestPoint: return the point of Coords closest to given point

        Examples
        --------
        >>> X = Coords([[0.,0.,0.],[2.,0.,0.],[1.,3.,0.],[-1.,0.,0.]])
        >>> print(X.distanceFromPoint([0.,0.,0.]))
        [ 0.    2.    3.16  1.  ]

        """
        p = checkArray(p,size=3,kind='f',allow='i').reshape(3)
        return length(self-p)


    def distanceFromLine(self, p, n):
        """Returns the distance of all points from the line (p,n).

        Parameters
        ----------
        p: float array-like with shape (3,) or (1,3)
            Coordinates of some point on the line.
        n: float array-like with shape (3,) or (1,3)
            Vector specifying the direction of the line.

        Returns
        -------
        float array
            Array with shape :meth:`pshape` holding the distance
            of each point to the line through p and having direction n.
            All values are positive or zero.

        Examples
        --------
        >>> X = Coords([[0.,0.,0.],[2.,0.,0.],[1.,3.,0.],[-1.,0.,0.]])
        >>> print(X.distanceFromLine([0.,0.,0.],[1.,1.,0.]))
        [ 0.   1.41 1.41 0.71]

        """
        p = checkArray(p,size=3,kind='f',allow='i').reshape(3)
        n = checkArray(n,size=3,kind='f',allow='i').reshape(3)
        n = normalize(n)
        xp = self-p
        xpt = dotpr(xp,n)
        a = dotpr(xp,xp)-xpt*xpt
        return sqrt(a.clip(0))


    def distanceFromPlane(self, p, n):
        """Return the distance of all points from the plane (p,n).

        Parameters
        ----------
        p: float array-like with shape (3,) or (1,3)
            Coordinates of some point in the plane.
        n: float array-like with shape (3,) or (1,3)
            The normal vector to the plane.

        Returns
        -------
        float array
            Array with shape :meth:`pshape` holding the distance
            of each point to the plane through p and having normal n.
            The values are positive if the point is on the side of the
            plane indicated by the positive normal.

        See Also
        --------
        directionalSize: find the most distant points at both sides of plane

        Examples
        --------
        >>> X = Coords([[0.,0.,0.],[2.,0.,0.],[1.,3.,0.],[-1.,0.,0.]])
        >>> print(X.distanceFromPlane([0.,0.,0.],[1.,0.,0.]))
        [ 0.  2.  1. -1.]

        """
        p = checkArray(p,size=3,kind='f',allow='i').reshape(3)
        n = checkArray(n,size=3,kind='f',allow='i').reshape(3)
        n = normalize(n)
        d = inner(self, n) - inner(p, n)
        return asarray(d)


    def closestToPoint(self, p, return_dist=False):
        """Returns the point closest to a given point p.

        Parameters
        ----------
        p: array-like (3,)
            Coordinates of a single point in space

        Returns
        -------
        :int
            Index of the point in the Coords that has the
            minimal Euclidean distance to the point `p`.
            Use this index with self.points() to get the
            coordinates of that point.

        Examples
        --------
        >>> X = Coords([[[0.,0.,0.],[3.,0.,0.],[0.,3.,0.]]])
        >>> X.closestToPoint([2.,0.,0.])
        1
        >>> X.closestToPoint([2.,0.,0.],True)
        (1, 1.0)

        """
        d = self.distanceFromPoint(p)
        i = d.argmin()
        if return_dist:
            return i,d.flat[i]
        else:
            return i


    def directionalSize(self,n,p=None,return_points=False):
        """Returns the extreme distances from the plane p,n.

        Parameters
        ----------
        n: a single int or a float array-like (3,)
            The direction of the normal to the plane. If an int, it is the
            number of a global axis. Else it is a vector with 3 components.
        p: array-like (3,), optional
            Coordinates of a point in the plane. If not provided,
            the :meth:`center` of the Coords is used.
        return_points: bool
           If True, also return a Coords with two points along the line
           (p,n) and at the extreme distances from the plane(p,n).

        Returns
        -------
        dmin: float
            The minimal (signed) distance of a point of the Coords to the
            plane (p,n). The value can be negative or positive.
        dmax: float
            The maximal (signed) distance of a point of the Coords to the
            plane (p,n). The value can be negative or positive.
        points: Coords (2,3), optional
            If `return_points=True` is provided, also returns a Coords
            holding two points on the line (p,n) with minimal and
            maximal distance from the plane (p,n). These two points
            together with the normal `n` define two parallel planes
            such that all points of `self` are between or on the planes.

        Notes
        -----
        The maximal size of `self` in the direction `n` is found from
        the difference `dmax` - dmin`. See also :meth:`directionalWidth`.

        See also
        --------
        directionalExtremes: return two points in the extreme planes
        directionalWidth: return the distance between the extreme planes
        distanceFromPlane: return distance of all points to a plane

        Examples
        --------
        >>> X = Coords([[[0.,0.,0.],[3.,0.,0.],[0.,3.,0.]]])
        >>> X.directionalSize([1,0,0])
        (-1.5, 1.5)
        >>> X.directionalSize([1,0,0],[1.,0.,0.])
        (-1.0, 2.0)
        >>> X.directionalSize([1,0,0],return_points=True)
        (-1.5, 1.5, Coords([[ 0. ,  1.5,  0. ],
                [ 3. ,  1.5,  0. ]]))

        """
        n = unitVector(n)
        if p is None:
            p = self.center()
        else:
            p = Coords(p)
        d = self.distanceFromPlane(p, n)
        dmin, dmax = d.min(), d.max()

        if return_points:
            return dmin, dmax, Coords([p+dmin*n, p+dmax*n])
        else:
            return dmin, dmax


    def directionalExtremes(self,n,p=None):
        """Returns extremal planes in the direction n.

        Parameters: see :meth:`directionalSize`.

        Returns
        -------
        :Coords (2,3)
            A Coords holding the two points on the line (p,n) with minimal
            and maximal distance from the plane (p,n). These two points
            together with the normal `n` define two parallel planes
            such that all points of `self` are between or on the planes.

        See also
        --------
        directionalSize: return minimal and maximal distance from plane

        Notes
        -----
        This is like directionalSize with the return_points options,
        but only returns the extreme points.

        Examples
        --------
        >>> X = Coords([[[0.,0.,0.],[3.,0.,0.],[0.,3.,0.]]])
        >>> X.directionalExtremes([1,0,0])
        Coords([[ 0. ,  1.5,  0. ],
                [ 3. ,  1.5,  0. ]])
        """
        return self.directionalSize(n,p,return_points=True)[2]


    def directionalWidth(self, n):
        """Returns the width of a Coords in the given direction.

        Parameters: see :meth:`directionalSize`.

        Returns
        -------
        :float
            The size of the Coords in the direction `n`. This is the
            distance between the extreme planes with normal `n` touching
            the Coords.

        See also
        --------
        directionalSize: return minimal and maximal distance from plane

        Notes
        -----
        This is like directionalSize but only returns the difference between
        `dmax` and `dmin`.

        Examples
        --------
        >>> X = Coords([[[0.,0.,0.],[3.,0.,0.],[0.,3.,0.]]])
        >>> print(X.directionalWidth([1,0,0]))
        3.0
        """
        dmin, dmax = self.directionalSize(n)
        return dmax-dmin


    # Test position

    def test(self,dir=0,min=None,max=None,atol=0.):
        """Flag points having coordinates between min and max.

        Test the position of the points of the :class:`Coords` with respect to
        one or two parallel planes. This method is very convenient in clipping
        a Coords in a specified direction. In most cases the clipping
        direction is one of the global coordinate axes, but a general
        direction may be used as well.

        Testing along global axis directions is highly efficient. It tests
        whether the corresponding coordinate is above or equal to the `min`
        value and/or below or equal to the `max` value. Testing in a general
        direction tests whether the distance to the `min` plane is positive
        or zero and/or the distance to the `max` plane is negative or zero.

        Parameters
        ----------
        dir: a single int or a float array-like (3,)
            The direction in which to measure distances. If an int, it is
            one of the global axes (0,1,2). Else it is a vector with 3
            components. The default direction is the global x-axis.
        min: float or point-like, optional
            Position of the minimal clipping plane.
            If `dir` is an int, this is a single float giving the coordinate
            along the specified global axis. If `dir` is a vector, this must
            be a point and the minimal clipping plane is defined by this point
            and the normal vector `dir`. If not provided, there is no clipping
            at the minimal side.
        max: float or point-like.
            Position of the maximal clipping plane.
            If `dir` is an int, this is a single float giving the coordinate
            along the specified global axis. If `dir` is a vector, this must
            be a point and the maximal clipping plane is defined by this point
            and the normal vector `dir`. If not provided, there is no clipping
            at the maximal side.
        atol: float
            Tolerance value added to the tests to account for accuracy
            and rounding errors.
            A `min` test will be ok if the point's distance from the
            `min` clipping plane is `> -atol` and/or the distance from the
            `max` clipping plane is `< atol`. Thus a positive atol widens the
            clipping planes.

        Returns
        -------
        : bool array with shape :meth:`pshape`
            Array flagging whether the points for the Coords pass the test(s)
            or not. The return value can directly be used as an index to
            `self` to obtain a :class:`Coords` with the points satisfying
            the test (or not).

        Raises
        ------
        ValueError: At least one of min or max have to be specified
            If neither `min` nor `max` are provided.

        Examples
        --------
        >>> x = Coords([[[0.,0.],[1.,0.]],[[0.,1.],[0.,2.]]])
        >>> print(x.test(min=0.5))
        [[False  True]
         [False False]]
        >>> t = x.test(dir=1,min=0.5,max=1.5)
        >>> print(x[t])
        [[ 0. 1. 0.]]
        >>> print(x[~t])
        [[ 0. 0. 0.]
         [ 1. 0. 0.]
         [ 0. 2. 0.]]

        """
        if min is None and max is None:
            raise ValueError("At least one of min or max have to be specified.")

        if array(dir).size == 1:
            if min is not None:
                T1 = self[..., dir] >= (min - atol)
            if max is not None:
                T2 = self[..., dir] <= (max + atol)
        else:
            if min is not None:
                T1 = self.distanceFromPlane(min, dir) > - atol
            if max is not None:
                T2 = self.distanceFromPlane(max, dir) < atol

        if min is None:
            T = T2
        elif max is None:
            T = T1
        else:
            T = T1 * T2
        return asarray(T)



##############################################################################

    def set(self, f):
        """Set the coordinates from those in the given array.

        Parameters
        ----------
        f: float array-like, broadcastable to self.shape.
            The coordinates to replace the current ones. This can not
            be used to chage the shape of the Coords.

        Raises
        ------
        ValueError:
            If the shape of `f` does not allow broadcasting to `self.shape`.

        Examples
        --------
        >>> x = Coords([[0],[1],[2]])
        >>> print(x)
        [[ 0.  0.  0.]
         [ 1.  0.  0.]
         [ 2.  0.  0.]]
        >>> x.set([0.,1.,0.])
        >>> print(x)
        [[ 0.  1.  0.]
         [ 0.  1.  0.]
         [ 0.  1.  0.]]

        """
        f = checkArray(f,kind='f',allow='i')
        if checkBroadcast(self.shape,f.shape) != self.shape:
            raise ValueError("Invalid array shape")
        self[...] = f      # do not be tempted to use self = f !

##############################################################################
    #
    #   Transformations that preserve the topology (but change coordinates)
    #
    #   A. Affine transformations
    #
    #      Scaling
    #      Translation
    #      Central Dilatation = Scaling + Translation
    #      Rotation
    #      Shear
    #      Reflection
    #      Affine
    #
    #  The following methods return transformed coordinates, but by default
    #  they do not change the original data. If the optional argument inplace
    #  is set True, however, the coordinates are changed inplace.


    def scale(self,scale,dir=None,center=None,inplace=False):
        """Return a scaled copy of the :class:`Coords` object.

        Parameters
        ----------
        scale: float or tuple of 3 floats
            Scaling factor(s). If it is a single value, and no `dir` is
            provided, scaling is uniformly applied to all axes; if `dir`
            is provided, only to the specified directions. If it is a tuple,
            the three scaling factors are applied to the respective
            global axes.
        dir: int or tuple of ints, optional
            One or more global axis numbers (0,1,2), indicating the direction(s)
            that should be scaled with the (single) value `scale`.
        center: point-like, optional
            If provided, use this point as the center of the scaling.
            The default is the global origin.
        inplace: bool,optional
            If True, the coordinates are change in-place.

        Returns
        -------
        : Coords
            The Coords scaled as specified.

        Notes
        -----
        If a `center` is provided,the operation is equivalent with
        ``self.translate(-center).scale(scale,dir).translate(center)``

        Examples
        --------
        >>> X = Coords([1.,1.,1.])
        >>> print(X.scale(2))
        [ 2. 2. 2.]
        >>> print(X.scale([2,3,4]))
        [ 2. 3. 4.]
        >>> print(X.scale(2,dir=(1,2)).scale(4,dir=0))
        [ 4. 2. 2.]
        >>> print(X.scale(2,center=[1.,0.5,0.]))
        [ 1.   1.5  2. ]

        """
        if center is not None:
            center = asarray(center)
            return self.trl(-center).scale(scale, dir).translate(center)

        if inplace:
            out = self
        else:
            out = self.copy()
        if dir is None:
            out *= scale
        else:
            out[..., dir] *= scale
        return out


    def translate(self,dir,step=1.0,inplace=False):
        """Return a translated copy of the :class:`Coords` object.

        Translate the Coords in the direction `dir` over a distance
        `step * length(dir)`.

        Parameters
        ----------
        dir: int (0,1,2) or float array-like (...,3)
            The translation vector. If an int, it specifies a global axis
            and the translation is in the direction of that axis.
            If an array-like, it specifies one or more translation vectors.
            If more than one, the array should be broadcastable to the Coords
            shape: this allows to translate different parts of the Coords over
            different vectors, all in one operation.
        step: float
            If ``dir`` is an int, this is the length of the translation.
            Else, it is a multiplying factor applied to ``dir`` prior to
            applying the translation.

        Returns
        -------
        : Coords
            The Coords translated over the specified vector(s).

        Note
        ----
        :meth:`trl` is a convenient shorthand for :meth:`translate`.

        See Also
        --------
        centered: translate to center around origin
        Coords.align: translate to align bounding box

        Examples
        --------
        >>> x = Coords([1.,1.,1.])
        >>> print(x.translate(1))
        [ 1. 2. 1.]
        >>> print(x.translate(1,1.))
        [ 1. 2. 1.]
        >>> print(x.translate([0,1,0]))
        [ 1. 2. 1.]
        >>> print(x.translate([0,2,0],0.5))
        [ 1. 2. 1.]
        >>> x = Coords(arange(4).reshape(2,2,1))
        >>> x
        Coords([[[ 0.,  0.,  0.],
                [ 1.,  0.,  0.]],
        <BLANKLINE>
               [[ 2.,  0.,  0.],
                [ 3.,  0.,  0.]]])
        >>> x.translate([[10.,-5.,0.],[20.,4.,0.]]) # translate with broadcasting
        Coords([[[ 10.,  -5.,   0.],
                [ 21.,   4.,   0.]],
        <BLANKLINE>
               [[ 12.,  -5.,   0.],
                [ 23.,   4.,   0.]]])

        """
        if inplace:
            out = self
        else:
            out = self.copy()
        if isInt(dir):
            out[...,dir] += step
        else:
            dir = Coords(dir, copy=True)
            if step != 1.:
                dir *= step
            out += dir
        return out


    def centered(self):
        """Return a centered copy of the Coords.

        Returns
        -------
        : Coords
            The Coords translated over tus that its :meth:`center`
            coincides with the origin of the global axes.

        Notes
        -----
        This is equivalent with ``self.translate(-self.center())``

        Examples
        --------
        >>> X = Coords('0123')
        >>> print(X)
        [[ 0.  0.  0.]
         [ 1.  0.  0.]
         [ 1.  1.  0.]
         [ 0.  1.  0.]]
        >>> print(X.centered())
        [[-0.5 -0.5  0. ]
         [ 0.5 -0.5  0. ]
         [ 0.5  0.5  0. ]
         [-0.5  0.5  0. ]]
       """
        return self.translate(-self.center())


    def align(self,alignment='---',point=[0., 0., 0.]):
        """Align a :class:`Coords` object on a given point.

        Alignment involves a translation such that the bounding box
        of the Coords object becomes aligned with a given point.
        The bounding box alignment is done by the translation of a
        to the target point.

        Parameters
        ----------
        alignment: str
            The requested alignment is a string of three characters,
            one for each of the coordinate axes. The character determines how
            the structure is aligned in the corresponding direction:

                - '-': aligned on the minimal value of the bounding box,
                - '+': aligned on the maximal value of the bounding box,
                - '0': aligned on the middle value of the bounding box.

            Any other value will make the alignment in that direction unchanged.

        point: point-like
            The target point of the alignment.

        Returns
        -------
        : Coords
            The Coords translated thus that the `alignment`
            :meth:`bboxPoint` is at `point`.

        Notes
        -----
        The default parameters translate the Coords thus that all points
        are in the octant with all positive coordinate values.

        ``Coords.align(alignment = '000')`` will center the object around
        the origin, just like the :meth:`centered` (which is slightly faster).
        This can however be used for centering around any point.

        See also
        --------
        align: aligning multiple objects with respect to each other.
        """
        return self.translate(point-self.bboxPoint(alignment))


    def rotate(self,angle,axis=2,around=None,angle_spec=DEG):
        """Return a copy rotated over angle around axis.

        Parameters
        ----------
        angle: float or float array-like (3,3)
            If a float, it is the rotation angle, by default in degrees,
            and the parameters (angle, axis, angel_spec) are passed to
            :func:`~arraytools.rotationMatrix` to produce a (3,3) rotation
            matrix.
            Alternatively, the rotation matrix may be directly provided
            in the `angle` parameter. The `axis` and `angle_spec` are then
            ignored.
        axis: int (0,1,2) or float array-like (3,)
            Only used if `angle` is a float.
            If provided, it specifies the direction of the rotation axis:
            either one of 0,1,2 for a global axis, or a vector with
            3 components for a general direction. The default (axis 2)
            is convenient for working with 2D-structures in the x-y plane.
        around: float array-like (3,)
            If provided, it species a point on the rotation axis. If not,
            the rotation axis goes through the origin of the global axes.
        angle_spec: float, DEG or RAD, optional
            Only used if `angle` is a float.
            The default (DEG) interpretes the angle in degrees. Use RAD to
            specify the angle in radians.

        Returns
        -------
        Coords
            The Coords rotated as specified by the parameters.

        Note
        ----
        :meth:`rot` is a convenient shorthand for :meth:`rotate`.

        See Also
        --------
        translate: translate a Coords
        affine: rotate and translate a Coords
        arraytools.rotationMatrix:
             create a rotation matrix for use in :meth:`rotate`

        Examples
        --------
        >>> X = Coords('0123')
        >>> print(X.rotate(30))
        [[ 0.    0.    0.  ]
         [ 0.87  0.5   0.  ]
         [ 0.37  1.37  0.  ]
         [-0.5   0.87  0.  ]]
        >>> print(X.rotate(30,axis=0))
        [[ 0.    0.    0.  ]
         [ 1.    0.    0.  ]
         [ 1.    0.87  0.5 ]
         [ 0.    0.87  0.5 ]]
        >>> print(X.rotate(30,axis=0,around=[0.,0.5,0.]))
        [[ 0.    0.07 -0.25]
         [ 1.    0.07 -0.25]
         [ 1.    0.93  0.25]
         [ 0.    0.93  0.25]]
        >>> m = rotationMatrix(30,axis=0)
        >>> print(X.rotate(m))
        [[ 0.    0.    0.  ]
         [ 1.    0.    0.  ]
         [ 1.    0.87  0.5 ]
         [ 0.    0.87  0.5 ]]
        """
        mat = asarray(angle)
        if mat.size == 1:
            mat = rotationMatrix(angle, axis=axis, angle_spec=angle_spec)
        mat = checkArray(mat,shape=(3,3),kind='f')
        if around is not None:
            around = asarray(around)
            out = self.translate(-around)
        else:
            out = self
        return out.affine(mat, around)


    def shear(self,dir,dir1,skew,inplace=False):
        """Return a copy skewed in the direction of a global axis.

        This translates points in the direction of a global axis,
        over a distance dependent on the coordinates along another axis.

        Parameters
        ----------
        dir: int (0,1,2)
            Global axis in which direction the points are translated.
        dir1: int (0,1,2)
            Global axis whose coordinates determine the length of the
            translation.
        skew: float
            Multiplication factor to the coordinates dir1 defining the
            translation distance.
        inplace: bool, optional
            If True, the coordinates are translated in-place.

        Notes
        -----
        This replaces the  coordinate ``dir`` with ``(dir + skew * dir1)``.
        If dir and dir1 are different, rectangular shapes in the plane
        (dir,dir1) are thus skewed along the direction dir into
        parallellogram shapes. If dir and dir1 are the same direction, the
        effect is that of scaling in the dir direction.

        Examples
        --------
        >>> X = Coords('0123')
        >>> print(X.shear(0,1,0.5))
        [[ 0.   0.   0. ]
         [ 1.   0.   0. ]
         [ 1.5  1.   0. ]
         [ 0.5  1.   0. ]]
        """
        if inplace:
            out = self
        else:
            out = self.copy()
        out[..., dir] += skew * out[..., dir1]
        return out


    # TODO: Add mirroring against any plane/axis/point
    #
    def reflect(self,dir=0,pos=0.,inplace=False):
        """Reflect the coordinates in the direction of a global axis.

        Parameters
        ----------
        dir: int (0,1,2)
            Global axis direction of the reflection (default 0 or x-axis).
        pos: float
            Offset of the mirror plane from origin (default 0.0)
        inplace: bool, optional
            If True, the coordinates are translated in-place.

        Returns
        -------
        :Coords
            A mirror copy with respect to the plane perpendicular to axis `dir`
            and placed at coordinate `pos` along the `dir` axis.

        Examples
        --------
        >>> X = Coords('012')
        >>> print(X)
        [[ 0.  0.  0.]
         [ 1.  0.  0.]
         [ 1.  1.  0.]]
        >>> print(X.reflect(0))
        [[ 0.  0.  0.]
         [-1.  0.  0.]
         [-1.  1.  0.]]
        >>> print(X.reflect(1,0.5))
        [[ 0.  1.  0.]
         [ 1.  1.  0.]
         [ 1.  0.  0.]]
        >>> print(X.reflect([0,1],[0.5,0.]))
        [[ 1.  0.  0.]
         [ 0.  0.  0.]
         [ 0. -1.  0.]]

        """
        if inplace:
            out = self
        else:
            out = self.copy()
        out[..., dir] = 2*asarray(pos) - out[..., dir]
        return out


    def affine(self,mat,vec=None):
        """Perform a general affine transformation.

        Parameters
        ----------
        mat: float array-like (3,3)
            Matrix used in post-multiplication on a row vector to produce
            a new vector. The matrix can express scaling and/or rotation
            or a more general (affine) transformation.
        vec: float array-like (3,)
            Translation vector to add after the transformation with `mat`.

        Returns
        -------
        : Coords
            A Coords with same shape as self, but with coordinates given by
            ``self * mat + vec``. If `mat` is a rotation matrix or a uniform
            scaling plus rotation, the full operation performs a rigid rotation
            plus translation of the object.

        Examples
        --------
        >>> X = Coords('0123')
        >>> S = array([[2.,0.,0.],[0.,3.,0.],[0.,0.,4.]]) # non-uniform scaling
        >>> R = rotationMatrix(90.,2) # rotation matrix
        >>> T = [20., 0., 2.] # translation
        >>> M = dot(S,R)  # combined scaling and rotation
        >>> print(X.affine(M,T))
        [[ 20.   0.   2.]
         [ 20.   2.   2.]
         [ 17.   2.   2.]
         [ 17.   0.   2.]]

        """
        out = dot(self, mat)
        if vec is not None:
            out += vec
        return out


    def toCS(self,cs):
        """Transform the coordinates to another CoordSys.

        Parameters
        ----------
        cs: :class:`~coordsys.CoordSys` object
            Cartesian coordinate system in which to take the coordinates of
            the current Coords object.

        Returns
        -------
        Coords
            A Coords object identical to self but having global coordinates
            equal to the coordinates of self in the `cs` CoordSys axes.

        Note
        ----
        This returns the coordinates of the original points in another
        CoordSys. If you use these coordinates as points in the global axes,
        the transformation of the original points to these new ones is the
        inverse transformation of the transformation of the global axes
        to the `cs` coordinate system.

        See Also
        --------
        fromCS: the inverse transformation

        Examples
        --------
        >>> X = Coords('01')
        >>> print(X)
        [[ 0.  0.  0.]
         [ 1.  0.  0.]]
        >>> from pyformex.coordsys import CoordSys
        >>> CS = CoordSys(oab=[[0.5,0.,0.],[1.,0.5,0.],[0.,1.,0.]])
        >>> print(CS)
        CoordSys: trl=[ 0.5  0.   0. ]; rot=[[ 0.71  0.71  0.  ]
                                             [-0.71  0.71  0.  ]
                                             [ 0.   -0.    1.  ]]
        >>> print(X.toCS(CS))
        [[-0.35  0.35  0.  ]
         [ 0.35 -0.35  0.  ]]
        >>> print(X.toCS(CS).fromCS(CS))
        [[ 0.  0.  0.]
         [ 1. -0.  0.]]

        """
        return self.trl(-cs.trl).rot(cs.rot.transpose())


    def fromCS(self,cs):
        """Transform the coordinates from another CoordSys to global axes.

        Parameters
        ----------
        cs: :class:`~coordsys.CoordSys` object
            Cartesian coordinate system in which the current coordinate
            values are taken.

        Returns
        -------
        Coords
            A Coords object with the global coordinates of the same points
            as the input coordinates represented in the `cs` CoordSys axes.

        See Also
        --------
        toCS: the inverse transformation


        Examples: see :meth:`toCS`
        """
        return self.rot(cs.rot).trl(cs.trl)


    def transformCS(self,cs,cs0=None):
        """Perform a coordinate system transformation on the Coords.

        This method transforms the Coords object by the transformation that
        turns one coordinate system into a another.

        Parameters
        ----------
        cs: :class:`~coordsys.CoordSys`
            The final coordinate system.
        cs0: :class:`~coordsys.CoordSys`, optional
            The initial coordinate system. If not provided, the global
            coordinate system is used.

        Returns
        -------
        Coords
            The input Coords transformed by the same affine transformation
            that turns the axes of the coordinate system `cs0`
            into those of the system `cs`.

        Notes
        -----
        For example, with the default `cs0` and a `cs` CoordSys created with
        the points ::

           0. 1. 0.
          -1. 0. 0.
           0. 0. 1.
           0. 0. 0.

        the transformCS results in a rotation of 90 degrees around the z-axis.

        See Also
        --------
        toCS: transform coordinates to another CS
        fromCS: transfrom coordinates from another CS

        """
        if cs0 is not None:
            f = self.toCS(cs0)
        else:
            f = self
        return f.fromCS(cs)


    def position(self, x, y):
        """Position a :class:`Coords` so that 3 points x are aligned with y.

        Aligning 3 points x with 3 points y is done by a rotation and
        translation in such way that

        - point x0 coincides with point y0,
        - line x0,x1 coincides with line y0,y1
        - plane x0,x1,x2 coincides with plane y0,y1,y2

        Parameters
        ----------
        x: float array-like (3,3)
            Original coordinates of three non-collinear points. These points
            can be be part of the Coords or not.
        y: float array-like (3,3)
            Final coordinates of the three points.

        Returns
        -------
        Coords
            The input Coords rotated and translated thus that the points
            x are aligned with y.

        Notes
        -----
        This is a convenient shorthand for ``self.affine(*trfmat(x, y))``.

        See Also
        --------
        arraytools.trfmat: compute the transformation matrices from points x to y
        affine: general transform using rotation and translation

        Examples
        --------
        >>> X = Coords([[0,0,0],[1,0,0],[1,1,0]])
        >>> Y = Coords([[1,1,1],[1,10,1],[1,1,100]])
        >>> print(X.position(X,Y))
        [[ 1.  1.  1.]
         [ 1.  2.  1.]
         [ 1.  2.  2.]]
        """
        return self.affine(*trfmat(x, y))

#
#
#   B. Non-Affine transformations.
#
#      These always return copies !
#
#        Cylindrical, Spherical, Isoparametric
#

    def cylindrical(self,dir=(0, 1, 2),scale=(1., 1., 1.),angle_spec=DEG):
        """Convert from cylindrical coordinates to cartesian.

        A cylindrical coordinate system is defined by a longitudinal axis
        axis (z) and a radial axis (r). The cylindrical coordinates of
        a point are:

        - r: the radial distance from the z-axis,
        - theta: the circumferential angle measured positively around
          the z-axis starting from zero at the (r-z) halfplane,
        - z: the axial distance along the z-axis,

        This function interpretes the 3 coordinates of the points
        as (r,theta,z) values and computes the corresponding global
        cartesian coordinates (x,y,z).

        Parameters
        ----------
        dir: tuple of 3 ints, optional
            If provided, it is a permutation of (0,1,2) and specifies
            which of the current coordinates are interpreted as resp.
            distance(r), angle(theta) and height(z).
            Default order is (r,theta,z).
            Beware that if the permutation is not conserving the order
            of the axes, a left-handed system results, and the Coords will
            appear mirrored in the right-handed systems exclusively used
            by pyFormex
        scale: tuple of 3 floats, optional
            Scaling factors that are applied on the values prior to make the
            conversion from cylindrical to cartesian coordinates. These factors
            are always given in the order (r,theta,z), irrespective of the
            permutation by `dir`.
        angle_spec: float, DEG or RAD, optional
            Multiplication factor for angle coordinates.
            The default (DEG) interpretes the angle in degrees. Use RAD to
            specify the angle in radians.

        Returns
        -------
        Coords
            The global coordinates of the points that were specified with
            cylindrical coordinates as input.

        Notes
        -----
        The scaling can also be applied independently prior to transforming.
        ``X.cylindrical(scale=s)`` is equivalent with
        ``X.scale(s).cylindrical()``. The scale option is provided here
        because in many cases you need  at least to scale the theta direction
        to have proper angle values.

        See Also
        --------
        hyperCylindrical: similar but allowing scaling as function of angle
        toCylindrical: inverse transformation (cartesian to cylindrical)

        Examples
        --------
        We want to create six points on a circle with radius 2. We start
        by creating the points in cylindrical coordinates with unit distances.

        >>> X = Coords('1'+'2'*5)
        >>> print(X)
        [[ 1.  0.  0.]
         [ 1.  1.  0.]
         [ 1.  2.  0.]
         [ 1.  3.  0.]
         [ 1.  4.  0.]
         [ 1.  5.  0.]]

        Remember these are (r,theta,z) coordinates of the points. So we
        will scale the r-direction with 2 (the target radius) and the
        angular direction theta with 360/6 = 60.
        Then we get the cartesian coordinates of the points from

        >>> Y = X.cylindrical(scale=(2.,60.,1.))
        >>> print(Y)
        [[ 2.    0.    0.  ]
         [ 1.    1.73  0.  ]
         [-1.    1.73  0.  ]
         [-2.   -0.    0.  ]
         [-1.   -1.73  0.  ]
         [ 1.   -1.73  0.  ]]

        Going back to cylindrical coordinates yields

        >>> print(Y.toCylindrical())
        [[   2.    0.    0.]
         [   2.   60.    0.]
         [   2.  120.    0.]
         [   2. -180.    0.]
         [   2. -120.    0.]
         [   2.  -60.    0.]]

        This differs from the original input X because of the scaling
        factors, and the wrapping around angles are reported in the
        range [-180,180].

        """
        f = zeros_like(self)
        theta = (scale[1]*angle_spec) * self[..., dir[1]]
        r = scale[0] * self[..., dir[0]]
        f[..., 0] = r*cos(theta)
        f[..., 1] = r*sin(theta)
        f[..., 2] = scale[2] * self[..., dir[2]]
        return f


    def hyperCylindrical(self,dir=(0, 1, 2),scale=(1., 1., 1.),rfunc=None,zfunc=None,angle_spec=DEG):
        """Convert cylindrical coordinates to cartesian with advanced scaling.

        This is similar to :meth:`cylindrical` but allows the specification
        of two functions defining extra scaling factors for the r and z
        directions that are dependent on the theta value.

        Parameters
        ----------
        (dir, scale, angle_spec): see :meth:`cylindrical`
        rfunc: callable, optional
            Function r(theta) taking one, float parameter and returning a float.
            Like scale[0] it is multiplied with the provided r values before
            converting them to cartesian coordinates.
        zfunc: callable, optional
            Function z(theta) taking one float parameter and returning a float.
            Like scale[2] it is multiplied with the provided z values before
            converting them to cartesian coordinates.

        See Also
        --------
        cylindrical: similar but without the rfunc and zfunc options.

        """
        if rfunc is None:
            rfunc = lambda x:1
        if zfunc is None:
            zfunc = lambda x:1
        f = zeros_like(self)
        theta = (scale[1]*angle_spec) * self[..., dir[1]]
        r = scale[0] * rfunc(theta) * self[..., dir[0]]
        f[..., 0] = r * cos(theta)
        f[..., 1] = r * sin(theta)
        f[..., 2] = scale[2] * zfunc(theta) * self[..., dir[2]]
        return f


    def toCylindrical(self,dir=(0, 1, 2),angle_spec=DEG):
        """Converts from cartesian to cylindrical coordinates.

        Returns a Coords where the values are the coordinates of the
        input points in a cylindrical coordinate system. The three axes
        of the Coords then correspond to (r, theta, z).

        Parameters
        ----------
        dir: tuple of ints
            A permutation of (0,1,2) specifying which of the global axes are
            the radial, circumferential and axial direction of the
            cylindrical coordinate system. Make sure to keep the axes ordering
            in order to get a right-handed system.
        angle_spec: float, DEG or RAD, optional
            Multiplication factor for angle coordinates.
            The default (DEG) returns angles in degrees. Use RAD to
            return angles in radians.

        Returns
        -------
        Coords
            The cylindrical coordinates of the input points.

        See Also
        --------
        cylindrical: conversion from cylindrical to cartesian coordinates

        Examples
        --------
        see :meth:`cylindrical`

        """
        f = zeros_like(self)
        x, y, z = ( self[..., i] for i in dir )
        f[..., 0] = sqrt(x*x+y*y)
        f[..., 1] = arctand2(y, x, angle_spec)
        f[..., 2] = z
        return f


    def spherical(self,dir=(0,1,2),scale=(1.,1.,1.),angle_spec=DEG,colat=False):
        """Convert spherical coordinates to cartesian coordinates.

        Consider a spherical coordinate system with the global xy-plane as
        its equatorial plane and the z-axis as axis. The zero meridional
        halfplane is taken along th positive x-axis.
        The spherical coordinates of a point are:

        - the longitude (theta): the circumferential angle, measured around
          the z-axis from the zero-meridional halfplane to the meridional
          plane containing the point: this angle normally ranges from -180 to
          +180 degrees (or from 0 to 360);
        - the latitude (phi): the elevation angle of the point's position
          vector, measured from the equatorial plane, positive when the point
          is at the positive side of the plane: this angle is normally
          restricted to the range from -90 (south pole) to +90 (north pole);
        - the distance (r): the radial distance of the point from the origin:
          this is normally positive.

        This function interpretes the 3 coordinates of the points
        as (theta,phi,r) values and computes the corresponding global
        cartesian coordinates (x,y,z).

        Parameters
        ----------
        dir: tuple of 3 ints, optional
            If provided, it is a permutation of (0,1,2) and specifies
            which of the current coordinates are interpreted as resp.
            longitude(theta), latitude(phi) and distance(r). This allows
            the axis to be aligned with any of the global axes.
            Default order is (0,1,2), with (0,1) the equatorial plane
            and 2 the axis.
            Beware that using a permutation that is not conserving the order
            of the globale axes (0,1,2), may lead to a confusing left-handed
            system.
        scale: tuple of 3 floats, optional
            Scaling factors that are applied on the coordinate values prior
            to making the conversion from spherical to cartesian coordinates.
            These factors are always given in the order (theta,phi,rz),
            irrespective of the permutation by `dir`.
        angle_spec: float, DEG or RAD, optional
            Multiplication factor for angle coordinates.
            The default (DEG) interpretes the angles in degrees. Use RAD to
            specify the angles in radians.
        colat: bool
            If True, the second coordinate is the colatitude instead. The
            colatitude is the angle measured from the north pole towards
            the south. In degrees, it is equal to ``90 - latitude`` and
            ranges from 0 to 180.
            Applications that deal with regions around the pole may
            benefit from using this option.

        Returns
        -------
        Coords
            The global coordinates of the points that were specified with
            spherical coordinates as input.

        See Also
        --------
        toSpherical: the inverse transformation (cartesian to spherical)
        cylindrical: similar function for spherical coordinates

        Examples
        --------
        >>> X = Coords('0123').scale(90).trl(2,1.)
        >>> X
        Coords([[  0.,   0.,   1.],
                [ 90.,   0.,   1.],
                [ 90.,  90.,   1.],
                [  0.,  90.,   1.]])
        >>> X.spherical()
        Coords([[ 1.,  0.,  0.],
                [-0.,  1.,  0.],
                [ 0., -0.,  1.],
                [-0., -0.,  1.]])

        Note that the last two points, though having different spherical
        coordinates, are coinciding at the north pole.
        """
        f = self.reshape((-1, 3))
        theta = (scale[0]*angle_spec) * f[:, dir[0]]
        phi = (scale[1]*angle_spec) * f[:, dir[1]]
        r = scale[2] * f[:, dir[2]]
        if colat:
            phi = 90.0*angle_spec - phi
        rc = r*cos(phi)
        f = column_stack([rc*cos(theta), rc*sin(theta), r*sin(phi)])
        return Coords(f.reshape(self.shape))


    def superSpherical(self,n=1.0,e=1.0,k=0.0, dir=(0,1,2),scale=(1.,1.,1.),angle_spec=DEG,colat=False):
        """Performs a superspherical transformation.

        superSpherical is much like :meth:`spherical`, but adds some extra
        parameters to enable the quick creation of a wide range of complicated
        shapes.
        Again, the input coordinates are interpreted as
        the longitude, latitude and distance in a spherical coordinate system.

        Parameters
        ----------
        n: float, >=0
            Exponent defining the variation of the distance in nort-south
            (latitude) direction. The default value 1 turns constant r-values
            into circular meridians. See notes.
        e: float, >=0
            Exponent defining the variation of the distance in nort-south
            (latitude) direction. The default value 1 turns constant r-values
            into a circular latitude lines. See notes.
        k: float, -1 < k < 1
            Eggness factor. If nonzero, creates asymmetric northern and
            southern hemisheres. Values > 0 enlarge the southern hemisphere
            and shrink the northern, while negative values yield the
            opposite.
        dir: tuple of 3 ints, optional
            If provided, it is a permutation of (0,1,2) and specifies
            which of the current coordinates are interpreted as resp.
            longitude(theta), latitude(phi) and distance(r). This allows
            the axis to be aligned with any of the global axes.
            Default order is (0,1,2), with (0,1) the equatorial plane
            and 2 the axis.
            Beware that using a permutation that is not conserving the order
            of the globale axes (0,1,2), may lead to a confusing left-handed
            system.
        scale: tuple of 3 floats, optional
            Scaling factors that are applied on the coordinate values prior
            to making the conversion from spherical to cartesian coordinates.
            These factors are always given in the order (theta,phi,rz),
            irrespective of the permutation by `dir`.
        angle_spec: float, DEG or RAD, optional
            Multiplication factor for angle coordinates.
            The default (DEG) interpretes the angles in degrees. Use RAD to
            specify the angles in radians.
        colat: bool
            If True, the second coordinate is the colatitude instead. The
            colatitude is the angle measured from the north pole towards
            the south. In degrees, it is equal to ``90 - latitude`` and
            ranges from 0 to 180.
            Applications that deal with regions around the pole may
            benefit from using this option.

        Raises
        ------
        ValueError
            If one of `n`, `e` or `k` is out of the acceptable range.

        Notes
        -----
        Values of `n` and `e` should not be negative. Values equal to 1
        create a circular shape. Other values keep the radius at angles
        corresponding to mmultiples of 90 degrees, while the radius at the
        intermediate 45 degree angles will be  maximally changed. Values
        larger than 1 shrink at 45 degrees directions, while lower values
        increase it. A value 2 creates a straight line between the 90 degrees
        points (the radius at 45 degrees being reduced to 1/sqrt(2).

        See also example SuperShape.

        Examples
        --------
        >>> X = Coords('02222').scale(22.5).trl(2,1.)
        >>> X
        Coords([[  0. ,   0. ,   1. ],
                [  0. ,  22.5,   1. ],
                [  0. ,  45. ,   1. ],
                [  0. ,  67.5,   1. ],
                [  0. ,  90. ,   1. ]])
        >>> X.superSpherical(n=3).toSpherical()
        Coords([[ 90.  ,   0.  ,   1.  ],
                [ 85.93,   0.  ,   0.79],
                [ 45.  ,   0.  ,   0.5 ],
                [  4.07,   0.  ,   0.79],
                [ -0.  ,  -0.  ,   1.  ]])

        The result is smaller radius at angle 45.


        """
        if n < 0. or e < 0. or k <= -1. or k >= 1.:
            raise ValueError("n, e or k out of acceptable range")

        def c(o, m):
            c = cos(o)
            return sign(c)*abs(c)**m
        def s(o, m):
            c = sin(o)
            return sign(c)*abs(c)**m

        f = self.reshape((-1, 3))
        theta = (scale[0]*angle_spec) * f[:, dir[0]]
        phi = (scale[1]*angle_spec) * f[:, dir[1]]
        r = scale[2] * f[:, dir[2]]
        if colat:
            phi = 90.0*angle_spec - phi
        rc = r*c(phi, n)
        if k != 0:   # k should be > -1.0 !!!!
            x = sin(phi)
            rc *= (1.-k*x)/(1.+k*x)
        f = column_stack([rc*c(theta, e), rc*s(theta, e), r*s(phi, n)])
        return Coords(f.reshape(self.shape))


    def toSpherical(self,dir=[0, 1, 2],angle_spec=DEG):
        """Converts from cartesian to spherical coordinates.

        Returns a Coords where the values are the coordinates of the
        input points in a spherical coordinate system. The three axes
        of the Coords then correspond to (theta, phi, r).

        Parameters
        ----------
        dir: tuple of ints
            A permutation of (0,1,2) specifying how the spherical coordinate
            system is oriented in the global axes. The last value is the axis
            of the system; the first two values are the equatorial plane;
            the first and last value define the meridional zero plane.
            Make sure to preserve the axes ordering in order to get a
            right-handed system.
        angle_spec: float, DEG or RAD, optional
            Multiplication factor for angle coordinates.
            The default (DEG) returns angles in degrees. Use RAD to
            return angles in radians.

        Returns
        -------
        Coords
            The spherical coordinates of the input points.

        See Also
        --------
        spherical: conversion from spherical to cartesian coordinates

        Examples
        --------
        See :meth:`superSpherical`

        """
        v = self[..., dir].reshape((-1, 3))
        dist = sqrt(sum(v*v,-1))
        long = arctand2(v[:,0], v[:,2], angle_spec)
        lat = where(dist <= 0.0, 0.0, arcsind(v[:,1]/dist, angle_spec))
        f = column_stack([long, lat, dist])
        return Coords(f.reshape(self.shape))


    def circulize(self,n):
        """Transform sectors of a regular polygon into circular sectors.

        Parameters
        ----------
        n: int
            Number of edges of the regular polygon.

        Returns
        -------
        Coords
            A Coords where the points inside each sector of a n-sided
            regular polygon around the origin are reposition to fill
            a circular sector. The polygon is in the x-y-plane and has a
            vertex on the x-axis.

        Notes
        -----
        Points on the x-axis and on radii at i * 360 / n degrees are not moved.
        Points on the bisector lines between these radii are move maximally
        outward. Points on a regular polygon will become points on a circle
        if circulized with parameter n equal to the number of sides of the
        polygon.

        Examples
        --------
        >>> Coords([[1.,0.],[0.5,0.5],[0.,1.]]).circulize(4)
        Coords([[ 1.  ,  0.  ,  0.  ],
               [ 0.71,  0.71,  0.  ],
               [-0.  ,  1.  ,  0.  ]])

        """
        if n < 3:
            raise ValueError("n should be at least 3")
        angle = 360./n
        X = self.toCylindrical()
        t = X.y / angle  # ranges from 0 to 1 for sector 0..angle degrees
        # Reduce values to a single sector
        while t.min() < 0.0:
            t[t<0.0] += 1.0
        while t.max() > 1.0:
            t[t>1.0] -= 1.0
        u = abs(0.5-t) * angle
        c = cosd(u) / cosd(angle/2)
        X.x *= c
        return X.cylindrical()


    def bump(self,dir,a,func=None,dist=None,xb=1.):
        """Create a 1-, 2-, or 3-dimensional bump in a Coords.

        A bump is a local modification of the coordinates of a collection
        of points. The bump can be 1-, 2- or 3-dimensional, meaning that
        the intensity of the coordinate modification varies in 1, 2 or 3
        axis directions. In all cases, the bump only changes one coordinate
        of the points.
        This method can produce various effects, but one of the most common
        uses is to force a surface to be indented at some point.

        Parameters
        ----------
        dir: int, one of (0,1,2)
            The axis of the coordinates to be modified.
        a: point (3,)
            The point that sets the bump location and intensity.
        func: callable, optional
            A function that returns the bump intensity in function
            of the distance from the bump point `a`. The distance is the
            Euclidean distance over all directions except `dir`.
            The function takes a single (positive) float value and
            returns a float (the bump intensity). Its value should
            not be zero at the origin. The function may include constants,
            which can be specified as `xb`. If no function is specified,
            the default function will be used:
            ``lambda x: where(x<xb,1.-(x/3)**2,0)``
            This makes the bump quadratically die out over a distance `xb`.
        dist: int or tuple of ints, optional
            Specifies how the distance from points to the bump point `a` is
            measured. It can be a single axis number (0,1,2) or a tuple of
            two or three axis numbers. If a single axis, the bump
            will vary only in one direction and distance is measured
            along that direction and is signed. The bump will only vary
            in that direction.
            If two or three axes, distance is the (always positive)
            euclidean distance over these directions and the bump will
            vary in these directions. Default value is the set of 3 axes
            minus the direction of modification `dir`.
        xb: float or list of floats
            Constant(s) to be used in func. Often, this includes the
            distance over which the bump will extend. The default
            bump function will reach zero at this distance.

        Returns
        -------
        Coords
            A Coords with same shape as input, but having a localized
            change of coordinates as specified by the parameters.

        Notes
        -----
        This function replaces the `bump1` and `bump2` functions in older
        pyFormex versions. The default value of `dist` makes it work like
        `bump2`. Specifyin a single axis for `dist` makes it work like `bump1`.

        See also examples BaumKuchen, Circle, Clip, Novation

        Examples
        --------
        One-dimensional bump in a linear set of points.

        >>> X = Coords(arange(6).reshape(-1,1))
        >>> X.bump1(1,[3.,5.,0.],dist=0)
        Coords([[ 0.,  0.,  0.],
                [ 1.,  0.,  0.],
                [ 2.,  0.,  0.],
                [ 3.,  5.,  0.],
                [ 4.,  0.,  0.],
                [ 5.,  0.,  0.]])
        >>> X.bump(1,[3.,5.,0.],dist=0,xb=3.)
        Coords([[ 0.  ,  0.  ,  0.  ],
                [ 1.  ,  2.78,  0.  ],
                [ 2.  ,  4.44,  0.  ],
                [ 3.  ,  5.  ,  0.  ],
                [ 4.  ,  4.44,  0.  ],
                [ 5.  ,  2.78,  0.  ]])

        Create a grid a points in xz-plane, with a bump in direction y
        with a maximum 5 at x=1.5, z=0., extending over a distance 2.5.

        >>> X = Coords(arange(4).reshape(-1,1)).replicate(4,dir=2)
        >>> X.bump(1,[1.5,5.,0.],xb=2.5)
        Coords([[[ 0.  ,  3.75,  0.  ],
                 [ 1.  ,  4.86,  0.  ],
                 [ 2.  ,  4.86,  0.  ],
                 [ 3.  ,  3.75,  0.  ]],
        <BLANKLINE>
                [[ 0.  ,  3.19,  1.  ],
                 [ 1.  ,  4.31,  1.  ],
                 [ 2.  ,  4.31,  1.  ],
                 [ 3.  ,  3.19,  1.  ]],
        <BLANKLINE>
                [[ 0.  ,  0.  ,  2.  ],
                 [ 1.  ,  2.64,  2.  ],
                 [ 2.  ,  2.64,  2.  ],
                 [ 3.  ,  0.  ,  2.  ]],
        <BLANKLINE>
                [[ 0.  ,  0.  ,  3.  ],
                 [ 1.  ,  0.  ,  3.  ],
                 [ 2.  ,  0.  ,  3.  ],
                 [ 3.  ,  0.  ,  3.  ]]])

        """
        if func is None:
            func = lambda x: where(abs(x)<xb,1.-(x/3)**2,0)
        if func(0.) == 0.:
            raise ValueError("Invalid func: f(0)=0")
        f = self.copy()
        if dist is None:
            dist = otherAxes(dir)
        if isInt(dist):
            dist = [dist]
        d = f[..., dist[0]] - a[dist[0]]
        if len(dist) > 1:
            # compute euclidean distance
            d = d*d
            for i in dist[1:]:
                d1 = f[..., i] - a[i]
                d += d1*d1
            d = sqrt(d)
        f[..., dir] += func(d)*a[dir]/func(0)
        return f


    def flare (self,xf,f,dir=(0,2),end=0,exp=1.):
        """Create a flare at the end of a :class:`Coords` block.

        A flare is a local change of geometry (widening, narrowing)
        at the end of a structure.

        Parameters
        ----------
        xf: float
            Length over which the local change occurs, measured along ``dir[0]``.
        f: float
            Maximum amplitude of the flare, in the direction ``dir[1]``.
        dir: tuple of two ints (0,1,2)
            Two axis designations. The first axis defines the direction along
            which the flare decays. The second is the direction of the
            coordinate modification.
        end: 0 or 1
            With end=0, the flare exists at the end with the smallest
            coordinates in ``dir[0]]`` direction; with end=1, at the
            end with the highest coordinates.
        exp: float
            Exponent setting the speed of decay of the flare. The default
            makes the flare change linearly over the length `f`.

        Returns
        -------
        Coords
            A Coords with same shape as the input, but having a localized
            change of coordinates at one end of the point set.

        Examples
        --------
        >>> Coords(arange(6).reshape(-1,1)).flare(3.,1.6,(0,1),0)
        Coords([[ 0.  ,  1.6 ,  0.  ],
                [ 1.  ,  1.07,  0.  ],
                [ 2.  ,  0.53,  0.  ],
                [ 3.  ,  0.  ,  0.  ],
                [ 4.  ,  0.  ,  0.  ],
                [ 5.  ,  0.  ,  0.  ]])
        """
        ix, iz = dir
        bb = self.bbox()
        if end == 0:
            xmin = bb[0][ix]
            endx = self.test(dir=ix, max=xmin+xf)
            func = lambda x: (1.-(x-xmin)/xf) ** exp
        else:
            xmax = bb[1][ix]
            endx = self.test(dir=ix, min=xmax-xf)
            func = lambda x: (1.-(xmax-x)/xf) ** exp
        x = self.copy()
        x[endx, iz] += f * func(x[endx, ix])
        return x


    def map(self, func):
        """Map a :class:`Coords` by a 3-D function.

        This allows any mathematical transformation being applied to the
        coordinates of the Coords.

        Parameters
        ----------
        func: callable
            A function taking three float arguments (x,y,z coordinates of a
            point) and returning a tuple of three float values: the new
            coordinate values to replace (x,y,z).

            The function must be applicable to NumPy arrays, so it should
            only include numerical operations and functions understood by the
            numpy module.

            Often an inline lambda function is used, but a normally defined
            function will work as well.

        Returns
        -------
        Coords object
            The input Coords mapped through the specified function

        See Also
        --------
        map1: apply a 1-dimensional mapping to one coordinate direction
        mapd: map one coordinate by a function of the distance to a point

        Notes
        -----
        See also examples Cones, Connect, HorseTorse, Manantiales, Mobius,
        ScallopDome

        Examples
        --------
        >>> print(Coords([[1.,1.,1.]]).map(lambda x,y,z: [2*x,3*y,4*z]))
        [[ 2. 3. 4.]]

        """
        # we flatten coordinate sets to ease use of complicated functions
        # we should probably do this for map1 and mapd too
        X = self.points()
        f = zeros_like(X)
        f[..., 0], f[..., 1], f[..., 2] = func(X.x, X.y, X.z)
        return f.reshape(self.shape)


    def map1(self,dir,func,x=None):
        """Map one coordinate by a 1-D function of one coordinate.

        Parameters
        ----------
        dir: int (0,1 or 2)
            The coordinate axis to be modified.
        func: callable
            Function taking a single float argument (the coordinate `x`)
            and returning a float value: the new coordinate to replace the
            `dir` coordinate.

            The function must be applicable to NumPy arrays, so it should
            only include numerical operations and functions understood by the
            numpy module.

            Often an inline lambda function is used, but a normally defined
            function will work as well.
        x: int(0,1,2), optional
            If provided, specifies the coordinate that is used as argument
            in `func`. Default is to use the same as `dir`.

        Returns
        -------
        Coords object
            The input Coords where the `dir` coordinate has been mapped
            through the specified function.

        See Also
        --------
        map: apply a general 3-dimensional mapping function
        mapd: map one coordinate by a function of the distance to a point

        Notes
        -----
        See also example SplineSurface

        Examples
        --------
        >>> Coords(arange(4).reshape(-1,1)).map1(1,lambda x:0.1*x,0)
        Coords([[ 0. ,  0. ,  0. ],
                [ 1. ,  0.1,  0. ],
                [ 2. ,  0.2,  0. ],
                [ 3. ,  0.3,  0. ]])

        """
        if x is None:
            x = dir
        f = self.copy()
        f[..., dir] = func(self[..., x])
        return f


    def mapd(self,dir,func,point=(0., 0., 0.),dist=None):
        """Map one coordinate by a function of the distance to a point.

        Parameters
        ----------
        dir: int (0, 1 or 2)
           The coordinate that will be replaced with ``func(d)``, where `d`
           is calculated as the distance to `point`.
        func: callable
            Function taking one float argument (distance to `point`) and
            returning a float: the new value for the `dist` coordinate.
            `dir` coordinate.

            The function must be applicable to NumPy arrays, so it should
            only include numerical operations and functions understood by the
            numpy module.

            Often an inline lambda function is used, but a normally defined
            function will work as well.
        point: float array-like (3,)
            The point to where the distance is computed.
        dist: int or tuple of ints (0, 1, 2)
             The coordinate directions that are used to compute the distance
             to `point`. The default is to use 3-D distances.

        Examples
        --------
        Map a regular 4x4 point grid in the xy-plane onto a sphere with
        radius 1.5 and center at the corner of the grid.

        >>> from .simple import regularGrid
        >>> X = Coords(regularGrid([0.,0.],[1.,1.],[3,3]))
        >>> X.mapd(2,lambda d:sqrt(1.5**2-d**2),X[0,0],[0,1])
        Coords([[[ 0.  ,  0.  ,  1.5 ],
                 [ 0.33,  0.  ,  1.46],
                 [ 0.67,  0.  ,  1.34],
                 [ 1.  ,  0.  ,  1.12]],
        <BLANKLINE>
                [[ 0.  ,  0.33,  1.46],
                 [ 0.33,  0.33,  1.42],
                 [ 0.67,  0.33,  1.3 ],
                 [ 1.  ,  0.33,  1.07]],
        <BLANKLINE>
                [[ 0.  ,  0.67,  1.34],
                 [ 0.33,  0.67,  1.3 ],
                 [ 0.67,  0.67,  1.17],
                 [ 1.  ,  0.67,  0.9 ]],
        <BLANKLINE>
                [[ 0.  ,  1.  ,  1.12],
                 [ 0.33,  1.  ,  1.07],
                 [ 0.67,  1.  ,  0.9 ],
                 [ 1.  ,  1.  ,  0.5 ]]])

        """
        f = self.copy()
        if dist is None:
            dist = [0, 1, 2]
        try:
            l = len(dist)
        except TypeError:
            l = 1
            dist = [dist]
        d = f[..., dist[0]] - point[dist[0]]
        if l==1:
            d = abs(d)
        else:
            d = d*d
            for i in dist[1:]:
                d1 = f[..., i] - point[i]
                d += d1*d1
            d = sqrt(d)
        f[..., dir] = func(d)
        return f


    def copyAxes(self,i,j,other=None):
        """Copy the coordinates along the axes j to the axes i.

        Parameters
        ----------
        i: int (0,1 2) or tuple of ints (0,1,2)
            One or more coordinate axes that should have replaced their
            coordinates by those along the axes `j`.
        j: int (0,1 2) or tuple of ints (0,1,2)
            One or more axes whose coordinates should be copied along the
            axes `i`. `j` should have the same type and length as `i`.
        other: Coords object, optional
            If provided, this is the source Coords for the coordinates. It
            should have the same shape as self. The default is to take the
            coords from self.

        Returns
        -------
        Coords object
            A Coords where the coordinates along axes `i` have been replaced
            by those along axes `j`.

        Examples
        --------
        >>> X = Coords([[1],[2]]).trl(2,5)
        >>> X
        Coords([[ 1.,  0.,  5.],
                [ 2.,  0.,  5.]])
        >>> X.copyAxes(1,0)
        Coords([[ 1.,  1.,  5.],
                [ 2.,  2.,  5.]])
        >>> X.copyAxes((0,1),(1,0))
        Coords([[ 0.,  1.,  5.],
                [ 0.,  2.,  5.]])
        >>> X.copyAxes((0,1,2),(1,2,0))
        Coords([[ 0.,  5.,  1.],
                [ 0.,  5.,  2.]])

        """
        if other is None:
            other = self
        f = self.copy()
        f[..., i] = other[..., j]
        return f


    def swapAxes(self, i, j):
        """Swap two coordinate axes.

        Parameters
        ----------
        i: int (0,1,2)
            First coordinate axis
        j: int (0,1,2)
            Second coordinate axis

        Returns
        -------
        Coords
            A Coords with interchanged `i` and `j` coordinates.

        Warning
        -------
        Coords.swapAxes merely changes the order of the elements
        along the last axis of the ndarray.
        This is quite different from :meth:`numpy.ndarray.swapaxes`,
        which is inherited  by the Coords class. The latter method
        interchanges the array axes of the ndarray, and will not yield
        a valid Coords object if the interchange involves the last axis.

        Notes
        -----
        This is equivalent with ``self.copyAxes((i,j),(j,i))``

        Swapping two coordinate axes has the same effect as mirroring
        against the bisector plane between the two axes.

        Examples
        --------
        >>> X = Coords(arange(6).reshape(-1,3))
        >>> X
        Coords([[ 0.,  1.,  2.],
               [ 3.,  4.,  5.]])
        >>> X.swapAxes(2,0)
        Coords([[ 2.,  1.,  0.],
               [ 5.,  4.,  3.]])
        >>> X.swapaxes(1,0)
        array([[ 0.,  3.],
               [ 1.,  4.],
               [ 2.,  5.]], dtype=float32)
        """
        order = [0,1,2]
        order[i],order[j] = j,i
        return self[...,order]


    def rollAxes(self,n=1):
        """Roll the coordinate axes over the given amount.

        Parameters
        ----------
        n: int
            Number of positions to roll the axes. With the default (1), the
            old axes (0,1,2) become the new axes (2,0,1).

        Returns
        -------
        Coords
            A Coords where the coordinate axes of the points have been rolled
            over `n` positions.

        Notes
        -----
        ``X.rollAxes(1)`` can also be obtained by
        ``X.copyAxes((0,1,2),(2,0,1))``. It is also equivalent with
        a rotation over -120 degrees around the trisectrice of the first
        quadrant.

        Examples
        --------
        >>> X = Coords('0123')
        >>> X
        Coords([[ 0.,  0.,  0.],
                [ 1.,  0.,  0.],
                [ 1.,  1.,  0.],
                [ 0.,  1.,  0.]])
        >>> X.rollAxes(1)
        Coords([[ 0.,  0.,  0.],
                [ 0.,  1.,  0.],
                [ 0.,  1.,  1.],
                [ 0.,  0.,  1.]])
        >>> X.rotate(120,axis=[1.,1.,1.])
        Coords([[ 0.,  0.,  0.],
                [-0.,  1., -0.],
                [-0.,  1.,  1.],
                [-0., -0.,  1.]])

        """
        return roll(self, int(n) % 3, axis=-1)


    def projectOnPlane(self,n=2,P=(0., 0., 0.)):
        """Project a :class:`Coords` on a plane.

        Creates a parallel projection of the Coords on a plane.

        Parameters
        ----------
        n: int (0,1,2) or float array-like (3,)
            The normal direction to the plane on which to project the Coords.
            If an int, it is a global axis.
        P: float array-like (3,)
            A point in the projection plane, by default the global origin.

        Returns
        -------
        Coords
            The points of the Coords projected on the specified plane.

        Notes
        -----
        For projection on a plane parallel to a coordinate plane,
        it is far more efficient to specify the normal by an axis
        number rather than by a three component vector.

        This method will also work if any or both of P and n have
        the same shape as self, or can be reshaped to the same shape.
        This will project each point on its individual plane.

        See also example BorderExtension

        Examples
        --------
        >>> X = Coords(arange(6).reshape(2,3))
        >>> X.projectOnPlane(0,P=[2.5,0.,0.])
        Coords([[ 2.5,  1. ,  2. ],
                [ 2.5,  4. ,  5. ]])
        >>> X.projectOnPlane([1.,1.,0.])
        Coords([[-0.5,  0.5,  2. ],
                [-0.5,  0.5,  5. ]])
        """
        x = self.reshape(-1,3).copy()
        P = Coords(P).reshape(-1, 3)
        if isInt(n):
            x[:, n] = P[:, n]
        else:
            n = normalize(Coords(n).reshape(-1, 3))
            d = dotpr(n, x-P).reshape(-1,1)
            x -= d * n
        return x.reshape(self.shape)


    def projectOnSphere(self,radius=1.,center=(0., 0., 0.)):
        """Project a :class:`Coords` on a sphere.

        Creates a central projection of a Coords on a sphere.

        Parameters
        ----------
        radius: float, optional
            The radius of the sphere, default 1.
        center: float array-like (3,), optional
            The center of the sphere. This point should not be part
            the the Coords. The default is the origin of the global axes.

        Returns
        -------
        Coords
            A Coords with the input points projected on the specified sphere.

        Notes
        -----
        This is a central projection from the center of the sphere. If you
        want a parallel projection on a spherical surface, you can use
        :meth:`map`. See the Examples there.

        Examples
        --------
        >>> X = Coords([[x,x,1.] for x in range(1,4)])
        >>> X
        Coords([[ 1.,  1.,  1.],
               [ 2.,  2.,  1.],
               [ 3.,  3.,  1.]])
        >>> X.projectOnSphere()
        Coords([[ 0.58,  0.58,  0.58],
               [ 0.67,  0.67,  0.33],
               [ 0.69,  0.69,  0.23]])

        """
        d = self.distanceFromPoint(center)
        s = radius / d
        f = self - center
        for i in range(3):
            f[..., i] *= s
        f += center
        return f


    def projectOnCylinder(self,radius=1.,dir=0,center=[0., 0., 0.]):
        """Project the Coords on a cylinder with axis parallel to a global axis.

        Given a cylinder with axis parallel to a global axis, the points
        of the Coords are projected from the axis onto the surface of the
        cylinder.
        The default cylinder has its axis along the x-axis and a unit radius.
        No points of the :class:`Coords` should belong to the axis.

        Parameters
        ----------
        radius: float, optional
            The radius of the sphere, default 1.
        dir: int (0,1,2), optional
            The global axis parallel to the cylinder's axis.
        center: float array-like (3,), optional
            A point on the axis of the cylinder. Default is the origin of
            the global axes.

        Returns
        -------
        Coords
            A Coords with the input points projected on the specified cylinder.

        Notes
        -----
        This is a projection from the axis of the cylinder. If you
        want a parallel projection on a cylindrical surface, you can use
        :meth:`map`.

        Examples
        --------
        >>> X = Coords([[x,x,1.] for x in range(1,4)])
        >>> X
        Coords([[ 1.,  1.,  1.],
               [ 2.,  2.,  1.],
               [ 3.,  3.,  1.]])
        >>> X.projectOnCylinder()
        Coords([[ 1.  ,  0.71,  0.71],
               [ 2.  ,  0.89,  0.45],
               [ 3.  ,  0.95,  0.32]])
        """
        d = self.distanceFromLine(center, unitVector(dir))
        s = radius / d
        c = resize(asarray(center), self.shape)
        c[..., dir] = self[..., dir]
        f = self - c
        for i in range(3):
            if i != dir:
                f[..., i] *= s
        f += c
        return f


    def projectOnSurface(self,S,dir=0,missing='e',return_indices=False):
        """Project a :class:`Coords` on a triangulated surface.

        The points of the Coords are projected in the specified direction `dir`
        onto the surface S. If a point has multiple projecions in the direction,
        the one nearest to the original is returned.

        Parameters
        ----------
        S: :class:`~trisurface.TriSurface`
            A triangulated surface
        dir: int (0,1,2) or float array-like (3,)
            The direction of the projection, either a global axis direction
            or specified as a vector with three components.
        missing: 'o', 'r' or 'e'
            Specifies how to treat cases where the projective line does not
            intersect the surface:

            - 'o': return the original point,
            - 'r': remove the point from the result.
              Use `return_indices` = True to find out which original
              points correspond with the projections.
            - 'e': raise an exception (default).
        return_indices: bool, optional
            If True, also returns the indices of the points that have a
            projection on the surface.

        Returns
        -------
        x: Coords
            A Coords with the projections of the input points on the surface.
            With `missing='o'`, this will have the same shape as the input,
            but some points might not actually lie on the surface.
            With `missing='r'`, the shape will be (npoints,3) and the number
            of points may be less than the input.
        ind: int array, optional
            Only returned if `return_indices` is True: an index in the input
            Coords of the points that have a projection on the surface.
            With `missing='r'`, this gives the indices of the orginal
            points corresponding with the projections.
            With `missing='o'`, this can be used to check which points
            are located on the surface.
            The index is sequential, no matter what the shape of the input
            Coords is.

        Examples
        --------
        >>> from pyformex import simple
        >>> S = simple.sphere().scale(2).trl([0.,0.,0.2])
        >>> x = pattern('0123')
        >>> print(x)
        [[ 0. 0. 0.]
         [ 1. 0. 0.]
         [ 1. 1. 0.]
         [ 0. 1. 0.]]
        >>> xp = x.projectOnSurface(S,[0.,0.,1.])
        >>> print(xp)
        [[ 0.    0.   -1.8 ]
         [ 1.    0.   -1.52]
         [ 1.    1.   -1.2 ]
         [ 0.    1.   -1.53]]
        """
        from pyformex import olist
        from pyformex.geomtools import anyPerpendicularVector
        try:
            missing = float(missing)
        except:
            if isinstance(missing, str) and len(missing) > 0:
                if missing[0] not in '+-':
                    missing = '+' + missing
                missing = missing[:2]
            else:
                missing = None

        if isinstance(dir, int):
            dir = unitVector()
        else:
            dir = asarray(dir)
        x = self.reshape(-1, 3)
        # Create planes through x in direction n
        # WE SHOULD MOVE THIS TO geomtools?
        v1 = anyPerpendicularVector(dir)
        v2 = cross(dir, v1)
        # Create set of cuts with set of planes
        cuts = [ S.intersectionWithPlane(xi, v1) for xi in x ]
        nseg = [ c.nelems() for c in cuts ]
        # remove the empty intersections
        cutid = [ i for i, n in enumerate(nseg) if n > 0 ]
        cuts = olist.select(cuts, cutid)

        # cut the cuts with second set of planes
        cuts = [ c.toFormex().intersectionWithPlane(xi, v2).coords for c, xi in zip(cuts, x[cutid]) ]
        npts = [ p.shape[0] for p in cuts ]
        okid = [ i for i, n in enumerate(npts) if n > 0 ]
        # remove the empty intersections
        cutid = olist.select(cutid, okid)
        cuts = olist.select(cuts, okid)

        # find the points closest to self
        cuts = [ p.points()[p.closestToPoint(xi)] for p, xi in zip(cuts, x[cutid]) ]
        cuts = Coords.concatenate(cuts)

        if cuts.shape[0] < x.shape[0]:
            if missing == 'e':
                raise ValueError("The projection of some point(s) in the specified direction does not cut the surface")
            elif missing == 'o':
                x = x.copy()
                x[cutid] = cuts
                cuts = x.reshape(self.shape)

        if return_indices:
            return cuts, cutid
        else:
            return cuts


    # Extra transformations implemented by plugins

    def isopar(self, eltype, coords, oldcoords):
        """Perform an isoparametric transformation on a Coords.

        This creates an isoparametric transformation :class:`plugins.isopar.Isopar` object
        and uses it to transform the input Coords. It is equivalent to::

          Isopar(eltype,coords,oldcoords).transform(self)

        See :class:`~plugins.isopar.Isopar` for parameters.
        """
        from pyformex.plugins.isopar import Isopar
        return Isopar(eltype, coords, oldcoords).transform(self)


    def addNoise(self,rsize=0.05,asize=0.0):
        """Add random noise to a Coords.

        A random amount is added to each individual coordinate of the Coords.
        The maximum difference of the coordinates from their original value
        is controled by two parameters `rsize` and `asize` and will not
        exceed ``asize+rsize*self.maxsize()``.

        Parameters
        ----------
        rsize: float
            Relative size of the noise compared with the maximum size
            of the input Coords.
        asize: float
            Absolute size of the noise

        Examples
        --------
        >>> X = Coords(arange(6).reshape(2,3))
        >>> print((abs(X.addNoise(0.1) - X) < 0.1 * X.sizes()).all())
        True
        """
        max = asize + rsize * self.maxsize()
        return self + randomNoise(self.shape, -max, +max)


############################################################################
    #
    #   Transformations that change the shape of the Coords array
    #


    def replicate(self,n,dir=0,step=1.):
        """Replicate a Coords n times with a fixed translation step.

        Parameters
        ----------
        n: int
            Number of times to replicate the Coords.
        dir: int (0,1,2) or float array-like (3,)
            The translation vector. If an int, it specifies a global axis
            and the translation is in the direction of that axis.
        step: float
            If ``dir`` is an int, this is the length of the translation.
            Else, it is a multiplying factor applied to the translation
            vector.

        Returns
        -------
        Coords
            A Coords with an extra first axis with length `n`. The new
            shape thus becomes ``(n,) + self.shape``.
            The first component along the axis 0 is identical to the
            original Coords. Each following component is equal to the
            previous translated over `(dir,step)`, where `dir` and `step` are
            interpreted just like in the :meth:`translate` method.

        Notes
        -----
        :meth:`rep` is a convenient shorthand for :meth:`replicate`.

        Examples
        --------
        >>> Coords([0.,0.,0.]).replicate(4,1,1.2)
        Coords([[ 0. ,  0. ,  0. ],
                [ 0. ,  1.2,  0. ],
                [ 0. ,  2.4,  0. ],
                [ 0. ,  3.6,  0. ]])
        >>> Coords([0.]).replicate(3,0).replicate(2,1)
        Coords([[[ 0.,  0.,  0.],
                 [ 1.,  0.,  0.],
                 [ 2.,  0.,  0.]],
        <BLANKLINE>
                [[ 0.,  1.,  0.],
                 [ 1.,  1.,  0.],
                 [ 2.,  1.,  0.]]])

        """
        n = int(n)
        f = np.resize(self, (n,)+self.shape)
        if isInt(dir):
            for i in range(1, n):
                f[i,...,dir] += i*step
        else:
            dir = Coords(dir, copy=True)
            if step != 1.:
                dir *= step
            for i in range(1, n):
                f[i] += i*dir
        return Coords(f)


    def split(self):
        """Split the Coords in blocks along first axis.

        Returns
        -------
        list of Coords objects
            A list of Coords objects being the subarrays takeb along the
            axis 0. The number of objects in the list is ``self.shape[0]``
            and each Coords has the shape ``self.shape[1:]``.

        Raises
        ------
        ValueError
            If ``self.ndim < 2``.

        Examples
        --------
        >>> Coords(arange(6).reshape(2,3)).split()
        [Coords([ 0.,  1.,  2.]), Coords([ 3.,  4.,  5.])]
        """
        if self.ndim < 2:
            raise ValueError("Can only split arrays with dim >= 2")
        return [ self[i] for i in range(self.shape[0]) ]


    def sort(self,order=(0, 1, 2)):
        """Sort points in the specified order of their coordinates.

        Parameters
        ----------
        order: int (0,1,2) or tuple of ints (0,1,2)
            The order in which the coordinates have to be taken into
            account during the sorting operation.

        Returns
        -------
        int array
            An index into the sequential point list ``self.points()``
            thus that the points are sorted in order of the specified
            coordinates.

        Examples
        --------
        >>> X = Coords([[5,3,0],[2,4,3],[2,3,3],[5,6,2]])
        >>> X.sort()
        array([2, 1, 0, 3])
        >>> X.sort((2,1,0))
        array([0, 3, 2, 1])
        >>> X.sort(1)
        array([0, 2, 1, 3])
        """
        if isInt(order):
            order = (order,)
        return sortByColumns(self.points()[:,order])


    def boxes(self,ppb=1,shift=0.5,minsize=1.e-5):
        """Create a grid of equally sized boxes spanning the :class:`Coords`.

        A regular 3D grid of equally sized boxes is created enclosing all
        the points of the Coords. The size, position and number of boxes
        are determined from the specified parameters.

        Parameters
        ----------
        ppb: int
            Average number of points per box. The box sizes and number
            of boxes will be determined to approximate this number.
        shift: float (0.0 .. 1.0)
            Relative shift value for the grid. Applying a shift of 0.5
            will make the lowest coordinate values fall at the center
            of the outer boxes.
        minsize: float
            Absolute minimal size of the boxes, in each coordinate direction.

        Returns
        -------
        ox: float array (3,)
            The minimal coordinates of the box grid.
        dx: float array (3,)
            The box size in the three global axis directions.
        nx: int array (3,)
            Number of boxes in each of the coordinate directions.

        Notes
        -----
        The primary purpose of this method is its use in the :meth:`fuse`
        method. The boxes allow to quickly label the points inside each
        box with an integer value (the box number), so that it becomes
        easy to find close points by their same label.

        Because of the possibility that two very close points fall in
        different boxes (if they happen to be close to a box border),
        procedures based on these boxes are often repeated twice, with
        a different shift value.

        Examples
        --------
        >>> X = Coords([[5,3,0],[2,4,3],[2,3,3],[5,6,2]])
        >>> print(*X.boxes())
        [ 0.5  1.5 -1.5] [ 3.  3.  3.] [2 2 2]
        >>> print(* X.boxes(shift=0.1))
        [ 1.7  2.7 -0.3] [ 3.  3.  3.] [2 2 2]
        >>> X = Coords([[1.,1.,0.],[1.001,1.,0.],[1.1,1.,0.]])
        >>> print(*X.boxes())
        [ 0.98  0.98 -0.02] [ 0.03  0.03  0.03] [4 1 1]
        """
        # serialize points
        x = self.reshape(-1, 3)
        nnod = x.shape[0]
        # Calculate box size
        lo, hi = x.bbox()
        sz = hi-lo
        esz = sz[sz > 0.0]  # only keep the nonzero dimensions

        if esz.size == 0:
            # All points are coincident
            ox = zeros(3, dtype=Float)
            dx = ones(3, dtype=Float)
            nx = ones(3, dtype=Int)

        else:
            nboxes = max(1, nnod // ppb) # ideal total number of boxes
            vol = esz.prod()
            # avoid error message on the global sz/nx calculation
            errh = seterr(all='ignore')
            # set ideal box size, but not smaller than minsize
            boxsz = max(minsize, (vol/nboxes) ** (1./esz.shape[0]))
            nx = (sz/boxsz).astype(int32)
            dx = where(nx>0, sz/nx, boxsz)
            seterr(**errh)

        # perform the origin shift and adjust nx to make sure we enclose all
        ox = lo - dx*shift
        ex = ox + dx*nx
        adj = ceil((hi-ex)/dx).astype(Int).clip(min=0)
        nx += adj
        return ox, dx, nx


    def fuse(self,ppb=1,shift=0.5,rtol=1.e-5,atol=1.e-8,repeat=True):
        """Find (almost) coinciding points and return a compressed set.

        This method finds the points that are very close to each other
        and replaces them with a single point. See Notes below for
        explanation about the method being used and the parameters being
        used. In most cases, `atol` and `rtol` are probably the only
        ones you want to change from the defaults. Two points are
        considered the same if all their coordinates differ less than
        the maximum of `atol` and `rtol * self.maxsize()`.

        Parameters
        ----------
        ppb: int, optional
            Average number of points per box. The box sizes and number
            of boxes will be determined to approximate this number.
        shift: float (0.0 .. 1.0), optional
            Relative shift value for the box grid. Applying a shift of 0.5
            will make the lowest coordinate values fall at the center
            of the outer boxes.
        rtol: float, optional
            Relative tolerance used when considering two points for fusing.
        atol: float, optional
            Absolute tolerance used when considering two points for fusing.
        repeat: bool, optional
            If True, repeat the procedure with a second shift value.

        Returns
        -------
        coords: Coords object (npts,3)
            The unique points obtained from merging the very close points
            of a Coords.
        index: int array
            An index in the unique coordinates array `coords` for each of
            the original points. The shape of the index array is equal to
            the point shape of the input Coords (``self.pshape()``). All
            the values are in the range 0..npts.


        Note
        ----
        From the return values ``coords[index]`` will restore the original
        Coords (with accuracy equal to the tolerance used in the fuse
        operation)

        Notes
        -----
        The procedure works by first dividing the 3D space in a number of
        equally sized boxes, with a average population of `ppb` points. The
        arguments `pbb` and `shift` are passed to :meth:`boxes` for this
        purpose.
        The boxes are identified by 3 integer coordinates, from which a
        unique integer scalar is computed, which is then used to sort the
        points.
        Finally only the points inside the same box need to be compared.
        Two points are considered equal if all their coordinates differ less
        than the maximum of `atol` and `rtol * self.maxsize()`.
        Points considered very close are replaced by a single one, and an index
        is kept from the original points to the new list of points.

        Running the procedure once does not guarantee finding all close nodes:
        two close nodes might be in adjacent boxes. The performance hit for
        testing adjacent boxes is rather high, and the probability of separating
        two close nodes with the computed box limits is very small.
        Therefore, the most sensible way is to run the procedure twice, with
        a different shift value (they should differ more than the tolerance).
        Specifying repeat=True will automatically do this with a second
        shift value equal to shift+0.25.

        Because fusing points is a very important and frequent step in many
        geometrical modeling and conversion procedures, the core part of
        this function is available in a C as well as a Python version,
        in the module pyformex.lib.misc. The much faster C version will be used
        if available.

        Examples
        --------
        >>> X = Coords([[1.,1.,0.],[1.001,1.,0.],[1.1,1.,0.]])
        >>> x,e = X.fuse(atol=0.01)
        >>> print(x)
        [[ 1.  1.  0. ]
         [ 1.1 1.  0. ]]
        >>> print(e)
        [0 0 1]
        >>> allclose(X,x[e],atol=0.01)
        True

        """
        from pyformex.lib import misc
        if self.size == 0:
            # allow empty coords sets
            return Coords(), array([], dtype=Int).reshape(self.pshape())

        if repeat:
            # Apply twice with different shift value
            coords, index = self.fuse(ppb, shift, rtol, atol, repeat=False)
            newshift = shift + (0.25 if shift <= 0.5 else -0.25)
            coords, index2 = coords.fuse(ppb, newshift, rtol, atol, repeat=False)
            index = index2[index]
            return coords, index

        #########################################
        # This is the single pass
        x = self.points()

        if (self.sizes()==0.).all():
            # All points are coincident
            e = zeros(x.shape[0], dtype=int32)
            x = x[:1]
            return x, e

        # Compute boxes
        ox, dx, nx = self.boxes(ppb=ppb, shift=shift, minsize=atol)

        # Create box coordinates for all nodes
        ind = floor((x-ox)/dx).astype(int32)
        # Create unique box numbers in smallest direction first
        o = argsort(nx)
        val = ( ind[:, o[2]] * nx[o[2]] + ind[:, o[1]] ) * nx[o[1]] + ind[:, o[0]]
        # sort according to box number
        srt = argsort(val)
        # rearrange the data according to the sort order
        val = val[srt]
        x = x[srt]
        # make sure we use int32 (for the fast C fuse function)
        # Using int32 limits this procedure to 10**9 points, which is more
        # than enough for all practical purposes
        x = x.astype(float32)
        val = val.astype(int32)
        tol = float32(max(rtol*self.maxsize(), atol))
        nnod = val.shape[0]
        flag = ones((nnod,), dtype=int32)   # 1 = new, 0 = existing node
        # new fusing algorithm
        sel = arange(nnod).astype(int32)      # replacement unique node nr
        misc._fuse(x, val, flag, sel, tol)     # fuse the close points
        x = x[flag>0]          # extract unique nodes
        s = sel[argsort(srt)]  # and indices for old nodes
        return (x, s.reshape(self.shape[:-1]))


    def unique(self,**kargs):
        """Returns the unique points after fusing.

        This is just like :meth:`fuse` and takes the same arguments,
        but only returns the first argument: the unique points in the Coords.
        """
        return self.fuse(**kargs)[0]


    def adjust(self,**kargs):
        """Find (almost) identical nodes and adjust them to be identical.

        This is like the :meth:`fuse` operation, but it does not fuse the
        close neigbours to a single point.
        Instead it adjust the coordinates of the points to be identical.

        The parameters are the same as for the :meth:`fuse` method.

        Returns
        -------
        Coords
            A Coords with the same shape as the input, but where close
            points now have identical coordinates.

        Examples
        --------
        >>> X = Coords([[1.,1.,0.],[1.001,1.,0.],[1.1,1.,0.]])
        >>> print(X.adjust(atol=0.01))
        [[ 1.  1.  0. ]
         [ 1.  1.  0. ]
         [ 1.1 1.  0. ]]

        """
        coords,index = self.fuse(**kargs)
        return coords[index]


    def match(self,coords,**kargs):
        """Match points in another :class:`Coords` object.

        Find the points from another Coords object that coincide with
        (or are very close to) points of ``self``.
        This method works by concatenating the serialized point sets of
        both Coords and then fusing them.

        Parameters
        ----------
        coords: Coords
            The Coords object to compare the points with.
        **kargs: keyword arguments
            Keyword arguments passed to the :meth:`fuse` method.

        Returns
        -------
        : 1-dim int array
            The array has a length of `coords.npoints()`. For each point
            in `coords` it holds the index of a point in `self` coinciding
            with it, or a value -1 if there is no matching point.
            If there are multiple matching points in `self`, it is undefined
            which one will be returned. To avoid this ambiguity, you can
            first fuse the points of `self`.

        See Also
        --------
        hasMatch
        fuse

        Examples
        --------
        >>> X = Coords([[1.],[2.],[3.],[1.]])
        >>> Y = Coords([[1.],[4.],[2.00001]])
        >>> print(X.match(Y))
        [ 0 -1  1]

        """
        if 'clean' in kargs:
            utils.warn('warn_coords_match_changed')
            del kargs['clean']

        x = Coords.concatenate([self.points(), coords.points()])
        c, e = x.fuse(**kargs)
        e0, e1 = e[:self.npoints()], e[self.npoints():]
        return findFirst(e0, e1)


    # TODO: This is incorrect: doubles in self are not reported
    # It is documented, but should we fix it or remove this method??
    def hasMatch(self,coords,**kargs):
        """Find out which points are also in another Coords object.

        Find the points from self that coincide with (or are very close to)
        some point of `coords`.
        This method is very similar to :meth:`match`, but does not give
        information about which point of `self` matches which point of
        `coords`.

        Parameters
        ----------
        coords: Coords
            The Coords object to compare the points with.
        **kargs: keyword arguments
            Keyword arguments passed to the :meth:`fuse` method.

        Returns
        -------
        int array
            A 1-dim int array with the unique sorted indices of the points
            in `self` that have a (nearly) matching point in `coords`.

        Warning
        -------
        If multiple points in `self` coincide with the same point in
        `coords`, only one index will be returned for this case. To avoid
        this, you can fuse `self` before using this method.

        See also
        --------
        match

        Examples
        --------
        >>> X = Coords([[1.],[2.],[3.],[1.]])
        >>> Y = Coords([[1.],[4.],[2.00001]])
        >>> print(X.hasMatch(Y))
        [0 1]

        """
        matches = self.match(coords,**kargs)
        return unique(matches[matches>-1])


    def append(self, coords):
        """Append more coords to a Coords object.

        The appended coords should have matching dimensions in all
        but the first axis.

        Parameters
        ----------
        coords: Coords object
            A Coords having a shape with ``shape[1:]`` equal to
            ``self.shape[1:]``.

        Returns
        -------
        Coords
            The concatenated Coords object (self,coords).

        Notes
        -----
        This is comparable to :func:`numpy.append`, but the result
        is a :class:`Coords` object, the default axis is the first one
        instead of the last, and it is a method rather than a function.

        See Also
        --------
        concatenate: concatenate a list of Coords

        Examples
        --------
        >>> X = Coords([[1],[2]])
        >>> Y = Coords([[3],[4]])
        >>> X.append(Y)
        Coords([[ 1.,  0.,  0.],
               [ 2.,  0.,  0.],
               [ 3.,  0.,  0.],
               [ 4.,  0.,  0.]])
        """
        return Coords(append(self, coords, axis=0))


    @classmethod
    def concatenate(clas,L,axis=0):
        """Concatenate a list of :class:`Coords` objects.

        Class method to concatenate a list of Coords along the given axis.

        Parameters
        ----------
        L: list of Coords objects
            The Coords objects to be concatenated. All should have the same
            shape except for the length of the specified axis.

        Returns
        -------
        Coords
            A Coords with at least two dimensions, even when the list contains
            only a single Coords with a single point, or is empty.

        Raises
        ------
        ValueError
            If the shape of the Coords in the list do not match or if
            concatenation along the last axis is attempted.

        Notes
        -----
        This is a class method. It is commonly invoked as
        ``Coords.concatenate``, and if used as a method on a Coords object,
        that object will not be included in the list.

        It is like :func:`numpy.concatenate` (which it uses internally),
        but makes sure to return :class:`Coords` object, and sets the first
        axis as default instead of the last (which would not make sense).

        See Also
        --------
        append: append a Coords to self

        Examples
        --------
        >>> X = Coords([1.,1.,0.])
        >>> Y = Coords([[2.,2.,0.],[3.,3.,0.]])
        >>> print(Coords.concatenate([X,Y]))
        [[ 1. 1. 0.]
         [ 2. 2. 0.]
         [ 3. 3. 0.]]
        >>> print(Coords.concatenate([X,X]))
        [[ 1. 1. 0.]
         [ 1. 1. 0.]]
        >>> print(Coords.concatenate([X]))
        [[ 1. 1. 0.]]
        >>> print(Coords.concatenate([Y]))
        [[ 2. 2. 0.]
         [ 3. 3. 0.]]
        >>> print(X.concatenate([Y]))
        [[ 2. 2. 0.]
         [ 3. 3. 0.]]
        >>> Coords.concatenate([])
        Coords([], shape=(0, 3))
        >>> Coords.concatenate([[Y],[Y]],axis=1)
        Coords([[[ 2.,  2.,  0.],
                [ 3.,  3.,  0.],
                [ 2.,  2.,  0.],
                [ 3.,  3.,  0.]]])
        """
        L2 = atleast_2d(*L)
        if len(L2) == 0 or max([len(x) for x in L2]) == 0:
            return Coords()
        if len(L)==1:
            return L2
        else:
            return Coords(data=concatenate(L2, axis=axis))


    @classmethod
    def fromstring(clas,s,sep=' ',ndim=3,count=-1):
        """Create a :class:`Coords` object with data from a string.

        This uses :func:`numpy.fromstring` to read coordinates from a
        string and creates a Coords object from them.

        Parameters
        ----------
        s: str
            A string containing a single sequence of float numbers separated
            by whitespace and a possible separator string.
        sep: str
            The separator used between the coordinates. If not a space,
            all extra whitespace is ignored.
        ndim: int,
            Number of coordinates per point. Should be 1, 2 or 3 (default).
            If 1, resp. 2, the coordinate string only holds x, resp. x,y
            values.
        count: int, optional
            Total number of coordinates to read. This should be a multiple
            of `ndim`. The default is to read all the coordinates in the
            string.

        Returns
        -------
        Coords
            A Coords object with the coordinates read from the string.

        Raises
        ------
        ValueError
            If count was provided and the string does not contain that exact
            number of coordinates.

        Notes
        -----
        For writing the coordinates to a string, :func:`numpy.tostring` can
        be used.

        Examples
        --------
        >>> Coords.fromstring('4 0 0 3 1 2 6 5 7')
        Coords([[ 4.,  0.,  0.],
               [ 3.,  1.,  2.],
               [ 6.,  5.,  7.]])
        >>> Coords.fromstring('1 2 3 4 5 6',ndim=2)
        Coords([[ 1.,  2.,  0.],
               [ 3.,  4.,  0.],
               [ 5.,  6.,  0.]])

        """
        x = fromstring(s, dtype=Float, sep=sep, count=count)
        if count > 0 and x.size != count :
            raise ValueError("Number of coordinates read: %s, expected %s!" % (x.size, count))
        if x.size % ndim != 0 :
            raise ValueError("Number of coordinates read: %s, expected a multiple of %s!" % (x.size, ndim))
        return Coords(x.reshape(-1, ndim))


    @classmethod
    def fromfile(clas,fil,**kargs):
        """Read a :class:`Coords` from file.

        This uses :func:`numpy.fromfile` to read coordinates from a file
        and create a Coords. Coordinates X, Y and Z for subsequent points
        are read from the file. The total number of coordinates on the file
        should be a multiple of 3.

        Parameters
        ----------
        fil: str or file
            If str, it is a file name. An open file object can also be passed
        **kargs:
            Arguments to be passed to :func:`numpy.fromfile`.

        Returns
        -------
        Coords
            A Coords formed by reading all coordinates from the specified file.

        Raises
        ------
        ValueError
            If the number of coordinates read is not a multiple of 3.

        See Also
        --------
        numpy.fromfile: read an array to file
        numpy.tofile: write an array to file

        """
        x = fromfile(fil,dtype=Float,**kargs)
        if x.size % 3 != 0 :
            raise ValueError("Number of coordinates read: %s, should be multiple of 3!" % x.size)
        return Coords(x.reshape(-1, 3))


    def interpolate(self, X, div):
        """Create linear interpolations between two Coords.

        A linear interpolation of two equally shaped Coords X and Y at
        parameter value t is a Coords with the same shape as X and Y
        and with coordinates given by ``X * (1.0-t) + Y * t``.

        Parameters
        ----------
        X: Coords object
            A Coords object with same shape as `self`.
        div: int or list of floats
            The list of parameter values for which to compute the
            interpolation. Usually, they are in the range 0.0 (self)
            to 1.0 (X). Values outside the range can be used however
            and result in linear extrapolations.

            If an int is provided, a list with ``(div+1)`` parameter
            values is used, obtained by dividing the interval [0..1] into
            `div` equal segments. This ``F.interpolate(G,n)`` is equivalent with
            ``F.interpolate(G,arange(0,n+1)/float(n))``.

        Returns
        -------
        Coords
            A Coords object with an extra (first) axis, containing the
            concatenation of the interpolations of `self` and `X` at all
            parameter values in `div`.
            Its shape is (n,) + self.shape, where n is the number of values
            in `div`.

        Examples
        --------
        >>> X = Coords([0])
        >>> Y = Coords([1])
        >>> X.interpolate(Y,4)
        Coords([[ 0.  ,  0.  ,  0.  ],
               [ 0.25,  0.  ,  0.  ],
               [ 0.5 ,  0.  ,  0.  ],
               [ 0.75,  0.  ,  0.  ],
               [ 1.  ,  0.  ,  0.  ]])
        >>> X.interpolate(Y,[-0.1, 0.5, 1.25])
        Coords([[-0.1 ,  0.  ,  0.  ],
                [ 0.5 ,  0.  ,  0.  ],
                [ 1.25,  0.  ,  0.  ]])

        """
        if self.shape != X.shape:
            raise RuntimeError("`X` should have same shape as `self`")
        div = unitDivisor(div)
        return self + outer(div, X-self).reshape((-1,)+self.shape)


    def convexHull(self,dir=None,return_mesh=False):
        """Return the 2D or 3D convex hull of a :class:`Coords`.

        Parameters
        ----------
        dir: int (0,1,2), optional
            If provided, it is one if the global axes and the 2D convex
            hull in the specified viewing direction will be computed.
            The default is to compute the 3D convex hull.
        return_mesh: bool, optional
            If True, returns the convex hull as a :class:`~mesh.Mesh` object
            instead of a :class:`~connectivity.Connectivity`.

        Returns
        -------
        :class:`~connectivity.Connectivity` or :class:`~mesh.Mesh`
            The default is to return a Connectivity table containing the
            indices of the points that constitute the convex hull of the
            Coords. For a 3D hull, the Connectivity has plexitude 3, and
            eltype 'tri3'; for a 2D hull these are respectively 2 and 'line2'.
            The values in the Connectivity refer to the flat points list
            as obtained from :meth:`points`.

            If `return_mesh` is True, a compacted Mesh is returned
            instead of the Connectivity. For a 3D hull, the Mesh will be a
            :class:`~trisurface.TriSurface`, otherwise
            it is a Mesh of 'line2' elements.

            The returned Connectivity or Mesh will be empty if all the points
            are in a plane for the 3D version, or an a line in the viewing
            direction for the 2D version.

        Notes
        -----
        This uses SciPy to compute the convex hull. You need to have
        SciPy version 0.12.0 or higher.

        See also example ConvexHull.

        """
        from pyformex.plugins import scipy_itf
        points = self.points()
        if dir is not None and isInt(dir):
            ind = list(range(3))
            ind.remove(dir)
            points = points[:,ind]
        hull = scipy_itf.convexHull(points)
        if return_mesh:
            from pyformex.mesh import Mesh
            hull = Mesh(self.points(),hull).compact()
            if dir is not None and isInt(dir):
                hull.coords[:,dir] = 0.0

        return hull


    def actor(self,**kargs):
        """This allows a Coords object to be drawn as Geometry"""

        if self.npoints() == 0:
            return None

        from pyformex.formex import Formex
        return Formex(self.reshape(-1, 3)).actor(**kargs)


    # Convenient shorter notations
    rot = rotate
    trl = translate
    rep = replicate


    # deprecated but kept for compatibility

    def bump1(self,dir,a,func=None,dist=0):
        return self.bump(dir,a,func,dist=dist)


###########################################################################
##
##   functions
##
#########################
#

def otherAxes(i):
    """Return all global axes except the specified one

    Parameters
    ----------
    i: int (0,1,2)
        One of the global axes.

    Returns
    -------
    tuple of ints
        Two ints (j,k) identifying the other global axes in such order
        that (i,j,k) is a right-handed coordinate system.
    """
    if i == 0:
        return (1,2)
    elif i == 1:
        return (2,0)
    elif i == 2:
        return (0,1)
    else:
        raise ValueError("Invalid axis number")


def bbox(objects):
    """Compute the bounding box of a list of objects.

    The bounding box of an object is the smallest rectangular cuboid
    in the global Cartesian coordinates, such that no points of the
    objects lie outside that cuboid. The resulting bounding box of the list
    of objects is the smallest bounding box that encloses all the objects
    in the list.

    Parameters
    ----------
    objects: object or list of objects
        One or more (list or tuple) objects that have a method :meth:`bbox`
        returning the object's bounding box as a Coords with two points.

    Returns
    -------
    Coords
        A Coords object with two points: the first contains the minimal
        coordinate values, the second has the maximal ones of the overall
        bounding box encompassing all objects.

    Notes
    -----
    Objects that do not have a :meth:`bbox` method or whose :meth:`bbox`
    method returns invalid values, are silently ignored.

    See Also
    --------
    Coords.bbox: compute the bounding box of a :class:`Coords` object.

    Examples
    --------
    >>> bbox((Coords([-1.,1.,0.]),Coords([2,-3])))
    Coords([[-1., -3.,  0.],
            [ 2.,  1.,  0.]])

    """
    if not isinstance(objects,(list,tuple)):
        objects = [ objects ]
    bboxes = [f.bbox() for f in objects if hasattr(f, 'bbox') and not isnan(f.bbox()).any()]
    bboxes = [bb for bb in bboxes if bb is not None]
    if len(bboxes) == 0:
        o = origin()
        bboxes = [ [o, o] ]
    return Coords(concatenate(bboxes)).bbox()


def bboxIntersection(A, B):
    """Compute the intersection of the bounding box of two objects.

    Parameters
    ----------
    A: first object
        An object having a bbox method returning its boundary box.
    B: second object
        Another object having a bbox method returning its boundary box.

    Returns
    -------
    Coords (2,3)
        A Coords specifying the intersection of the bounding boxes of
        the two objects. This again has the format of a bounding box:
        a coords with two points: one with the minimal and one with
        the maximal coordinates. If the two bounding boxes do not
        intersect, an empty Coords is returned.

    Notes
    -----
    Since bounding boxes are Coords objects, it is possible to pass
    computed bounding boxes as arguments. The bounding boxes are indeed
    their own bounding box.

    Examples
    --------
    >>> A = Coords([[-1.,1.],[2,-3]])
    >>> B = Coords([[0.,1.],[4,2]])
    >>> C = Coords([[0.,2.],[4,2]])
    >>> bbox((A,B))
    Coords([[-1., -3.,  0.],
            [ 4.,  2.,  0.]])

    The intersection of the bounding boxes of A and B
    degenerates into a line segment parallel to the x-axis:

    >>> bboxIntersection(A,B)
    Coords([[ 0.,  1.,  0.],
            [ 2.,  1.,  0.]])

    The bounding boxes of A and C do not intersect:

    >>> bboxIntersection(A,C)
    Coords([], shape=(0, 3))

    """
    Amin, Amax = A.bbox()
    Bmin, Bmax = B.bbox()
    min = where(Amin>Bmin, Amin, Bmin)
    max = where(Amax<Bmax, Amax, Bmax)
    if (min > max).any():
        bb = Coords()
    else:
        bb = Coords([min, max])
    return bb


def origin():
    """ReturnCreate a Coords holding the origin of the global coordinate system.

    Returns
    -------
    Coords (3,)
        A Coords holding a single point with coordinates (0.,0.,0.).

    Exmaples
    --------
    >>> origin()
    Coords([ 0.,  0.,  0.])

    """
    return Coords(zeros((3)))


def pattern(s,aslist=False):
    """Generate a sequence of points on a regular grid.

    This function creates a sequence of points that are on a regular grid
    with unit step. These points are created from a simple string input,
    interpreting each character as a code specifying how to move from the
    last to the next point.
    The start position is always the origin (0.,0.,0.).

    Currently the following codes are defined:

    - 0 or +: goto origin (0.,0.,0.)
    - 1..8: move in the x,y plane
    - 9 or .: remain at the same place (i.e. duplicate the last point)
    - A..I: same as 1..9 plus step +1. in z-direction
    - a..i: same as 1..9 plus step -1. in z-direction
    - /: do not insert the next point

    Any other character raises an error.

    When looking at the x,y-plane with the x-axis to the right and the
    y-axis up, we have the following basic moves:
    1 = East, 2 = North, 3 = West, 4 = South, 5 = NE, 6 = NW, 7 = SW, 8 = SE.

    Adding 16 to the ordinal of the character causes an extra move of +1. in
    the z-direction. Adding 48 causes an extra move of -1. This means that
    'ABCDEFGHI', resp. 'abcdefghi', correspond with '123456789' with an extra
    z +/-= 1. This gives the following schema::

                 z+=1             z unchanged            z -= 1

             F    B    E          6    2    5         f    b    e
                  |                    |                   |
                  |                    |                   |
             C----I----A          3----9----1         c----i----a
                  |                    |                   |
                  |                    |                   |
             G    D    H          7    4    8         g    d    h

    The special character '/' can be put before any character to make the
    move without inserting the new point. The string should start
    with a '0' or '9' to include the starting point (the origin) in the output.

    Parameters
    ----------
    s: str
        A string with characters generating subsequent points.
    aslist: bool, optional
        If True, the points are returned as lists of **integer**
        coordinates instead of a :class:`Coords` object.

    Returns
    -------
    Coords or list of ints
        The default is to return the generated points as a Coords.
        With ``aslist=True`` however, the points are returned as a list
        of tuples holding 3 integer grid coordinates.

    See Also
    --------
    xpattern

    Examples
    --------
    >>> pattern('0123')
    Coords([[ 0.,  0.,  0.],
            [ 1.,  0.,  0.],
            [ 1.,  1.,  0.],
            [ 0.,  1.,  0.]])
    >>> pattern('2'*4)
    Coords([[ 0.,  1.,  0.],
            [ 0.,  2.,  0.],
            [ 0.,  3.,  0.],
            [ 0.,  4.,  0.]])
    """
    x = y = z = 0
    l = []
    insert = True
    for c in s:
        if c == '/':
            insert = False
            continue
        elif c == '0' or c == '+':
            x = y = z = 0
        elif c == '.':
            pass
        else:
            j, i = divmod(ord(c), 16)
            if j == 3:
                pass
            elif j == 4:
                z += 1
            elif j == 6:
                z -= 1
            else:
                raise RuntimeError("Unknown character '%c' in pattern input" % c)
            if i == 1:
                x += 1
            elif i == 2:
                y += 1
            elif i == 3:
                x -= 1
            elif i == 4:
                y -= 1
            elif i == 5:
                x += 1
                y += 1
            elif i == 6:
                x -= 1
                y += 1
            elif i == 7:
                x -= 1
                y -= 1
            elif i == 8:
                x += 1
                y -= 1
            elif i == 9:
                pass
            else:
                raise RuntimeError("Unknown character '%c' in pattern input" % c)
        if insert:
            l.append((x, y, z))
        insert = True
    if not aslist:
        l = Coords(l)
    return l


def xpattern(s,nplex=1):
    """Create a Coords object from a string pattern.

    Create a sequence of points using :func:`pattern`, and groups the
    points by ``nplex`` to create a Coords with shape ``(-1,nplex,3)``.

    Parameters
    ----------
    s: str
        The string to pass to :func:`pattern` to produce the sequence
        of points.
    nplex: int
        The number of subsequent points to group together to create
        the structured Coords.

    Returns
    -------
    Coords
        A Coords with shape (-1,nplex,3).

    Raises
    ------
    ValueError
        If the number of points produced by the input string `s` is not
        a multiple of `nplex`.

    Examples
    --------
    >>> print(xpattern('.12.34',3))
    [[[ 0. 0. 0.]
      [ 1. 0. 0.]
      [ 1. 1. 0.]]
    <BLANKLINE>
     [[ 1. 1. 0.]
      [ 0. 1. 0.]
      [ 0. 0. 0.]]]
    """
    x = Coords(pattern(s))
    try:
        return x.reshape(-1, nplex, 3)
    except:
        raise ValueError("Could not reshape points list to plexitude %s" % nplex)


def align(L,align,offset=(0.,0.,0.)):
    """Align a list of geometrical objects.

    Parameters
    ----------
    L: list of Coords or Geometry objects
        A list of objects that have an appropriate ``align`` method,
        like the :class:`Coords` and :class:`~geometry.Geometry`
        (and its subclasses).
    align: str
        A string of three characters, one for each coordinate direction,
        that define how the subsequent objects have to be aligned in each of
        the global axis directions:

            - '-' : align on the minimal coordinate value
            - '+' : align on the maximal coordinate value
            - '0' : align on the middle coordinate value
            - '|' : align the minimum value on the maximal value of the
                    previous item

        Thus the string ``'|--'`` will juxtapose the objects in the x-direction,
        while aligning them on their minimal coordinates in the y- and
        z- direction.
    offset: float array-like (3,), optional
        An extra translation to be given to each subsequent object. This can
        be used to create a space between the objects, instead of
        juxtaposing them.

    Returns
    -------
    list of objects
        A list with the aligned objects.

    Notes
    -----
    See also example Align.

    See Also
    --------
    Coords.align: align a single object with respect to a point.

    """
    r = L[:1]
    al = am =''
    for i in range(3):
        if align[i] == '|':
            al += '-'
            am += '+'
        else:
            al += align[i]
            am += align[i]
    for o in L[1:]:
        r.append(o.align(al, r[-1].bboxPoint(am)+offset))
    return r


### End
