# Source code for plugins.isosurface

    # $Id$
    ##
    ##  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/.
    ##
    #
    """Isosurface: surface reconstruction algorithms

    This module contains the marching cube algorithm.

    Some of the code is based on the example by Paul Bourke from
    http://paulbourke.net/geometry/polygonise/

    """
    from __future__ import absolute_import, division, print_function
    from pyformex import zip

    import numpy as np
    from pyformex.multi import multitask, cpu_count, splitar

    def isosurface(data,level,nproc=-1,tet=0):
        """Create an isosurface through data at given level.

        - `data`: (nx,ny,nz) shaped array of data values at points with
          coordinates equal to their indices. This defines a 3D volume
          [0,nx-1], [0,ny-1], [0,nz-1]
        - `level`: data value at which the isosurface is to be constructed
        - `nproc`: number of parallel processes to use. On multiprocessor machines
          this may be used to speed up the processing. If <= 0 , the number of
          processes will be set equal to the number of processors, to achieve
          a maximal speedup.

        Returns an (ntr,3,3) array defining the triangles of the isosurface.
        The result may be empty (if level is outside the data range).
        """
        if nproc is None:
            nproc = -1
        if nproc < 1:
            nproc = cpu_count()

        if nproc == 1:
            # Perform single process isosurface (accelerated)
            from pyformex.lib import misc
            data = data.astype(np.float32)
            level = np.float32(level)
            tri = misc.isosurface(data, level, tet)

        else:
            # Perform parallel isosurface
            # 1. Split in blocks (and remember shift)
            datablocks = splitar(data, nproc, close=True)
            shift = (np.array([d.shape[0] for d in datablocks]) - 1).cumsum()
            # 2. Solve blocks independently
            tasks = [(isosurface, (d, level, 1, tet)) for d in datablocks]
            tri = multitask(tasks, nproc)
            # 3. Shift and merge blocks
            for t, s in zip(tri[1:], shift[:-1]):
                t[:,:, 2] += s
            tri = np.concatenate(tri, axis=0)

        return tri

    def isoline(data,level,nproc=-1):
        """Create an isocontour through data at given level.

        - `data`: (nx,ny,nz) shaped array of data values at points with
          coordinates equal to their indices. This defines a 2D area
          [0,nx-1], [0,ny-1]
        - `level`: data value at which the isocontour is to be constructed
        - `nproc`: number of parallel processes to use. On multiprocessor machines
          this may be used to speed up the processing. If <= 0 , the number of
          processes will be set equal to the number of processors, to achieve
          a maximal speedup.

        Returns an (nseg,2,2) array defining the 2D coordinates of the
        segments of the isocontour.
        The result may be empty (if level is outside the data range).
        """
        if nproc is None:
            nproc = -1
        if nproc < 1:
            nproc = cpu_count()

        if nproc == 1:
            # Perform single process isoline (accelerated)
            from pyformex.lib import misc
            data = data.astype(np.float32)
            level = np.float32(level)
            seg = misc.isoline(data, level)

        else:
            # Perform parallel isoline
            # 1. Split in blocks (and remember shift)
            datablocks = splitar(data, nproc, close=True)
            shift = (np.array([d.shape[0] for d in datablocks]) - 1).cumsum()
            # 2. Solve blocks independently
            tasks = [(isoline, (d, level, 1)) for d in datablocks]
            seg = multitask(tasks, nproc)
            # 3. Shift and merge blocks
            for t, s in zip(seg[1:], shift[:-1]):
                t[:,:,1] += s
            seg = np.concatenate(seg, axis=0)

        return seg

    # End