levelSetSegmentationObject.cpp

/*---------------------------------------------------------------------------*\
 *
 *  bitpit
 *
 *  Copyright (C) 2015-2021 OPTIMAD engineering Srl
 *
 *  -------------------------------------------------------------------------
 *  License
 *  This file is part of bitpit.
 *
 *  bitpit is free software: you can redistribute it and/or modify it
 *  under the terms of the GNU Lesser General Public License v3 (LGPL)
 *  as published by the Free Software Foundation.
 *
 *  bitpit 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 Lesser General Public
 *  License for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public License
 *  along with bitpit. If not, see <http://www.gnu.org/licenses/>.
 *
\*---------------------------------------------------------------------------*/
 
#include "levelSetSegmentationObject.hpp"
 
#include "bitpit_CG.hpp"
#include "levelSetObject.hpp"
 
namespace bitpit {


const double LevelSetSegmentationSurfaceInfo::DEFAULT_FEATURE_ANGLE = 2. * BITPIT_PI;

LevelSetSegmentationSurfaceInfo::LevelSetSegmentationSurfaceInfo()
    : m_surface(nullptr),
      m_featureAngle(0),
      m_surfaceSmoothing(LevelSetSurfaceSmoothing::LOW_ORDER)
{
}

LevelSetSegmentationSurfaceInfo::LevelSetSegmentationSurfaceInfo(const LevelSetSegmentationSurfaceInfo &other)
    : m_surface(other.m_surface),
      m_featureAngle(other.m_featureAngle),
      m_surfaceSmoothing(other.m_surfaceSmoothing),
      m_segmentVertexOffset(other.m_segmentVertexOffset),
      m_segmentNormalsValid(other.m_segmentNormalsValid),
      m_segmentNormalsStorage(other.m_segmentNormalsStorage),
      m_unlimitedVertexNormalsValid(other.m_unlimitedVertexNormalsValid),
      m_unlimitedVertexNormalsStorage(other.m_unlimitedVertexNormalsStorage),
      m_limitedSegmentVertexNormalValid(other.m_limitedSegmentVertexNormalValid),
      m_limitedSegmentVertexNormalStorage(other.m_limitedSegmentVertexNormalStorage)
{
    if (other.m_ownedSurface) {
        m_ownedSurface = std::unique_ptr<SurfUnstructured>(new SurfUnstructured(*(other.m_ownedSurface)));
    } else {
        m_ownedSurface = nullptr;
    }
 
    if (other.m_searchTree) {
        m_searchTree = std::unique_ptr<SurfaceSkdTree>(new SurfaceSkdTree(m_surface));
        m_searchTree->build();
    } else {
        m_searchTree = nullptr;
    }
}

LevelSetSegmentationSurfaceInfo::LevelSetSegmentationSurfaceInfo(LevelSetSurfaceSmoothing surfaceSmoothing)
    : m_surface(nullptr),
      m_featureAngle(0),
      m_surfaceSmoothing(surfaceSmoothing)
{
}

LevelSetSegmentationSurfaceInfo::LevelSetSegmentationSurfaceInfo(std::unique_ptr<const SurfUnstructured> &&surface, double featureAngle, LevelSetSurfaceSmoothing surfaceSmoothing)
    : m_surfaceSmoothing(surfaceSmoothing)
{
    setSurface(std::move(surface), featureAngle);
}

LevelSetSegmentationSurfaceInfo::LevelSetSegmentationSurfaceInfo(const SurfUnstructured *surface, double featureAngle, LevelSetSurfaceSmoothing surfaceSmoothing)
{
    setSurface(surface, featureAngle, surfaceSmoothing);
}

const SurfUnstructured & LevelSetSegmentationSurfaceInfo::getSurface() const {
    return *m_surface;
}

void LevelSetSegmentationSurfaceInfo::setSurface(std::unique_ptr<const SurfUnstructured> &&surface, double featureAngle){
 
    m_ownedSurface = std::move(surface);
    setSurface(m_ownedSurface.get(), featureAngle);
}

void LevelSetSegmentationSurfaceInfo::setSurface(const SurfUnstructured *surface, double featureAngle, LevelSetSurfaceSmoothing surfaceSmoothing){
 
    // Check if adjacencies are built
    if (surface->getAdjacenciesBuildStrategy() == SurfUnstructured::ADJACENCIES_NONE) {
        throw std::runtime_error ("Segmentation needs adjacencies!") ;
    }
 
    // Surface information
    m_surface          = surface;
    m_featureAngle     = featureAngle;
    m_surfaceSmoothing = surfaceSmoothing;
 
    // Segment vertices information
    m_segmentVertexOffset.unsetKernel();
    m_segmentVertexOffset.setStaticKernel(&m_surface->getCells());
 
    std::size_t nTotalSegmentVertices = 0;
    for (auto itr = m_segmentVertexOffset.begin(); itr != m_segmentVertexOffset.end(); ++itr) {
        *itr = nTotalSegmentVertices;
        nTotalSegmentVertices += m_surface->getCells().rawAt(itr.getRawIndex()).getVertexCount();
    }
 
    // Normals
    m_segmentNormalsValid.unsetKernel();
    m_segmentNormalsValid.setStaticKernel(&m_surface->getCells());
    m_segmentNormalsValid.fill(false);
    m_segmentNormalsStorage.unsetKernel();
    m_segmentNormalsStorage.setStaticKernel(&m_surface->getCells());
 
    m_unlimitedVertexNormalsValid.unsetKernel();
    m_unlimitedVertexNormalsValid.setStaticKernel(&m_surface->getVertices());
    m_unlimitedVertexNormalsValid.fill(false);
    m_unlimitedVertexNormalsStorage.unsetKernel();
    m_unlimitedVertexNormalsStorage.setStaticKernel(&m_surface->getVertices());
 
    m_limitedSegmentVertexNormalValid.assign(nTotalSegmentVertices, false);
 
    // Initialize search tree
    m_searchTree = std::unique_ptr<SurfaceSkdTree>(new SurfaceSkdTree(surface));
    m_searchTree->build();
}

void LevelSetSegmentationSurfaceInfo::setSurface(const SurfUnstructured *surface, double featureAngle){
    setSurface(surface, featureAngle, LevelSetSurfaceSmoothing::LOW_ORDER);
}

const SurfaceSkdTree & LevelSetSegmentationSurfaceInfo::getSearchTree() const {
    return *m_searchTree;
}

double LevelSetSegmentationSurfaceInfo::getFeatureAngle() const {
    return m_featureAngle;
}

bitpit::LevelSetSurfaceSmoothing LevelSetSegmentationSurfaceInfo::getSurfaceSmoothing() const {
    return m_surfaceSmoothing;
}

double LevelSetSegmentationSurfaceInfo::evalDistance(const std::array<double, 3> &point,
                                                     const SegmentConstIterator &segmentItr,
                                                     bool signedDistance,
                                                     std::array<double, 3> *distanceVector) const
{
    // Project the point on the surface and evaluate the point-projection vector
    int nSegmentVertices = segmentItr->getVertexCount();
    BITPIT_CREATE_WORKSPACE(lambda, double, nSegmentVertices, ReferenceElementInfo::MAX_ELEM_VERTICES);
    std::array<double, 3> pointProjection = evalProjection(point, segmentItr, lambda);
    (*distanceVector) = point - pointProjection;
 
    // Evaluate unsigned distance
    double unsignedDistance = norm2(*distanceVector);
    if (!signedDistance) {
        return unsignedDistance;
    }
 
    // Signed distance
    //
    // If the sign is null and the point doesn't lie on the segmentation, it lies on the normal
    // plane. This case is not supported, because it would require to evaluate the sign taking
    // into account the the curvature of the surface.
    std::array<double, 3> pseudoNormal = computePseudoNormal(segmentItr, lambda);
    double pointProjectionNormalComponent = dotProduct(*distanceVector, pseudoNormal);
 
    double distanceTolerance = m_surface->getTol();
    if (utils::DoubleFloatingEqual()(pointProjectionNormalComponent, 0., distanceTolerance, distanceTolerance)) {
        bool pointOnSegmentation = utils::DoubleFloatingEqual()(unsignedDistance, 0., distanceTolerance, distanceTolerance);
        if (!pointOnSegmentation) {
            throw std::runtime_error("Unable to evaluate point sign: the point lies on the normal plane!");
        }
    }
 
    return sign(pointProjectionNormalComponent) * unsignedDistance;
}

double LevelSetSegmentationSurfaceInfo::evalDistance(const std::array<double, 3> &point,
                                                     const SegmentConstIterator &segmentItr,
                                                     bool signedDistance) const
{
    std::array<double, 3> distanceVector;
    return evalDistance(point, segmentItr, signedDistance, &distanceVector);
}

std::array<double, 3> LevelSetSegmentationSurfaceInfo::evalDistanceVector(const std::array<double, 3> &point,
                                                                          const SegmentConstIterator &segmentItr) const
{
    int nSegmentVertices = segmentItr->getVertexCount();
    BITPIT_CREATE_WORKSPACE(lambda, double, nSegmentVertices, ReferenceElementInfo::MAX_ELEM_VERTICES);
    std::array<double, 3> pointProjection = evalProjection(point, segmentItr, lambda);
 
    return (point - pointProjection);
}

std::array<double, 3> LevelSetSegmentationSurfaceInfo::evalNormal(const std::array<double, 3> &point,
                                                                  const SegmentConstIterator &segmentItr) const
{
    // Project the point on the surface and evaluate the point-projection vector
    std::array<double, 3> projectionPoint;
    std::array<double, 3> projectionNormal;
    evalProjection(point, segmentItr, &projectionPoint, &projectionNormal);
 
    return projectionNormal;
}

std::array<double,3> LevelSetSegmentationSurfaceInfo::evalPseudoNormal(const SegmentConstIterator &segmentItr,
                                                                       const double *lambda ) const
{
    return computePseudoNormal(segmentItr, lambda);
}

void LevelSetSegmentationSurfaceInfo::evalProjectionOnVertex(const std::array<double,3> &point,
                                                             const SegmentConstIterator &segmentItr,
                                                             std::array<double, 3> *projectionPoint,
                                                             std::array<double, 3> *projectionNormal) const
{
    BITPIT_UNUSED(point);
 
    // Get segment
    const Cell &segment = *segmentItr;
    assert(segment.getType() == ElementType::VERTEX);
 
    // Get vertex id
    ConstProxyVector<long> segmentVertexIds = segment.getVertexIds();
    long id = segmentVertexIds[0];
 
    // Compute projection point and normal 
    (*projectionPoint)  = m_surface->getVertexCoords(id);
    (*projectionNormal) = {0., 0., 0.};
}

void LevelSetSegmentationSurfaceInfo::evalHighOrderProjectionOnLine(const std::array<double,3> &point,
                                                                    const SegmentConstIterator &segmentItr,
                                                                    std::array<double, 3> *projectionPoint,
                                                                    std::array<double, 3> *projectionNormal) const
{
    // Get segment
    const Cell &segment = *segmentItr;
    assert(segment.getType() == ElementType::LINE);
 
    // Get segment vertices
    ConstProxyVector<long> segmentVertexIds = segment.getVertexIds();
 
    const std::array<double,3> &point0 = m_surface->getVertexCoords(segmentVertexIds[0]);
    const std::array<double,3> &point1 = m_surface->getVertexCoords(segmentVertexIds[1]);
 
    // Get projection point on segment
    std::array<double,2> t;
    std::array<double, 3> point_s = CGElem::projectPointSegment(point, point0, point1, &(t[0]));
 
    // Early return if point_s consides with segment's node
    double distanceTolerance = m_surface->getTol();
    for (int i = 0; i < 2; ++i) {
        if (utils::DoubleFloatingEqual()(t[i], 1., distanceTolerance, distanceTolerance)) {
            (*projectionPoint)  = m_surface->getVertexCoords(segmentVertexIds[i]);
            (*projectionNormal) = computeSegmentVertexNormal(segmentItr, i, true);
            return;
        }
    }
 
    // Get normal on segment
    std::array<double, 3> facetNormal = computeSegmentNormal(segmentItr);
    std::array<double, 3> normal_s = point - point_s;
 
    double distance = norm2(normal_s);
    if (utils::DoubleFloatingEqual()(distance, 0., distanceTolerance, distanceTolerance)) {
        normal_s = facetNormal;
    } else {
        normal_s /= distance;
        if (dotProduct(normal_s, facetNormal) < 0.0) {
            normal_s *= -1.0;
        }
    }
 
    // Get normals on segment's vertices
    std::array<double, 3> normal0 = computeSegmentVertexNormal(segmentItr, 0, true);
    std::array<double, 3> normal1 = computeSegmentVertexNormal(segmentItr, 1, true);
 
    // Compute lambda coefficients
    std::array<double, 3> edge = point1 - point0;
    double lambda01 = -dotProduct(edge, normal0) / dotProduct(normal_s, normal0); 
    double lambda10 =  dotProduct(edge, normal1) / dotProduct(normal_s, normal1); 
 
    // Eval polynomial parameterizing the curve
    double product1 = t[0] * t[1];
    double product2 = lambda01 * t[0] + lambda10 * t[1];
    double f        = product1 * product2;
    double df       = product1 * (lambda01 - lambda10) + product2 * (t[1] - t[0]);
 
    // Eval projection point on curve
    (*projectionPoint) = point_s + f * normal_s;
 
    // Eval normal vector on curve
    double length = norm2(edge);
    (*projectionNormal)  = df * edge + length * length * normal_s;
    (*projectionNormal) /= norm2((*projectionNormal));
 
    if (dotProduct(normal_s, (*projectionNormal)) < 0.0) {
       (*projectionNormal) *= -1.0;
    }
}

void LevelSetSegmentationSurfaceInfo::evalHighOrderProjectionOnTriangle(const std::array<double,3> &point,
                                                                        const SegmentConstIterator &segmentItr,
                                                                        std::array<double, 3> *projectionPoint,
                                                                        std::array<double, 3> *projectionNormal) const
{
    // Get segment
    const Cell &segment = *segmentItr;
    assert(segment.getType() == ElementType::TRIANGLE);
 
    // Get segment vertices
    ConstProxyVector<long> segmentVertexIds = segment.getVertexIds();
 
    const std::array<double,3> &point0 = m_surface->getVertexCoords(segmentVertexIds[0]);
    const std::array<double,3> &point1 = m_surface->getVertexCoords(segmentVertexIds[1]);
    const std::array<double,3> &point2 = m_surface->getVertexCoords(segmentVertexIds[2]);
 
    // Get projection point on segment
    std::array<double,3> tau;
    std::array<double, 3> point_s = CGElem::projectPointTriangle(point, point0, point1, point2, &(tau[0]));
    std::array<double, 3> d0_point_s = point0 - point2;
    std::array<double, 3> d1_point_s = point1 - point2;
 
    // Early return if point_s consides with segment's node
    double distanceTolerance = m_surface->getTol();
    for (int i = 0; i < 3; ++i) {
        if (utils::DoubleFloatingEqual()(tau[i], 1., distanceTolerance, distanceTolerance)) {
            (*projectionPoint)  = m_surface->getVertexCoords(segmentVertexIds[i]);
            (*projectionNormal) = computeSegmentVertexNormal(segmentItr, i, true);
            return;
        }
    }
 
    // Get normal on segment
    std::array<double, 3> facetNormal = computeSegmentNormal(segmentItr);
    std::array<double, 3> normal_s = point - point_s;
 
    double distance = norm2(normal_s);
    if (utils::DoubleFloatingEqual()(distance, 0., distanceTolerance, distanceTolerance)) {
        normal_s = facetNormal;
    } else {
        normal_s /= distance;
        if (dotProduct(normal_s, facetNormal) < 0.0) {
            normal_s *= -1.0;
        }
    }
 
    // Compute projection point and normal at each edge "e"
    std::array<std::array<double, 3>,3 > edgePoint_s;
    std::array< std::array<double, 3>, 3 > d0_edgePoint_s;
    std::array< std::array<double, 3>, 3 > d1_edgePoint_s;
    std::array< std::array<double, 3>, 3 > edgePoint_p;
    std::array< std::array<double, 3>, 3 > d0_edgePoint_p;
    std::array< std::array<double, 3>, 3 > d1_edgePoint_p;
    std::array< std::array<double, 3>, 3 > edgeNormal_p;
    std::array< std::array<double, 3>, 3 > d0_edgeNormal_p;
    std::array< std::array<double, 3>, 3 > d1_edgeNormal_p;
    for (int e = 0; e < 3; ++e) {
        // Get segment
        const Cell &segment = *segmentItr;
 
        // Get edge's local and global node ids
        ConstProxyVector<int> edgeLocalVertexIds = segment.getFaceLocalVertexIds(e);
        int localId0 = edgeLocalVertexIds[0]; 
        int localId1 = edgeLocalVertexIds[1]; 
        int id0      = segment.getVertexId(localId0);
        int id1      = segment.getVertexId(localId1);
 
        // Get vertex information
        const std::array<double,3> &node0 = m_surface->getVertexCoords(id0);
        const std::array<double,3> &node1 = m_surface->getVertexCoords(id1);
 
        std::array<double, 3> normal0 = computeSegmentVertexNormal(segmentItr, localId0, true);
        std::array<double, 3> normal1 = computeSegmentVertexNormal(segmentItr, localId1, true);
 
        // Get edge normal
        std::array<double, 3> edgeNormal_s = computeSegmentEdgeNormal(segmentItr, e, true);
 
        // Compute projection coefficients
        std::array<double, 3> edge = node1 - node0;
        double p01 = dotProduct(point_s - node0, edge) / dotProduct(edge, edge);
        double d0_p01;
        double d1_p01;
        if (p01 > 1.0) {
            p01 = 1.0;
            d0_p01 = 0.0;
            d1_p01 = 0.0;
        } else if (p01 < 0.0) {
            p01 = 0.0;
            d0_p01 = 0.0;
            d1_p01 = 0.0;
        } else {
            d0_p01 = dotProduct(d0_point_s - node0, edge) / dotProduct(edge, edge);
            d1_p01 = dotProduct(d1_point_s - node0, edge) / dotProduct(edge, edge);
        }
        double p10 = 1.0 - p01;
        double d0_p10 = - d0_p01;
        double d1_p10 = - d1_p01;
 
        // Compute lambda coefficients
        double lambda01 = -dotProduct(edge, normal0) / dotProduct(edgeNormal_s, normal0); 
        double lambda10 =  dotProduct(edge, normal1) / dotProduct(edgeNormal_s, normal1); 
 
        // Edge projection point
        double product1   = p10 * p01;
        double d_product1 = p01 - p10; // derivative wrt p10
 
        double product2   = lambda01 * p10 + lambda10 * p01;
        double d_product2 = lambda01 - lambda10;
 
        double f01     = product1 * product2;
        double d_f01   = d_product1 * product2 + product1 * d_product2;
        double d_d_f01 = 2.0 * (d_product1 * d_product2 - product2);
 
        edgePoint_s[e] = p10 * node0 + p01 * node1;
        std::array<double, 3> d_edgePoint_s = node0 - node1;
 
        d0_edgePoint_s[e] = d_edgePoint_s  * d0_p10;
        d1_edgePoint_s[e] = d_edgePoint_s  * d1_p10;
 
        edgePoint_p[e] = edgePoint_s[e] + f01 * edgeNormal_s;
        std::array<double, 3> d_edgePoint_p = d_edgePoint_s + d_f01 * edgeNormal_s;
 
        d0_edgePoint_p[e] = d_edgePoint_p * d0_p10;
        d1_edgePoint_p[e] = d_edgePoint_p * d1_p10;
 
        // Edge surface normal
        double length   = norm2(edge);
        edgeNormal_p[e] = d_f01 * edge + length * length * edgeNormal_s;
        std::array<double, 3> d_edgeNormal_p = d_d_f01 * edge;
 
        double norm = norm2(edgeNormal_p[e]);
        double d_norm = dotProduct(edgeNormal_p[e], d_edgeNormal_p) / norm;
        edgeNormal_p[e] /= norm;
        d_edgeNormal_p = (d_edgeNormal_p - edgeNormal_p[e] * d_norm) / norm;
 
        d0_edgeNormal_p[e] = d_edgeNormal_p * d0_p10;
        d1_edgeNormal_p[e] = d_edgeNormal_p * d1_p10;
    }
 
    // Project surface point to the intermediator polygon
    double sum_t = 0.0;
    double d0_sum_t = 0.0;
    double d1_sum_t = 0.0;
    std::array<double, 3> t;
    std::array<double, 3> d0_t;
    std::array<double, 3> d1_t;
    for (int i = 0; i < 3; i ++) { // loop in vertices
        int j = (i + 1) % 3;
        int k = (i + 2) % 3;
        std::array<double, 3> v1 = edgePoint_s[j] - point_s;
        std::array<double, 3> d0_v1 = d0_edgePoint_s[j]- d0_point_s;
        std::array<double, 3> d1_v1 = d1_edgePoint_s[j]- d1_point_s;
 
        std::array<double, 3> v2 = edgePoint_s[k] - point_s;
        std::array<double, 3> d0_v2 = d0_edgePoint_s[k]- d0_point_s;
        std::array<double, 3> d1_v2 = d1_edgePoint_s[k]- d1_point_s;
 
        std::array<double, 3> normal = crossProduct(v1, v2);
        std::array<double, 3> d0_normal = crossProduct(d0_v1, v2) + crossProduct(v1, d0_v2);
        std::array<double, 3> d1_normal = crossProduct(d1_v1, v2) + crossProduct(v1, d1_v2);
 
        t[i] = norm2(normal);
        d0_t[i] = dotProduct(normal, d0_normal) / t[i];
        d1_t[i] = dotProduct(normal, d1_normal) / t[i];
        sum_t += t[i];
        d0_sum_t += d0_t[i];
        d1_sum_t += d1_t[i];
    }
 
    std::array<double, 3> point_m    = {0.0, 0.0, 0.0};
    std::array<double, 3> d0_point_m = {0.0, 0.0, 0.0};
    std::array<double, 3> d1_point_m = {0.0, 0.0, 0.0};
    for (int i = 0; i < 3; i ++) {
        t[i] /= sum_t;
        d0_t[i] = (d0_t[i] - t[i] * d0_sum_t) / sum_t;
        d1_t[i] = (d1_t[i] - t[i] * d1_sum_t) / sum_t;
        point_m += t[i] * edgePoint_p[i];
        d0_point_m += d0_t[i] * edgePoint_p[i] + t[i] * d0_edgePoint_p[i];
        d1_point_m += d1_t[i] * edgePoint_p[i] + t[i] * d1_edgePoint_p[i];
    }
 
    // Early return if point_s coinsides with a triangle
    // edge (or a node of the intermediate triangle)
    for (int i = 0; i < 3; ++i) {
        if (utils::DoubleFloatingEqual()(t[i], 1., distanceTolerance, distanceTolerance)) {
            (*projectionPoint)  = edgePoint_p[i];
            (*projectionNormal) = edgeNormal_p[i];
            return;
        }
    }
 
    // Compute 2D f function
    double f    = 0.0;
    double d0_f = 0.0;
    double d1_f = 0.0;
    for (int e = 0; e < 3; ++e) {
        int i = e;           // previous vertex
        int j = (e + 1) % 3; // next vertex
 
        // Get vertex information
        std::array<double,3> &node_i    = edgePoint_p[i];
        std::array<double,3> &d0_node_i = d0_edgePoint_p[i];
        std::array<double,3> &d1_node_i = d1_edgePoint_p[i];
 
        std::array<double,3> &node_j     = edgePoint_p[j];
        std::array<double, 3> &d0_node_j = d0_edgeNormal_p[j];
        std::array<double, 3> &d1_node_j = d1_edgeNormal_p[j];
 
        std::array<double, 3> &normal_i    = edgeNormal_p[i];
        std::array<double, 3> &d0_normal_i = d0_edgeNormal_p[i];
        std::array<double, 3> &d1_normal_i = d1_edgeNormal_p[i];
 
        std::array<double, 3> &normal_j    = edgeNormal_p[j];
        std::array<double, 3> &d0_normal_j = d0_edgeNormal_p[j];
        std::array<double, 3> &d1_normal_j = d1_edgeNormal_p[j];
 
        // Compute lambda parameters
        std::array<double, 3> edge    = node_j - node_i;
        std::array<double, 3> d0_edge = d0_node_j - d0_node_i;
        std::array<double, 3> d1_edge = d1_node_j - d1_node_i;
 
        double lambda_ij    = -dotProduct(edge, normal_i) / dotProduct(normal_s, normal_i);
        double d0_lambda_ij = -(dotProduct(edge - lambda_ij * normal_s, d0_normal_i) + dotProduct(d0_edge, normal_i)) / dotProduct(normal_s, normal_i);
        double d1_lambda_ij = -(dotProduct(edge - lambda_ij * normal_s, d1_normal_i) + dotProduct(d1_edge, normal_i)) / dotProduct(normal_s, normal_i);
 
        double lambda_ji    = dotProduct(edge, normal_j) / dotProduct(normal_s, normal_j);
        double d0_lambda_ji = (dotProduct(edge - lambda_ji * normal_s, d0_normal_j) + dotProduct(d0_edge, normal_j)) / dotProduct(normal_s, normal_j);
        double d1_lambda_ji = (dotProduct(edge - lambda_ji * normal_s, d1_normal_j) + dotProduct(d1_edge, normal_j)) / dotProduct(normal_s, normal_j);
 
        // Compute contribution to the f function
        f += t[i] * t[j] * (lambda_ij * t[i] + lambda_ji * t[j]);
        d0_f += (d0_t[i] * t[j] + t[i] * d0_t[j]) * (lambda_ij * t[i] + lambda_ji * t[j])
              + t[i] * t[j] * (d0_lambda_ij * t[i] + lambda_ij * d0_t[i] + d0_lambda_ji * t[j] + lambda_ji * d0_t[j]);
        d1_f += (d1_t[i] * t[j] + t[i] * d1_t[j]) * (lambda_ij * t[i] + lambda_ji * t[j])
              + t[i] * t[j] * (d1_lambda_ij * t[i] + lambda_ij * d1_t[i] + d1_lambda_ji * t[j] + lambda_ji * d1_t[j]);
    }
 
    // Compute projection point on curved surface
    (*projectionPoint) = point_m + f * normal_s;
 
    // Compute projection normal on curved surface
    std::array<double, 3> d0_point_p = d0_point_m + d0_f * normal_s;
    std::array<double, 3> d1_point_p = d1_point_m + d1_f * normal_s;
 
    (*projectionNormal) = crossProduct(d0_point_p, d1_point_p);
    (*projectionNormal) /= norm2((*projectionNormal));
 
    if (dotProduct(normal_s, (*projectionNormal)) < 0.0) {
       (*projectionNormal) *= -1.0;
    }
}

void LevelSetSegmentationSurfaceInfo::evalHighOrderProjectionOnPolygon(const std::array<double,3> &point,
                                                                       const SegmentConstIterator &segmentItr,
                                                                       std::array<double, 3> *projectionPoint,
                                                                       std::array<double, 3> *projectionNormal) const
{
    // Get segment
    const Cell &segment = *segmentItr;
 
    // Get projection point on segment
    ConstProxyVector<long> segmentVertexIds = m_surface->getFacetOrderedVertexIds(segment);
 
    std::size_t nSegmentVertices = segmentVertexIds.size();
    BITPIT_CREATE_WORKSPACE(segmentVertexCoords, std::array<double BITPIT_COMMA 3>, nSegmentVertices, ReferenceElementInfo::MAX_ELEM_VERTICES);
    m_surface->getVertexCoords(segmentVertexIds.size(), segmentVertexIds.data(), segmentVertexCoords);
 
    BITPIT_CREATE_WORKSPACE(tau, double, nSegmentVertices, ReferenceElementInfo::MAX_ELEM_VERTICES);
    std::array<double, 3> point_s = CGElem::projectPointPolygon(point, nSegmentVertices, segmentVertexCoords, tau);
 
    std::array<double, 3> lowOrderProjectionNormal = tau[0] * computeSegmentVertexNormal(segmentItr, 0, true);
    for (std::size_t i = 1; i < nSegmentVertices; ++i) {
        lowOrderProjectionNormal += tau[i] * computeSegmentVertexNormal(segmentItr, i, true);
    }
 
    // Early return if point_s consides with segment's node
    double distanceTolerance = m_surface->getTol();
    for (std::size_t i = 0; i < nSegmentVertices; ++i) {
        if (utils::DoubleFloatingEqual()(tau[i], 1., distanceTolerance, distanceTolerance)) {
            (*projectionPoint)  = segmentVertexCoords[i];
            (*projectionNormal) = computeSegmentVertexNormal(segmentItr, i, true);
            return;
        }
    }
 
    std::array<double, 3> d0_point_s;
    std::array<double, 3> d1_point_s;
 
    // Define parameterization of point_s
    // Definition: point_s - x_central = t0 * (x_prev - x_central) + t1 * (x_next - x_central)
    // where x_prev, x_central, x_next are subsequent polygon vertices
    // Vertex x_central is the one corresponding to polygon angle closer to 90 degrees
    double cosMin = std::numeric_limits<double>::max();
    for (std::size_t i = 0; i < nSegmentVertices; ++i) {
        int i_p = i;                          // previous
        int i_c = (i + 1) % nSegmentVertices; // central
        int i_n = (i + 2) % nSegmentVertices; // next
 
        const std::array<double,3> &point_p = m_surface->getVertexCoords(segmentVertexIds[i_p]);
        const std::array<double,3> &point_c = m_surface->getVertexCoords(segmentVertexIds[i_c]);
        const std::array<double,3> &point_n = m_surface->getVertexCoords(segmentVertexIds[i_n]);
 
        std::array<double,3> x_p = point_p - point_c;
        std::array<double,3> x_n = point_n - point_c;
 
        double cos = std::abs(dotProduct(x_p, x_n)) / (norm2(x_p) * norm2(x_n));
        if (cos < cosMin) {
            cosMin = cos;
            d0_point_s = x_p;
            d1_point_s = x_n;
        }
    }
 
    // Get normal on segment
    std::array<double, 3> facetNormal = computeSegmentNormal(segmentItr);
    std::array<double, 3> normal_s = point - point_s;
    double distance = norm2(normal_s);
    if (utils::DoubleFloatingEqual()(distance, 0., distanceTolerance, distanceTolerance)) {
        normal_s = facetNormal;
    } else {
        normal_s /= distance;
        if (dotProduct(normal_s, facetNormal) < 0.0) {
            normal_s *= -1.0;
        }
    }
 
    // Compute projection point and normal at each edge "e"
    int nSegmentEdges = nSegmentVertices;
    BITPIT_CREATE_WORKSPACE(edgePoint_s, std::array<double BITPIT_COMMA 3>, nSegmentEdges, ReferenceElementInfo::MAX_ELEM_EDGES);
    BITPIT_CREATE_WORKSPACE(d0_edgePoint_s, std::array<double BITPIT_COMMA 3>, nSegmentEdges, ReferenceElementInfo::MAX_ELEM_EDGES);
    BITPIT_CREATE_WORKSPACE(d1_edgePoint_s, std::array<double BITPIT_COMMA 3>, nSegmentEdges, ReferenceElementInfo::MAX_ELEM_EDGES);
    BITPIT_CREATE_WORKSPACE(edgePoint_p, std::array<double BITPIT_COMMA 3>, nSegmentEdges, ReferenceElementInfo::MAX_ELEM_EDGES);
    BITPIT_CREATE_WORKSPACE(d0_edgePoint_p, std::array<double BITPIT_COMMA 3>, nSegmentEdges, ReferenceElementInfo::MAX_ELEM_EDGES);
    BITPIT_CREATE_WORKSPACE(d1_edgePoint_p, std::array<double BITPIT_COMMA 3>, nSegmentEdges, ReferenceElementInfo::MAX_ELEM_EDGES);
    BITPIT_CREATE_WORKSPACE(edgeNormal_p, std::array<double BITPIT_COMMA 3>, nSegmentEdges, ReferenceElementInfo::MAX_ELEM_EDGES);
    BITPIT_CREATE_WORKSPACE(d0_edgeNormal_p, std::array<double BITPIT_COMMA 3>, nSegmentEdges, ReferenceElementInfo::MAX_ELEM_EDGES);
    BITPIT_CREATE_WORKSPACE(d1_edgeNormal_p, std::array<double BITPIT_COMMA 3>, nSegmentEdges, ReferenceElementInfo::MAX_ELEM_EDGES);
    for (int e = 0; e < nSegmentEdges; ++e) {
        // Get edge's local and global node ids
        ConstProxyVector<int> edgeLocalVertexIds = segment.getFaceLocalVertexIds(e);
        int localId0 = edgeLocalVertexIds[0]; 
        int localId1 = edgeLocalVertexIds[1]; 
        int id0      = segment.getVertexId(localId0);
        int id1      = segment.getVertexId(localId1);
 
        // Get vertex information
        const std::array<double,3> &node0 = m_surface->getVertexCoords(id0);
        const std::array<double,3> &node1 = m_surface->getVertexCoords(id1);
 
        std::array<double, 3> normal0 = computeSegmentVertexNormal(segmentItr, localId0, true);
        std::array<double, 3> normal1 = computeSegmentVertexNormal(segmentItr, localId1, true);
 
        // Get edge normal
        std::array<double, 3> edgeNormal_s = computeSegmentEdgeNormal(segmentItr, e, true);
 
        // Compute projection coefficients
        std::array<double, 3> edge = node1 - node0;
        double p01 = dotProduct(point_s - node0, edge) / dotProduct(edge, edge);
        double d0_p01;
        double d1_p01;
        if (p01 > 1.0) {
            p01 = 1.0;
            d0_p01 = 0.0;
            d1_p01 = 0.0;
        } else if (p01 < 0.0) {
            p01 = 0.0;
            d0_p01 = 0.0;
            d1_p01 = 0.0;
        } else {
            d0_p01 = dotProduct(d0_point_s - node0, edge) / dotProduct(edge, edge);
            d1_p01 = dotProduct(d1_point_s - node0, edge) / dotProduct(edge, edge);
        }
        double p10 = 1.0 - p01;
        double d0_p10 = - d0_p01;
        double d1_p10 = - d1_p01;
 
        // Compute lambda coefficients
        double lambda01 = -dotProduct(edge, normal0) / dotProduct(edgeNormal_s, normal0); 
        double lambda10 =  dotProduct(edge, normal1) / dotProduct(edgeNormal_s, normal1); 
 
        // Edge projection point
        double product1 = p10 * p01;
        double d_product1 = p01 - p10; // derivative wrt p10
 
        double product2 = lambda01 * p10 + lambda10 * p01;
        double d_product2 = lambda01 - lambda10;
 
        double f01     = product1 * product2;
        double d_f01   = d_product1 * product2 + product1 * d_product2;
        double d_d_f01 = 2.0 * (d_product1 * d_product2 - product2);
 
        edgePoint_s[e] = p10 * node0 + p01 * node1;
        std::array<double, 3> d_edgePoint_s = node0 - node1;
 
        d0_edgePoint_s[e] = d_edgePoint_s  * d0_p10;
        d1_edgePoint_s[e] = d_edgePoint_s  * d1_p10;
 
        edgePoint_p[e] = edgePoint_s[e] + f01 * edgeNormal_s;
        std::array<double, 3> d_edgePoint_p = d_edgePoint_s + d_f01 * edgeNormal_s;
 
        d0_edgePoint_p[e] = d_edgePoint_p * d0_p10;
        d1_edgePoint_p[e] = d_edgePoint_p * d1_p10;
 
        // Edge surface normal
        double length   = norm2(edge);
        edgeNormal_p[e] = d_f01 * edge + length * length * edgeNormal_s;
        std::array<double, 3> d_edgeNormal_p = d_d_f01 * edge;
 
        double norm = norm2(edgeNormal_p[e]);
        double d_norm = dotProduct(edgeNormal_p[e], d_edgeNormal_p) / norm;
        edgeNormal_p[e] /= norm;
        d_edgeNormal_p = (d_edgeNormal_p - edgeNormal_p[e] * d_norm) / norm;
 
        d0_edgeNormal_p[e] = d_edgeNormal_p * d0_p10;
        d1_edgeNormal_p[e] = d_edgeNormal_p * d1_p10;
    }
 
    // Early return if point_s coinsides with a polygon's
    // edge (or a node of the intermediate polygon)
    for (int e = 0; e < nSegmentEdges; ++e) {
        double distance = norm2(edgePoint_s[e] - point_s);
        if (utils::DoubleFloatingEqual()(distance, 0.0, distanceTolerance, distanceTolerance)) {
            (*projectionPoint)  = edgePoint_p[e];
            (*projectionNormal) = edgeNormal_p[e];
            return;
        }
    }
 
    // Project surface point to the intermediator polygon
    // The projection is done by following formula:
    // point_m - point_s = sum i=0,vertices-1 { t[i] * (edgePoint_p[i] - edgePoint_s[i]) }
    // After various numerical experiments it is concluded that a quite smooth interpolation is
    // resulted from the following choice of the interpolation weights "t":
    // t[i] = w[i] / sum j=0,vertices-1 { w[j] }, where
    // w[i] = product j=0,vertices-1, j!=i { d[j]^4 }, where
    // d[i] is the distance between point_s and edgePoint_s[i]
    BITPIT_CREATE_WORKSPACE(d, double, nSegmentVertices, ReferenceElementInfo::MAX_ELEM_VERTICES);
    BITPIT_CREATE_WORKSPACE(d0_d, double, nSegmentVertices, ReferenceElementInfo::MAX_ELEM_VERTICES);
    BITPIT_CREATE_WORKSPACE(d1_d, double, nSegmentVertices, ReferenceElementInfo::MAX_ELEM_VERTICES);
    for (std::size_t i = 0; i < nSegmentVertices; ++i) {
        d[i] = (edgePoint_s[i][0] - point_s[0]) * (edgePoint_s[i][0] - point_s[0])
        + (edgePoint_s[i][1] - point_s[1]) * (edgePoint_s[i][1] - point_s[1])
        + (edgePoint_s[i][2] - point_s[2]) * (edgePoint_s[i][2] - point_s[2]);
 
        d[i] = d[i] * d[i];
 
        d0_d[i] = 4.0 * d[i] * dotProduct((edgePoint_s[i] - point_s), (d0_edgePoint_s[i] - d0_point_s));
        d1_d[i] = 4.0 * d[i] * dotProduct((edgePoint_s[i] - point_s), (d1_edgePoint_s[i] - d1_point_s));
    }
 
    BITPIT_CREATE_WORKSPACE(t, double, nSegmentVertices, ReferenceElementInfo::MAX_ELEM_VERTICES);
    BITPIT_CREATE_WORKSPACE(d0_t, double, nSegmentVertices, ReferenceElementInfo::MAX_ELEM_VERTICES);
    BITPIT_CREATE_WORKSPACE(d1_t, double, nSegmentVertices, ReferenceElementInfo::MAX_ELEM_VERTICES);
    double sum = 0.0;
    double d0_sum = 0.0;
    double d1_sum = 0.0;
    for (std::size_t i = 0; i < nSegmentVertices; ++i) {
        t[i]    = 1.0;
        d0_t[i] = 0.0;
        d1_t[i] = 0.0;
        for (std::size_t j = 0; j < nSegmentVertices; ++j) {
            if (i != j) {
                t[i] *= d[j];
 
                double product = 1.0;
                for (std::size_t k = 0; k < nSegmentVertices; ++k) {
                    if (j != k) {
                        product *= d[k];
                    }
                }
                d0_t[i] += d0_d[j] * product;
                d1_t[i] += d1_d[j] * product;
            }
        }
        sum    += t[i];
        d0_sum += d0_t[i];
        d1_sum += d1_t[i];
    }
    for (std::size_t i = 0; i < nSegmentVertices; ++i) {
        t[i] /= sum;
        d0_t[i] = (d0_t[i] - t[i] * d0_sum) / sum;
        d1_t[i] = (d1_t[i] - t[i] * d1_sum) / sum;
    }
 
    std::array<double, 3> point_m = {0.0, 0.0, 0.0};
    std::array<double, 3> d0_point_m = {0.0, 0.0, 0.0};
    std::array<double, 3> d1_point_m = {0.0, 0.0, 0.0};
    for (std::size_t i = 0; i < nSegmentVertices; ++i) {
        point_m += (edgePoint_p[i] - edgePoint_s[i]) * t[i];
        d0_point_m += d0_t[i] * (edgePoint_p[i]- edgePoint_s[i]) + t[i] * (d0_edgePoint_p[i] - d0_edgePoint_s[i]);
        d1_point_m += d1_t[i] * (edgePoint_p[i]- edgePoint_s[i]) + t[i] * (d1_edgePoint_p[i] - d1_edgePoint_s[i]);
    }
    point_m += point_s;
    d0_point_m += d0_point_s;
    d1_point_m += d1_point_s;
 
    // Compute 2D f function
    double f    = 0.0;
    double d0_f = 0.0;
    double d1_f = 0.0;
    for (int e = 0; e < nSegmentEdges; ++e) {
        int i = e;                       // previous vertex
        int j = (e + 1) % nSegmentEdges; // next vertex
 
        // Get vertex information
        std::array<double,3> &node_i    = edgePoint_p[i];
        std::array<double,3> &d0_node_i = d0_edgePoint_p[i];
        std::array<double,3> &d1_node_i = d1_edgePoint_p[i];
 
        std::array<double, 3> &node_j     = edgePoint_p[j];
        std::array<double, 3> &d0_node_j = d0_edgeNormal_p[j];
        std::array<double, 3> &d1_node_j = d1_edgeNormal_p[j];
 
        std::array<double, 3> &normal_i    = edgeNormal_p[i];
        std::array<double, 3> &d0_normal_i = d0_edgeNormal_p[i];
        std::array<double, 3> &d1_normal_i = d1_edgeNormal_p[i];
 
        std::array<double, 3> &normal_j    = edgeNormal_p[j];
        std::array<double, 3> &d0_normal_j = d0_edgeNormal_p[j];
        std::array<double, 3> &d1_normal_j = d1_edgeNormal_p[j];
 
        // Compute lambda parameters
        std::array<double, 3> edge    = node_j - node_i;
        std::array<double, 3> d0_edge = d0_node_j - d0_node_i;
        std::array<double, 3> d1_edge = d1_node_j - d1_node_i;
 
        double lambda_ij    = -dotProduct(edge, normal_i) / dotProduct(normal_s, normal_i);
        double d0_lambda_ij = -(dotProduct(edge - lambda_ij * normal_s, d0_normal_i) + dotProduct(d0_edge, normal_i)) / dotProduct(normal_s, normal_i);
        double d1_lambda_ij = -(dotProduct(edge - lambda_ij * normal_s, d1_normal_i) + dotProduct(d1_edge, normal_i)) / dotProduct(normal_s, normal_i);
 
        double lambda_ji    =  dotProduct(edge, normal_j) / dotProduct(normal_s, normal_j);
        double d0_lambda_ji = (dotProduct(edge - lambda_ji * normal_s, d0_normal_j) + dotProduct(d0_edge, normal_j)) / dotProduct(normal_s, normal_j);
        double d1_lambda_ji = (dotProduct(edge - lambda_ji * normal_s, d1_normal_j) + dotProduct(d1_edge, normal_j)) / dotProduct(normal_s, normal_j);
 
        // Compute contribution to the f function
        f += t[i] * t[j] * (lambda_ij * t[i] + lambda_ji * t[j]);
        d0_f += (d0_t[i] * t[j] + t[i] * d0_t[j]) * (lambda_ij * t[i] + lambda_ji * t[j])
              + t[i] * t[j] * (d0_lambda_ij * t[i] + lambda_ij * d0_t[i] + d0_lambda_ji * t[j] + lambda_ji * d0_t[j]);
        d1_f += (d1_t[i] * t[j] + t[i] * d1_t[j]) * (lambda_ij * t[i] + lambda_ji * t[j])
              + t[i] * t[j] * (d1_lambda_ij * t[i] + lambda_ij * d1_t[i] + d1_lambda_ji * t[j] + lambda_ji * d1_t[j]);
    }
 
    // Compute projection point on curved surface
    (*projectionPoint) = point_m + f * normal_s;
 
    // Compute projection normal on curved surface
    std::array<double, 3> d0_point_p = d0_point_m + d0_f * normal_s;
    std::array<double, 3> d1_point_p = d1_point_m + d1_f * normal_s;
 
    (*projectionNormal) = crossProduct(d0_point_p, d1_point_p);
    (*projectionNormal) /= norm2((*projectionNormal));
 
    if (dotProduct(normal_s, (*projectionNormal)) < 0.0) {
       (*projectionNormal) *= -1.0;
    }
}

void LevelSetSegmentationSurfaceInfo::evalHighOrderProjection(const std::array<double,3> &point,
                                                              const SegmentConstIterator &segmentItr,
                                                              std::array<double, 3> *projectionPoint,
                                                              std::array<double, 3> *projectionNormal) const
{
 
    const Cell &segment = *segmentItr;
    ElementType segmentType = segment.getType();
    switch (segmentType) {
 
    case ElementType::VERTEX:
    {
        evalProjectionOnVertex(point, segmentItr, projectionPoint, projectionNormal);
        return;
    }
 
    case ElementType::LINE:
    {
        evalHighOrderProjectionOnLine(point, segmentItr, projectionPoint, projectionNormal);
        return;
    }
 
    case ElementType::TRIANGLE:
    {
        evalHighOrderProjectionOnTriangle(point, segmentItr, projectionPoint, projectionNormal);
        return;
    }
 
    default:
    {
        evalHighOrderProjectionOnPolygon(point, segmentItr, projectionPoint, projectionNormal);
        return;
    }
 
    }
}

void LevelSetSegmentationSurfaceInfo::evalLowOrderProjectionOnLine(const std::array<double, 3> &point,
                                                                   const SegmentConstIterator &segmentItr,
                                                                   std::array<double, 3> *projectionPoint,
                                                                   std::array<double, 3> *projectionNormal) const
{
    const Cell &segment = *segmentItr;
    assert(segment.getType() == ElementType::LINE);
 
    ConstProxyVector<long> segmentVertexIds = segment.getVertexIds();
 
    const std::array<double,3> &point0 = m_surface->getVertexCoords(segmentVertexIds[0]);
    const std::array<double,3> &point1 = m_surface->getVertexCoords(segmentVertexIds[1]);
 
    int nSegmentVertices = segment.getVertexCount();
    BITPIT_CREATE_WORKSPACE(lambda, double, nSegmentVertices, ReferenceElementInfo::MAX_ELEM_VERTICES);
    (*projectionPoint) = CGElem::projectPointSegment(point, point0, point1, lambda);
 
    std::array<double, 3> normal0 = computeSegmentVertexNormal(segmentItr, 0, true);
    std::array<double, 3> normal1 = computeSegmentVertexNormal(segmentItr, 1, true);
 
    (*projectionNormal) = lambda[0] * normal0 + lambda[1] * normal1;
    (*projectionNormal) /= norm2((*projectionNormal));
}

void LevelSetSegmentationSurfaceInfo::evalLowOrderProjectionOnTriangle(const std::array<double, 3> &point,
                                                                       const SegmentConstIterator &segmentItr,
                                                                       std::array<double, 3> *projectionPoint,
                                                                       std::array<double, 3> *projectionNormal) const
{
    const Cell &segment = *segmentItr;
    assert(segment.getType() == ElementType::TRIANGLE);
 
    ConstProxyVector<long> segmentVertexIds = segment.getVertexIds();
 
    const std::array<double,3> &point0 = m_surface->getVertexCoords(segmentVertexIds[0]);
    const std::array<double,3> &point1 = m_surface->getVertexCoords(segmentVertexIds[1]);
    const std::array<double,3> &point2 = m_surface->getVertexCoords(segmentVertexIds[2]);
 
    int nSegmentVertices = segment.getVertexCount();
    BITPIT_CREATE_WORKSPACE(lambda, double, nSegmentVertices, ReferenceElementInfo::MAX_ELEM_VERTICES);
    (*projectionPoint) = CGElem::projectPointTriangle(point, point0, point1, point2, lambda);
 
    std::array<double, 3> normal0 = computeSegmentVertexNormal(segmentItr, 0, true);
    std::array<double, 3> normal1 = computeSegmentVertexNormal(segmentItr, 1, true);
    std::array<double, 3> normal2 = computeSegmentVertexNormal(segmentItr, 2, true);
 
    (*projectionNormal) = lambda[0] * normal0 + lambda[1] * normal1 + lambda[2] * normal2;
    (*projectionNormal) /= norm2((*projectionNormal));
}

void LevelSetSegmentationSurfaceInfo::evalLowOrderProjectionOnPolygon(const std::array<double, 3> &point,
                                                                      const SegmentConstIterator &segmentItr,
                                                                      std::array<double, 3> *projectionPoint,
                                                                      std::array<double, 3> *projectionNormal) const
{
    const Cell &segment = *segmentItr;
 
    ConstProxyVector<long> segmentVertexIds = m_surface->getFacetOrderedVertexIds(segment);
 
    std::size_t nSegmentVertices = segmentVertexIds.size();
    BITPIT_CREATE_WORKSPACE(segmentVertexCoords, std::array<double BITPIT_COMMA 3>, nSegmentVertices, ReferenceElementInfo::MAX_ELEM_VERTICES);
    m_surface->getVertexCoords(segmentVertexIds.size(), segmentVertexIds.data(), segmentVertexCoords);
 
    BITPIT_CREATE_WORKSPACE(lambda, double, nSegmentVertices, ReferenceElementInfo::MAX_ELEM_VERTICES);
    (*projectionPoint) = CGElem::projectPointPolygon(point, nSegmentVertices, segmentVertexCoords, lambda);
 
    (*projectionNormal) = lambda[0] * computeSegmentVertexNormal(segmentItr, 0, true);
    for (std::size_t i = 1; i < nSegmentVertices; ++i) {
        (*projectionNormal) += lambda[i] * computeSegmentVertexNormal(segmentItr, i, true);
    }
    (*projectionNormal) /= norm2(*projectionNormal);
}

void LevelSetSegmentationSurfaceInfo::evalLowOrderProjection(const std::array<double, 3> &point,
                                                             const SegmentConstIterator &segmentItr,
                                                             std::array<double, 3> *projectionPoint,
                                                             std::array<double, 3> *projectionNormal) const
{
    const Cell &segment = *segmentItr;
    ElementType segmentType = segment.getType();
    switch (segmentType) {
 
    case ElementType::VERTEX:
    {
        evalProjectionOnVertex(point, segmentItr, projectionPoint, projectionNormal);
        return;
    }
 
    case ElementType::LINE:
    {
        evalLowOrderProjectionOnLine(point, segmentItr, projectionPoint, projectionNormal);
        return;
    }
 
    case ElementType::TRIANGLE:
    {
        evalLowOrderProjectionOnTriangle(point, segmentItr, projectionPoint, projectionNormal);
        return;
    }
 
    default:
    {
        evalLowOrderProjectionOnPolygon(point, segmentItr, projectionPoint, projectionNormal);
        return;
    }
 
    }
}

void LevelSetSegmentationSurfaceInfo::evalProjection(const std::array<double, 3> &point,
                                                     const SegmentConstIterator &segmentItr,
                                                     std::array<double, 3> *projectionPoint,
                                                     std::array<double, 3> *projectionNormal) const
{
    if (m_surfaceSmoothing == LevelSetSurfaceSmoothing::HIGH_ORDER) {
        evalHighOrderProjection(point, segmentItr, projectionPoint, projectionNormal);
    } else {
        evalLowOrderProjection(point, segmentItr, projectionPoint, projectionNormal);
    }
}

void LevelSetSegmentationSurfaceInfo::evalProjection(const std::array<double, 3> &point,
                                                     const SegmentConstIterator &segmentItr,
                                                     std::array<double, 3> *projectionPoint) const
{
    std::array<double, 3> projectionNormal;
    evalProjection(point, segmentItr, projectionPoint, &projectionNormal);
 
    BITPIT_UNUSED(projectionNormal);
}

std::array<double, 3> LevelSetSegmentationSurfaceInfo::evalProjection(const std::array<double, 3> &point,
                                                                      const SegmentConstIterator &segmentItr,
                                                                      double *lambda) const
{
    std::array<double, 3> projectionPoint;
    evalProjection(point, segmentItr, &projectionPoint);
 
    const Cell &segment = *segmentItr;
    m_surface->evalBarycentricCoordinates(segment.getId(), point, lambda);
    return projectionPoint;
}

std::array<double,3> LevelSetSegmentationSurfaceInfo::computePseudoNormal(const SegmentConstIterator &segmentItr,
                                                                          const double *lambda ) const
{
    // Early return if the segment is a point
    const Cell &segment = *segmentItr;
    ElementType segmentType = segment.getType();
    if (segmentType == ElementType::VERTEX) {
        return {{0., 0., 0.}};
    }
 
    // Evaluate pseudo normal
    int positionFlag;
    if (segmentType == ElementType::LINE) {
        positionFlag = CGElem::convertBarycentricToFlagSegment(lambda, m_surface->getTol());
    } else {
        int nSegmentVertices = segment.getVertexCount();
        positionFlag = CGElem::convertBarycentricToFlagPolygon(nSegmentVertices, lambda, m_surface->getTol());
        if (positionFlag > 0) {
            int polygonVertex = positionFlag - 1;
            int elementVertex = m_surface->getFacetOrderedLocalVertex(segment, polygonVertex);
 
            positionFlag = elementVertex + 1;
        } else if (positionFlag < 0) {
            int polygonEdge = (- positionFlag) - 1;
            int elementEdge = m_surface->getFacetOrderedLocalEdge(segment, polygonEdge);
 
            positionFlag = - (elementEdge + 1);
        }
    }
 
    std::array<double,3> pseudoNormal;
    if (positionFlag == 0) {
        pseudoNormal = computeSegmentNormal(segmentItr);
    } else if (positionFlag > 0) {
        int vertex = positionFlag - 1;
        pseudoNormal = computeSegmentVertexNormal(segmentItr, vertex, false);
    } else {
        int edge = (- positionFlag) - 1;
        pseudoNormal = computeSegmentEdgeNormal(segmentItr, edge, false);
    }
 
    return pseudoNormal;
}

std::array<double,3> LevelSetSegmentationSurfaceInfo::computeSurfaceNormal(const SegmentConstIterator &segmentItr,
                                                                           const double *lambda ) const
{
    // Early return if the segment is a point
    const Cell &segment = *segmentItr;
    ElementType segmentType = segment.getType();
    if (segmentType == ElementType::VERTEX) {
        return {{0., 0., 0.}};
    }
 
    // Evaluate surface normal
    std::size_t nSegmentVertices = segment.getVertexCount();
    std::array<double,3> surfaceNormal = lambda[0] * computeSegmentVertexNormal(segmentItr, 0, true);
    for (std::size_t i = 1; i < nSegmentVertices; ++i) {
        surfaceNormal += lambda[i] * computeSegmentVertexNormal(segmentItr, i, true);
    }
    surfaceNormal /= norm2(surfaceNormal);
 
    return surfaceNormal;
}

std::array<double,3> LevelSetSegmentationSurfaceInfo::computeSegmentNormal(const SegmentConstIterator &segmentItr ) const {
 
    std::size_t segmentRawId = segmentItr.getRawIndex();
    std::array<double, 3> *segmentNormal = m_segmentNormalsStorage.rawData(segmentRawId);
    if (!m_segmentNormalsValid.rawAt(segmentRawId)) {
        *segmentNormal = m_surface->evalFacetNormal(segmentItr->getId());
        m_segmentNormalsValid.rawAt(segmentRawId) = true;
    }
 
    return *segmentNormal;
}

std::array<double,3> LevelSetSegmentationSurfaceInfo::computeSegmentEdgeNormal(const SegmentConstIterator &segmentItr, int edge, bool limited ) const {
 
    long segmentId = segmentItr.getId();
 
    std::array<double, 3> limitedEdgeNormal;
    std::array<double, 3> unlimitedEdgeNormal;
    m_surface->evalEdgeNormals(segmentId, edge, m_featureAngle, &unlimitedEdgeNormal, &limitedEdgeNormal) ;
 
    if (limited) {
        return limitedEdgeNormal;
    }
 
    return unlimitedEdgeNormal;
}

std::array<double,3> LevelSetSegmentationSurfaceInfo::computeSegmentVertexNormal(const SegmentConstIterator &segmentItr, int vertex, bool limited ) const {
 
    // Segment information
    long segmentId = segmentItr.getId();
    long segmentRawId = segmentItr.getRawIndex();
    const Cell &segment = *segmentItr;
 
    // Update the cache
    //
    // Cache is updated for all the vertices of the segment.
    long vertexId = segment.getVertexId(vertex);
    std::size_t vertexRawId = m_surface->getVertices().getRawIndex(vertexId);
    bool hasUnlimitedNormal = m_unlimitedVertexNormalsValid.rawAt(vertexRawId);
 
    bool hasLimitedNormal = m_limitedSegmentVertexNormalValid[m_segmentVertexOffset.rawAt(segmentRawId) + vertex];
 
    if (!hasUnlimitedNormal || !hasLimitedNormal) {
        static std::vector<long> vertexNeighbours;
        vertexNeighbours.clear();
        m_surface->findCellVertexNeighs(segmentId, vertex, &vertexNeighbours);
 
        std::array<double, 3> limitedVertexNormal;
        std::array<double, 3> unlimitedVertexNormal;
        if (hasUnlimitedNormal) {
            limitedVertexNormal   = m_surface->evalLimitedVertexNormal(segmentId, vertex, vertexNeighbours.size(), vertexNeighbours.data(), m_featureAngle) ;
            unlimitedVertexNormal = m_unlimitedVertexNormalsStorage.rawAt(vertexRawId);
        } else {
            m_surface->evalVertexNormals(segmentId, vertex, vertexNeighbours.size(), vertexNeighbours.data(), m_featureAngle, &unlimitedVertexNormal, &limitedVertexNormal) ;
        }
 
        // Store vertex limited normal
        //
        // Both limited and unlimited normal are evaluated, however limited
        // normal is only stored if its misalignment with respect to the
        // unlimited normal is greater than a defined tolerance.
        if (!hasLimitedNormal) {
            double misalignment = norm2(unlimitedVertexNormal - limitedVertexNormal) ;
            if (misalignment >= m_surface->getTol()) {
                std::pair<long, int> segmentVertexKey = std::make_pair(segmentId, vertex);
                m_limitedSegmentVertexNormalStorage.insert({segmentVertexKey, std::move(limitedVertexNormal)}) ;
            }
            m_limitedSegmentVertexNormalValid[m_segmentVertexOffset.rawAt(segmentRawId) + vertex] = true;
        }
 
        // Store vertex unlimited normal
        if (!hasUnlimitedNormal ) {
            m_unlimitedVertexNormalsStorage.rawAt(vertexRawId) = std::move(unlimitedVertexNormal);
            m_unlimitedVertexNormalsValid.rawAt(vertexRawId) = true ;
        }
    }
 
    // Get the normal from the cache
    if (limited) {
        std::pair<long, int> segmentVertexKey = std::make_pair(segmentId, vertex);
        auto limitedSegmentVertexNormal = m_limitedSegmentVertexNormalStorage.find(segmentVertexKey);
        if (limitedSegmentVertexNormal != m_limitedSegmentVertexNormalStorage.end()) {
            return limitedSegmentVertexNormal->second;
        }
    }
 
    return m_unlimitedVertexNormalsStorage.rawAt(vertexRawId);
}


const double LevelSetSegmentationBaseObject::AUTOMATIC_SEARCH_RADIUS = -1;

std::size_t LevelSetSegmentationBaseObject::createFieldCellCache(LevelSetField field, std::size_t cacheId)
{
    switch(field) {
 
    case LevelSetField::SUPPORT:
        return createFieldCellCache<long>(field, cacheId);
 
    default:
        return LevelSetObject::createFieldCellCache(field, cacheId);
 
    }
}

LevelSetFieldset LevelSetSegmentationBaseObject::getSupportedFields() const
{
    LevelSetFieldset supportedFields = LevelSetObject::getSupportedFields();
    supportedFields.push_back(LevelSetField::PART);
    supportedFields.push_back(LevelSetField::NORMAL);
    supportedFields.push_back(LevelSetField::SUPPORT);
 
    return supportedFields;
}

void LevelSetSegmentationBaseObject::fillFieldCellCache(LevelSetField field, long id)
{
    switch (field) {
 
    case LevelSetField::SUPPORT:
        evalCellSupport(id);
        break;
 
    default:
        LevelSetObject::fillFieldCellCache(field, id);
 
    }
}

const SurfUnstructured & LevelSetSegmentationBaseObject::evalCellSurface(long id) const
{
    return _evalCellSurface(id);
}

int LevelSetSegmentationBaseObject::evalCellPart(long id) const
{
    auto evaluator = [this] (long id)
        {
            return _evalCellPart(id);
        };
 
    auto fallback = [] (long id)
        {
            BITPIT_UNUSED(id);
 
            return levelSetDefaults::PART;
        };
 
    LevelSetField field = LevelSetField::PART;
    int part = evalCellFieldCached<int>(field, id, evaluator, fallback);
 
    return part;
}

LevelSetIntersectionStatus LevelSetSegmentationBaseObject::_isCellIntersected(long id, double distance, LevelSetIntersectionMode mode) const
{
    // Get surface information
    const SurfUnstructured &surface = evalCellSurface(id);
 
    // Early return if the surface is empty
    if (surface.empty()) {
        return LevelSetIntersectionStatus::FALSE;
    }
 
    // Evaluate intersection using base class
    return LevelSetObject::_isCellIntersected(id, distance, mode);
}

LevelSetIntersectionStatus LevelSetSegmentationBaseObject::_isInterfaceIntersected(long id, bool invert,
                                                                                   std::array<std::array<double, 3>, 2> *intersection,
                                                                                   std::vector<std::array<double, 3>> *polygon) const
{
 
    // Get the tolerance used for distance comparisons
    double tolerance = getKernel()->getDistanceTolerance();
 
    // Evaluate projection of interface centroid on surface
    const Interface &interface = m_kernel->getMesh()->getInterface(id);
    std::array<double,3> interfaceCentroid = m_kernel->getMesh()->evalElementCentroid(interface);
 
    std::array<double,3> projectionPoint;
    std::array<double,3> projectionNormal;
    evalProjection(interfaceCentroid, true, &projectionPoint, &projectionNormal);
    if (!invert) {
        projectionNormal *= -1.;
    }
 
    // Early return if projection point is not the closest point of the
    // plane (defined by the projection point and projection normal) to
    // the interface centroid
    long support = evalSupport(interfaceCentroid);
    const SurfUnstructured &surface = evalSurface(interfaceCentroid);
    LevelSetSegmentationSurfaceInfo::SegmentConstIterator segmentItr = surface.getCellConstIterator(support);
 
    const Cell &segment = *segmentItr;
    ElementType segmentType = segment.getType();
    if ((segmentType == ElementType::VERTEX)
    || (m_kernel->getMesh()->getDimension() == 3 && segmentType == ElementType::LINE)) {
        return LevelSetIntersectionStatus::FALSE;
    }
 
    // Detect the intersection by imposing a Newton algorithm
    std::array<double, 3> intersectionCentroid = interfaceCentroid;
    std::array<double, 3> intersectionCentroidNew;
 
    double residual    = tolerance + 1.0;
    int iter           = 0;
    int iterMax        = 5;
    bool isIntersected = true;
    while (residual > tolerance && iter < iterMax && isIntersected) {
        ++ iter;
 
        // Eval intersection between arbitrary 2D surface and planar interface
        isIntersected = m_kernel->getMesh()->intersectInterfacePlane(id, projectionPoint, projectionNormal, intersection, polygon);
 
        intersectionCentroidNew = 0.5 * ((*intersection)[0] + (*intersection)[1]);
 
        evalProjection(intersectionCentroidNew, true, &projectionPoint, &projectionNormal);
        if (!invert) {
            projectionNormal *= -1.;
        }
 
        residual             = norm2(intersectionCentroid - intersectionCentroidNew);
        intersectionCentroid = intersectionCentroidNew;
    }
 
    // Eval intersection corresponding to the latest projection point and normal
    if (isIntersected) {
        isIntersected = m_kernel->getMesh()->intersectInterfacePlane(id, projectionPoint, projectionNormal, intersection, polygon);
        if (isIntersected) {
            return LevelSetIntersectionStatus::TRUE;
        }
    }
    return LevelSetIntersectionStatus::FALSE;
}

std::array<double,3> LevelSetSegmentationBaseObject::evalCellNormal(long id, bool signedLevelSet) const
{
    // Evaluate signed normal
    //
    // The normal stored in the cache is unsigned.
    auto evaluator = [this] (long id) -> std::array<double,3>
        {
            return _evalCellNormal(id, false);
        };
 
    auto fallback = [] (long id) -> const std::array<double,3> &
        {
            BITPIT_UNUSED(id);
 
            return levelSetDefaults::NORMAL;
        };
 
    LevelSetField field = LevelSetField::NORMAL;
    std::array<double, 3> normal = evalCellFieldCached<std::array<double, 3>>(field, id, evaluator, fallback);
 
    // Evaluate the normal with the correct signdness
    //
    // If an unsigned evaluation is requested, the orientation of the surface should be discarded
    // and in order to have a normal that is agnostic with respect the two sides of the surface.
    if (signedLevelSet) {
        short cellSign = evalCellSign(id);
        if (cellSign <= 0) {
            normal *= static_cast<double>(cellSign);;
        }
    }
 
    return normal;
}

long LevelSetSegmentationBaseObject::evalCellSupport(long id, double searchRadius) const
{
    // Evaluate support
    //
    // If we are using a limited search radius, it is not possible to use the support stored
    // in the cache, because it may be outside the search radius.
    auto evaluator = [this, searchRadius] (long id)
        {
            // Evaluate the search radius for support evaluation
            //
            // For the automatic evaluation of the search range, it is possible to use the
            // information about the cell zone:
            //  - cells that are inside the narrow band because their distance from the surface is
            //    less than the narrow band size can use a search radius equal to the narrow band
            //    size
            //  - cells that are inside the narrow band because they intersects the surface can use
            //    a search radius equal to their bounding radius;
            //  - cells that are inside the narrow band because their one of their face neighbours
            //    intersect the surface can use a search radius equal to their bounding radius plus
            //    the bounding diameter of the smallest neighbour that intersects the surface.
            double supportSearchRadius;
            if (searchRadius == AUTOMATIC_SEARCH_RADIUS) {
                LevelSetCellLocation cellLocation = getCellLocation(id);
                if (cellLocation == LevelSetCellLocation::NARROW_BAND_DISTANCE) {
                    supportSearchRadius = this->m_narrowBandSize;
                } else if (cellLocation == LevelSetCellLocation::NARROW_BAND_INTERSECTED) {
                    supportSearchRadius = m_kernel->computeCellBoundingRadius(id);
                } else if (cellLocation == LevelSetCellLocation::NARROW_BAND_NEIGHBOUR) {
                    // Evaluate the bounding diameter of the smallest intersected face neighbour
                    supportSearchRadius = std::numeric_limits<double>::max();
                    auto neighProcessor = [this, &supportSearchRadius](long neighId, int layer) {
                        BITPIT_UNUSED(layer);
 
                        // Discard neighbours that are not intersected
                        LevelSetCellLocation neighLocation = getCellLocation(neighId);
                        if (neighLocation != LevelSetCellLocation::NARROW_BAND_INTERSECTED) {
                            return false;
                        }
 
                        // Consider the bounding diameter of the smallest intersected neighbour
                        supportSearchRadius = std::min(2 * m_kernel->computeCellBoundingRadius(neighId), supportSearchRadius);
 
                        // Continue processing the other neighbours
                        return false;
                    };
 
                    const VolumeKernel &mesh = *(m_kernel->getMesh()) ;
                    mesh.processCellFaceNeighbours(id, 1, neighProcessor);
 
                    // Ad the bounding radius of the cell
                    supportSearchRadius += m_kernel->computeCellBoundingRadius(id);
                } else {
                    supportSearchRadius = std::numeric_limits<double>::max();
                }
            } else {
                supportSearchRadius = searchRadius;
            }
 
            // Evaluate cell support
            return _evalCellSupport(id, supportSearchRadius);
        };
 
    auto fallback = [] (long id)
        {
            BITPIT_UNUSED(id);
 
            return levelSetDefaults::SUPPORT;
        };
 
    LevelSetField field = LevelSetField::SUPPORT;
 
    long support;
    if (searchRadius < std::numeric_limits<double>::max() && searchRadius != AUTOMATIC_SEARCH_RADIUS) {
        // Evaluate the support from scratch
        support = evalCellField<long>(field, id, evaluator, fallback);
 
        // Update the cache if the support is valid
        if (support >= 0) {
            fillFieldCellCache(field, id, support);
        }
    } else {
        // Evaluate the support
        support = evalCellFieldCached<long>(field, id, evaluator, fallback);
    }
 
    return support;
}

const SurfUnstructured & LevelSetSegmentationBaseObject::evalSurface(const std::array<double,3> &point) const
{
    return _evalSurface(point);
}

int LevelSetSegmentationBaseObject::evalPart(const std::array<double,3> &point) const
{
    return _evalPart(point);
}

std::array<double,3> LevelSetSegmentationBaseObject::evalNormal(const std::array<double,3> &point, bool signedLevelSet) const
{
    // Project the point on the surface and evaluate the point-projection vector
    std::array<double, 3> projectionPoint;
    std::array<double, 3> projectionNormal;
    evalProjection(point, signedLevelSet, &projectionPoint, &projectionNormal);
 
    return projectionNormal;
}

long LevelSetSegmentationBaseObject::evalSupport(const std::array<double,3> &point) const
{
    return _evalSupport(point);
}

long LevelSetSegmentationBaseObject::evalSupport(const std::array<double,3> &point, double searchRadius) const
{
    return _evalSupport(point, searchRadius);
}

void LevelSetSegmentationBaseObject::evalProjection(const std::array<double,3> &point,
                                                    bool signedLevelSet,
                                                    std::array<double, 3> *projectionPoint,
                                                    std::array<double, 3> *projectionNormal) const
{
    _evalProjection(point, signedLevelSet, projectionPoint, projectionNormal);
}

int LevelSetSegmentationBaseObject::_evalCellPart(long id) const
{
    long support = evalCellSupport(id);
    const SurfUnstructured &surface = evalCellSurface(id);
    int part = surface.getCell(support).getPID();
 
    return part;
}

int LevelSetSegmentationBaseObject::_evalPart(const std::array<double,3> &point) const
{
    long support = evalSupport(point);
    const SurfUnstructured &surface = evalSurface(point);
    int part = surface.getCell(support).getPID();
 
    return part;
}

void LevelSetSegmentationBaseObject::addVTKOutputData( LevelSetField field, const std::string &objectName)
{
    VTK &vtkWriter = this->m_kernel->getMesh()->getVTK() ;
    std::string name = this->getVTKOutputDataName(field, objectName);
 
    switch (field) {
 
        case LevelSetField::SUPPORT:
            vtkWriter.addData<long>( name, VTKFieldType::SCALAR, VTKLocation::CELL, this);
            break;
 
        case LevelSetField::PART:
            vtkWriter.addData<int>( name, VTKFieldType::SCALAR, VTKLocation::CELL, this);
            break;
 
        case LevelSetField::NORMAL:
            vtkWriter.addData<double>( name, VTKFieldType::VECTOR, VTKLocation::CELL, this);
            break;
 
        default:
            LevelSetObject::addVTKOutputData(field, objectName);
            break;
 
    }
}

std::string LevelSetSegmentationBaseObject::getVTKOutputFieldName( LevelSetField field) const
{
    switch (field) {
 
        case LevelSetField::SUPPORT:
            return "SupportId";
 
        case LevelSetField::PART:
            return "PartId";
 
        case LevelSetField::NORMAL:
            return "Normal";
 
        default:
            return LevelSetObject::getVTKOutputFieldName(field);
 
    }
}

void LevelSetSegmentationBaseObject::flushVTKOutputData(std::fstream &stream, VTKFormat format,
                                                        LevelSetField field) const
{
    switch (field) {
 
    case LevelSetField::SUPPORT:
    {
        auto evaluator = [this] (long id) -> long { return evalCellSupport(id); };
        auto fallback = [] (long id) -> long { BITPIT_UNUSED(id); return levelSetDefaults::SUPPORT; };
        flushVTKOutputData<double>(stream, format, field, evaluator, fallback);
        break;
    }
 
    case LevelSetField::PART:
    {
        auto evaluator = [this] (long id) -> int { return evalCellPart(id); };
        auto fallback = [] (long id) -> int { BITPIT_UNUSED(id); return levelSetDefaults::PART; };
        flushVTKOutputData<double>(stream, format, field, evaluator, fallback);
        break;
    }
 
    case LevelSetField::NORMAL:
    {
        auto evaluator = [this] (long id) -> std::array<double,3> { return evalCellNormal(id, true); };
        auto fallback = [] (long id) -> const std::array<double,3> & { BITPIT_UNUSED(id); return levelSetDefaults::NORMAL; };
        flushVTKOutputData<double>(stream, format, field, evaluator, fallback);
        break;
    }
 
    default:
        LevelSetObject::flushVTKOutputData(stream, format, field);
        break;
 
    }
}

int LevelSetSegmentationBaseObject::getPart(long cellId) const
{
    return evalCellPart(cellId);
}

std::array<double,3> LevelSetSegmentationBaseObject::getNormal(long cellId) const
{
    return evalCellNormal(cellId, m_defaultSignedLevelSet);
}

long LevelSetSegmentationBaseObject::getSupport(long cellId) const
{
    return evalCellSupport(cellId);
}

double LevelSetSegmentationBaseObject::getSurfaceFeatureSize(long cellId) const
{
    long support = evalCellSupport(cellId);
    if (support < 0) {
        return (- levelSetDefaults::SIZE);
    }
 
    const SurfUnstructured &surface = evalCellSurface(cellId);
 
    return surface.evalCellSize(support);
}

int LevelSetSegmentationBaseObject::getPart(const std::array<double,3> &point) const
{
    return evalPart(point);
}

std::array<double,3> LevelSetSegmentationBaseObject::getNormal(const std::array<double,3> &point) const
{
    return evalNormal(point, m_defaultSignedLevelSet);
}

long LevelSetSegmentationBaseObject::getSupport(const std::array<double,3> &point) const
{
    return evalSupport(point);
}

double LevelSetSegmentationBaseObject::getSurfaceFeatureSize(const std::array<double,3> &point) const
{
    long support = evalSupport(point);
    if (support < 0) {
        return (- levelSetDefaults::SIZE);
    }
 
    const SurfUnstructured &surface = evalSurface(point);
 
    return surface.evalCellSize(support);
}


LevelSetSegmentationObject::LevelSetSegmentationObject(int id)
    : LevelSetSegmentationBaseObject(id),
      m_surfaceInfo(nullptr)
{
}

LevelSetSegmentationObject::LevelSetSegmentationObject(int id, std::unique_ptr<const SurfUnstructured> &&surface, double featureAngle)
    : LevelSetSegmentationBaseObject(id)
{
    setSurface(std::move(surface), featureAngle);
}

LevelSetSegmentationObject::LevelSetSegmentationObject(int id, const SurfUnstructured *surface, double featureAngle)
    : LevelSetSegmentationBaseObject(id)
{
    setSurface(surface, featureAngle);
}

LevelSetSegmentationObject::LevelSetSegmentationObject(int id, const SurfUnstructured *surface, double featureAngle, LevelSetSurfaceSmoothing surfaceSmoothing)
    : LevelSetSegmentationBaseObject(id)
{
    setSurface(surface, featureAngle, surfaceSmoothing);
}

LevelSetSegmentationObject::LevelSetSegmentationObject(const LevelSetSegmentationObject &other)
    : LevelSetSegmentationBaseObject(other)
{
    if (other.m_surfaceInfo) {
        m_surfaceInfo = std::unique_ptr<LevelSetSegmentationSurfaceInfo>(new LevelSetSegmentationSurfaceInfo(*(other.m_surfaceInfo)));
    } else {
        m_surfaceInfo = nullptr;
    }
}

bool LevelSetSegmentationObject::empty() const
{
    return getSurface().empty();
}

LevelSetSegmentationObject * LevelSetSegmentationObject::clone() const {
    return new LevelSetSegmentationObject(*this );
}

const SurfUnstructured & LevelSetSegmentationObject::getSurface() const {
    return m_surfaceInfo->getSurface();
}

void LevelSetSegmentationObject::setSurface(std::unique_ptr<const SurfUnstructured> &&surface, bool force){
    setSurface(std::move(surface), LevelSetSegmentationSurfaceInfo::DEFAULT_FEATURE_ANGLE, force);
}

void LevelSetSegmentationObject::setSurface(std::unique_ptr<const SurfUnstructured> &&surface, double featureAngle, bool force){
    if (m_surfaceInfo) {
        // Check if replacing an existing surface is allowed
        if (!force) {
            throw std::runtime_error ("The surface can only be set once.");
        }
 
        // Replace the surface
        m_surfaceInfo->setSurface(std::move(surface), featureAngle);
    } else {
        // Set surface
        //
        // Since this is the first time we set the surface, there is no need
        // to clear the caches.
        m_surfaceInfo = std::unique_ptr<LevelSetSegmentationSurfaceInfo>(new LevelSetSegmentationSurfaceInfo(std::move(surface), featureAngle));
    }
}

void LevelSetSegmentationObject::setSurface(const SurfUnstructured *surface, bool force){
    setSurface(surface, LevelSetSegmentationSurfaceInfo::DEFAULT_FEATURE_ANGLE, force);
}

void LevelSetSegmentationObject::setSurface(const SurfUnstructured *surface, double featureAngle, bool force){
    if (m_surfaceInfo) {
        // Check if replacing an existing surface is allowed
        if (!force) {
            throw std::runtime_error ("The surface can only be set once.");
        }
 
        // Replace the surface
        m_surfaceInfo->setSurface(surface, featureAngle);
    } else {
        // Set surface
        //
        // Since this is the first time we set the surface, there is no need
        // to clear the caches.
        m_surfaceInfo = std::unique_ptr<LevelSetSegmentationSurfaceInfo>(new LevelSetSegmentationSurfaceInfo(surface, featureAngle));
    }
}

void LevelSetSegmentationObject::setSurface(const SurfUnstructured *surface, double featureAngle, LevelSetSurfaceSmoothing surfaceSmoothing, bool force){
    if (m_surfaceInfo) {
        // Check if replacing an existing surface is allowed
        if (!force) {
            throw std::runtime_error ("The surface can only be set once.");
        }
 
        // Replace the surface
        m_surfaceInfo->setSurface(surface, featureAngle, surfaceSmoothing);
    } else {
        // Set surface
        //
        // Since this is the first time we set the surface, there is no need
        // to clear the caches.
        m_surfaceInfo = std::unique_ptr<LevelSetSegmentationSurfaceInfo>(new LevelSetSegmentationSurfaceInfo(surface, featureAngle, surfaceSmoothing));
    }
}

const SurfaceSkdTree & LevelSetSegmentationObject::getSearchTree() const {
    return m_surfaceInfo->getSearchTree();
}

double LevelSetSegmentationObject::getFeatureAngle() const {
    return m_surfaceInfo->getFeatureAngle();
}

bitpit::LevelSetSurfaceSmoothing LevelSetSegmentationObject::getSurfaceSmoothing() const {
    return m_surfaceInfo->getSurfaceSmoothing();
}

void LevelSetSegmentationObject::fillCellLocationCache()
{
    // Cartesian patches are handled separately
    if (dynamic_cast<LevelSetCartesianKernel*>(m_kernel)) {
        fillCartesianCellZoneCache();
        return;
    }
 
    // All other patches are handled with the base method.
    LevelSetObject::fillCellLocationCache();
}

void LevelSetSegmentationObject::fillCellLocationCache(const std::vector<adaption::Info> &adaptionData)
{
    // Cartesian patches are handled separately
    //
    // Update is not implemented for Cartesian patches, the cells that should be inserted in
    // the cache will be evaluated considering all the mash not just the elements newly added.
    if (dynamic_cast<LevelSetCartesianKernel*>(m_kernel)) {
        fillCartesianCellZoneCache();
        return;
    }
 
    // All other patches are handled with the base method
    LevelSetObject::fillCellLocationCache(adaptionData);
}

void LevelSetSegmentationObject::fillCartesianCellZoneCache()
{
    // The function needs a Cartesian kernel
    const LevelSetCartesianKernel *cartesianKernel = dynamic_cast<LevelSetCartesianKernel*>(m_kernel);
    if (!dynamic_cast<LevelSetCartesianKernel*>(m_kernel)) {
        throw std::runtime_error("The function needs a Cartesian kernels.");
    }
 
    // Get mesh information
    const VolCartesian &mesh = *(cartesianKernel->getMesh() ) ;
    int meshDimension = mesh.getDimension();
    long nCells = mesh.getCellCount();
 
    std::array<double,3> meshMinPoint;
    std::array<double,3> meshMaxPoint;
    mesh.getBoundingBox(meshMinPoint, meshMaxPoint) ;
 
    // Get surface information
    const SurfUnstructured &surface = getSurface();
 
    // Initialize process list
    //
    // Process list is initialized with cells that are certainly inside the
    // narrow band. Those cells are the ones that contain the vertices of the
    // segments or the intersection between the segments and the bounding box
    // of the patch.
    std::unordered_set<long> processList;
 
    std::vector<std::array<double,3>> intersectionPoints;
    std::vector<std::array<double,3>> segmentVertexCoords;
    for (const Cell &segment : surface.getCells()) {
        // Get segment info
        //
        // Since vertex information will be passed to the CG module, we need to get the
        // vertices in the proper order.
        ConstProxyVector<long> segmentVertexIds = surface.getFacetOrderedVertexIds(segment);
        std::size_t nSegmentVertices = segmentVertexIds.size();
 
        // Get segment coordinates
        segmentVertexCoords.resize(nSegmentVertices);
        surface.getVertexCoords(nSegmentVertices, segmentVertexIds.data(), segmentVertexCoords.data());
 
        // Add to the process list the cells that contain the vertices of the
        // segment or the intersection between the segment and the bounding box
        // of the patch.
        int nInnerVertices = 0;
        for (const std::array<double,3> &vertexPoint : segmentVertexCoords) {
            long cellId = mesh.locatePoint(vertexPoint);
            if (cellId < 0) {
                continue;
            }
 
            processList.insert(cellId);
            ++nInnerVertices;
        }
 
        if (nInnerVertices == 0) {
            if (CGElem::intersectBoxPolygon(meshMinPoint, meshMaxPoint, segmentVertexCoords, false, true, true, intersectionPoints, meshDimension)) {
                for (const std::array<double,3> &intersectionPoint : intersectionPoints){
                    long cellId = mesh.locateClosestCell(intersectionPoint);
                    processList.insert(cellId);
                }
            }
        }
    }
 
    // Get cache for zone identification
    CellCacheCollection::ValueCache<char> *locationCache = getCellCache<char>(m_cellLocationCacheId);
 
    // Cells with an unknown region are in the bulk
    for (long cellId = 0; cellId < nCells; ++cellId) {
        locationCache->insertEntry(cellId, static_cast<char>(LevelSetCellLocation::UNKNOWN));
    }
 
    // Start filling the list of cells within the narrow band
    //
    // The initial process list is gradually expanded considering all the neighbours
    // inside the narrow band.
    std::vector<long> intersectedCellIds;
    std::vector<bool> alreadyProcessed(nCells, false);
    while (!processList.empty()) {
        // Get the cell to process
        long cellId = *(processList.begin());
        processList.erase(processList.begin());
 
        // Skip cells already processed
        if (alreadyProcessed[cellId]) {
            continue;
        }
        alreadyProcessed[cellId] = true;
 
        // Fill location cache for cells geometrically inside the narrow band
        LevelSetCellLocation cellLocation = fillCellGeometricNarrowBandLocationCache(cellId);
 
        // Track intersected cells
        if (cellLocation == LevelSetCellLocation::NARROW_BAND_INTERSECTED) {
            intersectedCellIds.push_back(cellId);
        }
 
        // Add unprocessed neighbours of cells inside the narrow band to the process list
        if (cellLocation != LevelSetCellLocation::UNKNOWN) {
            auto neighProcessor = [&processList, &alreadyProcessed](long neighId, int layer) {
                BITPIT_UNUSED(layer);
 
                // Skip neighbours already processed
                if (alreadyProcessed[neighId]) {
                    return false;
                }
 
                // Update the process list
                processList.insert(neighId);
 
                // Continue processing the other neighbours
                return false;
            };
 
            mesh.processCellFaceNeighbours(cellId, 1, neighProcessor);
        }
    }
 
    // Identify neighbours of cells that intersect the surface
    for (std::size_t cellId : intersectedCellIds) {
        // Process face neighbours
        auto neighProcessor = [this, &locationCache](long neighId, int layer) {
            BITPIT_UNUSED(layer);
 
            // Skip neighbours whose region has already been identified
            if (getCellLocation(neighId) != LevelSetCellLocation::UNKNOWN) {
                return false;
            }
 
            // The neighbour is inside the narrow band
            locationCache->insertEntry(neighId, static_cast<char>(LevelSetCellLocation::NARROW_BAND_NEIGHBOUR));
 
            // Continue processing the other neighbours
            return false;
        };
 
        mesh.processCellFaceNeighbours(cellId, 1, neighProcessor);
    }
 
    // Cells with an unknown region are in the bulk
    for (long cellId = 0; cellId < nCells; ++cellId) {
        CellCacheCollection::ValueCache<char>::Entry locationCacheEntry = locationCache->findEntry(cellId);
        if (*locationCacheEntry == static_cast<char>(LevelSetCellLocation::UNKNOWN)) {
            locationCache->insertEntry(cellId, static_cast<char>(LevelSetCellLocation::BULK));
        }
    }
 
#if BITPIT_ENABLE_MPI==1
    // Exchange ghost data
    if (mesh.isPartitioned()) {
        std::unique_ptr<DataCommunicator> dataCommunicator = m_kernel->createDataCommunicator();
        startCellCacheExchange(mesh.getGhostCellExchangeSources(), m_cellLocationCacheId, dataCommunicator.get());
        completeCellCacheExchange(mesh.getGhostCellExchangeTargets(), m_cellLocationCacheId, dataCommunicator.get());
    }
#endif
}

LevelSetCellLocation LevelSetSegmentationObject::fillCellGeometricNarrowBandLocationCache(long id)
{
    // Early return if the cell is geometrically outside the narrow band
    double searchRadius = std::max(m_kernel->computeCellBoundingRadius(id), m_narrowBandSize);
    long cellSupport = evalCellSupport(id, searchRadius);
    if (cellSupport < 0) {
        return LevelSetCellLocation::UNKNOWN;
    }
 
    // Evaluate levelset value
    std::array<double,3> cellCentroid = m_kernel->computeCellCentroid(id);
 
    double cellCacheValue    = _evalValue(cellCentroid, cellSupport, CELL_CACHE_IS_SIGNED);
    double cellUnsigendValue = std::abs(cellCacheValue);
 
    // Update the cell location cache
    //
    // First we need to check if the cell intersects the surface, and only if it
    // doesn't we  should check if its distance is lower than the narrow band size.
    //
    // When high order smoothing is active, there may be cells whose support is within the search
    // radius, but they are not intersected and their distance is less than the narrow band size.
    // These cells are not geometrically inside the narrow band, they are neighbours of cells
    // geometrically inside the narrow band and as such it's up to the caller of this function to
    // identify their cell location.
    LevelSetCellLocation cellLocation = LevelSetCellLocation::UNKNOWN;
    if (_isCellIntersected(id, cellUnsigendValue, CELL_LOCATION_INTERSECTION_MODE) == LevelSetIntersectionStatus::TRUE) {
        cellLocation = LevelSetCellLocation::NARROW_BAND_INTERSECTED;
    } else if (cellUnsigendValue <= m_narrowBandSize) {
        cellLocation = LevelSetCellLocation::NARROW_BAND_DISTANCE;
    }
    assert((getSurfaceSmoothing() == LevelSetSurfaceSmoothing::HIGH_ORDER) || (cellLocation != LevelSetCellLocation::UNKNOWN));
 
    if (cellLocation != LevelSetCellLocation::UNKNOWN) {
        CellCacheCollection::ValueCache<char> *locationCache = getCellCache<char>(m_cellLocationCacheId);
        locationCache->insertEntry(id, static_cast<char>(cellLocation));
    }
 
    // Fill the cache of the evaluated fields
    //
    // Now that the cell location has been identified, we can fill the cache of the
    // evaluated fields.
    fillFieldCellCache(LevelSetField::SUPPORT, id, cellSupport);
    fillFieldCellCache(LevelSetField::VALUE, id, cellCacheValue);
 
    // Return the location
    return cellLocation;
}

const SurfUnstructured & LevelSetSegmentationObject::_evalCellSurface(long id) const
{
    BITPIT_UNUSED(id);
 
    return getSurface();
}

short LevelSetSegmentationObject::_evalCellSign(long id) const
{
    long support = evalCellSupport(id);
    std::array<double,3> centroid = m_kernel->computeCellCentroid(id);
 
    return _evalSign(centroid, support);
}

double LevelSetSegmentationObject::_evalCellValue(long id, bool signedLevelSet) const
{
    long support = evalCellSupport(id);
    std::array<double,3> centroid = m_kernel->computeCellCentroid(id);
 
    return _evalValue(centroid, support, signedLevelSet);
}

std::array<double,3> LevelSetSegmentationObject::_evalCellGradient(long id, bool signedLevelSet) const
{
    long support = evalCellSupport(id);
    std::array<double,3> centroid = m_kernel->computeCellCentroid(id);
 
    return _evalGradient(centroid, support, signedLevelSet);
}

long LevelSetSegmentationObject::_evalCellSupport(long id, double searchRadius) const
{
    std::array<double,3> centroid = m_kernel->computeCellCentroid(id);
 
    return _evalSupport(centroid, searchRadius);
}

std::array<double,3> LevelSetSegmentationObject::_evalCellNormal(long id, bool signedLevelSet) const
{
    long support = evalCellSupport(id);
    std::array<double,3> centroid = m_kernel->computeCellCentroid(id);
 
    std::array<double, 3> projectionPoint;
    std::array<double, 3> projectionNormal;
    _evalProjection(centroid, support, signedLevelSet, &projectionPoint, &projectionNormal);
 
    return projectionNormal;
}

const SurfUnstructured & LevelSetSegmentationObject::_evalSurface(const std::array<double,3> &point) const
{
    BITPIT_UNUSED(point);
 
    return getSurface();
}

short LevelSetSegmentationObject::_evalSign(const std::array<double,3> &point) const
{
    long support = evalSupport(point);
 
    return _evalSign(point, support);
}

double LevelSetSegmentationObject::_evalValue(const std::array<double,3> &point, bool signedLevelSet) const
{
    long support = evalSupport(point);
 
    return _evalValue(point, support, signedLevelSet);
}

std::array<double,3> LevelSetSegmentationObject::_evalGradient(const std::array<double,3> &point, bool signedLevelSet) const
{
    long support = evalSupport(point);
 
    return _evalGradient(point, support, signedLevelSet);
}

void LevelSetSegmentationObject::_evalProjection(const std::array<double,3> &point,
                                                 bool signedLevelSet,
                                                 std::array<double, 3> *projectionPoint,
                                                 std::array<double, 3> *projectionNormal) const
{
    // Get closest segment
    long support = evalSupport(point);
 
    _evalProjection(point, support, signedLevelSet, projectionPoint, projectionNormal);
}

short LevelSetSegmentationObject::_evalSign(const std::array<double,3> &point, long support) const
{
    // Throw an error if the support is not valid
    if (support < 0) {
        if (empty()) {
            return levelSetDefaults::SIGN;
        }
 
        throw std::runtime_error("Unable to evaluate the sign: the support is not valid.");
    }
 
    LevelSetSegmentationSurfaceInfo::SegmentConstIterator supportItr = getSurface().getCellConstIterator(support);
 
    return evalValueSign(m_surfaceInfo->evalDistance(point, supportItr, true));
}

double LevelSetSegmentationObject::_evalValue(const std::array<double,3> &point, long support,
                                              bool signedLevelSet) const
{
    // Early return if the support is not valid
    //
    // With an invalid support, only the unsigend levelset can be evaluated.
    if (support < 0) {
        if (!signedLevelSet || empty()) {
            return levelSetDefaults::VALUE;
        }
 
        throw std::runtime_error("With an invalid support, only the unsigend levelset can be evaluated.");
    }
 
    // Evaluate the distance of the point from the surface
    LevelSetSegmentationSurfaceInfo::SegmentConstIterator supportItr = getSurface().getCellConstIterator(support);
    double distance = m_surfaceInfo->evalDistance(point, supportItr, signedLevelSet);
 
    // Early return if the point lies on the surface
    if (evalValueSign(distance) == 0) {
        return 0.;
    }
 
    // Evaluate levelset value
    double value;
    if (signedLevelSet) {
        value = distance;
    } else {
        value = std::abs(distance);
    }
 
    return value;
}

std::array<double,3> LevelSetSegmentationObject::_evalGradient(const std::array<double,3> &point, long support,
                                                               bool signedLevelSet) const
{
    // Early return if the support is not valid
    //
    // With an invalid support, only the unsigend levelset can be evaluated.
    if (support < 0) {
        if (!signedLevelSet || empty()) {
            return levelSetDefaults::GRADIENT;
        }
 
        throw std::runtime_error("With an invalid support, only the unsigend levelset can be evaluated.");
    }
 
    // Evaluate the distance of the point from the surface
    LevelSetSegmentationSurfaceInfo::SegmentConstIterator supportItr = getSurface().getCellConstIterator(support);
    std::array<double,3> distanceVector;
    double distance = m_surfaceInfo->evalDistance(point, supportItr, signedLevelSet, &distanceVector);
 
    // Early return if the point lies on the surface
    if (evalValueSign(distance) == 0) {
        if (signedLevelSet) {
            return m_surfaceInfo->evalNormal(point, supportItr);
        } else {
            return {{0., 0., 0.}};
        }
    }
 
    // Evaluate levelset gradient
    std::array<double,3> gradient = distanceVector / distance;
 
    return gradient;
}

long LevelSetSegmentationObject::_evalSupport(const std::array<double,3> &point) const
{
    return _evalSupport(point, std::numeric_limits<double>::max());
}

long LevelSetSegmentationObject::_evalSupport(const std::array<double,3> &point, double searchRadius) const
{
    long closestSegmentId;
    double closestDistance;
    m_surfaceInfo->getSearchTree().findPointClosestCell(point, searchRadius, &closestSegmentId, &closestDistance);
 
    return closestSegmentId;
}

void LevelSetSegmentationObject::_evalProjection(const std::array<double,3> &point,
                                                 long support,
                                                 bool signedLevelSet,
                                                 std::array<double, 3> *projectionPoint,
                                                 std::array<double, 3> *projectionNormal) const
{
    // Early return if projection normal and projection point pointers are null
    if (!(projectionNormal) && !(projectionPoint)) {
       return;
    }
 
    // Early return if the support is not valid
    //
    // With an invalid support, only the unsigend levelset can be evaluated.
    if (support < 0) {
        if (!signedLevelSet || empty()) {
            (*projectionPoint) = levelSetDefaults::POINT;
            (*projectionNormal) = levelSetDefaults::NORMAL;
        }
 
        throw std::runtime_error("With an invalid support, only the unsigend levelset can be evaluated.");
    }
 
    // Get closest segment
    LevelSetSegmentationSurfaceInfo::SegmentConstIterator segmentItr = getSurface().getCellConstIterator(support);
 
    // Eval projection point and normal
    m_surfaceInfo->evalProjection(point, segmentItr, projectionPoint, projectionNormal);
 
    // Early return if a projection normal null pointer is given
    if (!(projectionNormal)) {
       return;
    }
 
    // If an unsigned evaluation is requested, the orientation of the surface should be discarded
    // in order to have a normal that is agnostic with respect the two sides of the surface.
    if (!signedLevelSet) {
        short sign = evalSign(point);
        (*projectionNormal) *= static_cast<double>(sign);
    }
}

LevelSetIntersectionStatus LevelSetSegmentationObject::_isInterfaceIntersected(long id, bool invert,
                                                                               std::array<std::array<double, 3>, 2> *intersection,
                                                                               std::vector<std::array<double, 3>> *polygon) const
{
    // Early return if low order smoothing is chosen
    if (m_surfaceInfo->getSurfaceSmoothing() == LevelSetSurfaceSmoothing::LOW_ORDER) {
        return LevelSetObject::_isInterfaceIntersected(id, invert, intersection, polygon);
    }
    return LevelSetSegmentationBaseObject::_isInterfaceIntersected(id, invert, intersection, polygon);
}

double LevelSetSegmentationObject::getMinSurfaceFeatureSize() const {
 
    const SurfUnstructured &surface = getSurface();
 
    bool   minimumValid = false;
    double minimumSize  = levelSetDefaults::SIZE;
    for (const Cell &cell : surface.getCells()) {
        double segmentSize = surface.evalCellSize(cell.getId());
        if (segmentSize < 0) {
            continue;
        }
 
        minimumValid = true;
        minimumSize  = std::min(segmentSize, minimumSize);
    }
 
    if (!minimumValid) {
        minimumSize = - levelSetDefaults::SIZE;
    }
 
    return minimumSize;
}

double LevelSetSegmentationObject::getMaxSurfaceFeatureSize() const {
 
    const SurfUnstructured &surface = getSurface();
 
    double maximumSize = - levelSetDefaults::SIZE;
    for (const Cell &cell : surface.getCells()) {
        double segmentSize = surface.evalCellSize(cell.getId());
        maximumSize = std::max(segmentSize, maximumSize);
    }
 
    return maximumSize;
}

LevelSetBooleanObject<LevelSetSegmentationBaseObject>::LevelSetBooleanObject( int id, LevelSetBooleanOperation op, const LevelSetSegmentationBaseObject *source1, const LevelSetSegmentationBaseObject *source2  )
    : LevelSetBooleanBaseObject<LevelSetSegmentationBaseObject>(id, op, source1, source2) {
}

LevelSetBooleanObject<LevelSetSegmentationBaseObject>::LevelSetBooleanObject( int id, LevelSetBooleanOperation op, const std::vector<const LevelSetSegmentationBaseObject *> &sourceObjects )
    : LevelSetBooleanBaseObject<LevelSetSegmentationBaseObject>(id, op, sourceObjects) {
}

LevelSetBooleanObject<LevelSetSegmentationBaseObject> * LevelSetBooleanObject<LevelSetSegmentationBaseObject>::clone() const {
    return new LevelSetBooleanObject<LevelSetSegmentationBaseObject>(*this );
}

const SurfUnstructured & LevelSetBooleanObject<LevelSetSegmentationBaseObject>::_evalCellSurface(long id) const
{
    return getCellReferenceObject(id)->evalCellSurface(id);
}

long LevelSetBooleanObject<LevelSetSegmentationBaseObject>::_evalCellSupport(long id, double searchRadius) const
{
    return getCellReferenceObject(id)->evalCellSupport(id, searchRadius);
}

int LevelSetBooleanObject<LevelSetSegmentationBaseObject>::_evalCellPart(long id) const
{
    return getCellReferenceObject(id)->evalCellPart(id);
}

std::array<double,3> LevelSetBooleanObject<LevelSetSegmentationBaseObject>::_evalCellNormal(long id, bool signedLevelSet) const
{
    return _evalCellFunction<std::array<double,3>>(id, signedLevelSet, [&id, signedLevelSet] (const LevelSetBooleanResult<LevelSetSegmentationBaseObject> &result) -> std::array<double,3>
        {
            const LevelSetSegmentationBaseObject *resultObject = result.getObject();
            if ( !resultObject ) {
                return levelSetDefaults::NORMAL;
            }
 
            std::array<double,3> normal = resultObject->evalCellNormal(id, signedLevelSet);
            if (signedLevelSet) {
                normal *= static_cast<double>(result.getObjectSign());
            }
 
            return normal;
        });
}

const SurfUnstructured & LevelSetBooleanObject<LevelSetSegmentationBaseObject>::_evalSurface(const std::array<double,3> &point) const
{
    return getReferenceObject(point)->evalSurface(point);
}

long LevelSetBooleanObject<LevelSetSegmentationBaseObject>::_evalSupport(const std::array<double,3> &point) const
{
    return getReferenceObject(point)->evalSupport(point);
}

long LevelSetBooleanObject<LevelSetSegmentationBaseObject>::_evalSupport(const std::array<double,3> &point, double searchRadius) const
{
    return getReferenceObject(point)->evalSupport(point, searchRadius);
}

void LevelSetBooleanObject<LevelSetSegmentationBaseObject>::_evalProjection(const std::array<double,3> &point,
                                                                            bool signedLevelSet,
                                                                            std::array<double, 3> *projectionPoint,
                                                                            std::array<double, 3> *projectionNormal) const
{
    return _evalFunction<void>(point, signedLevelSet, [&point, projectionPoint, projectionNormal, signedLevelSet] (const LevelSetBooleanResult<LevelSetSegmentationBaseObject> &result)
        {
            const LevelSetSegmentationBaseObject *resultObject = result.getObject();
            if ( !resultObject ) {
                (*projectionNormal) = levelSetDefaults::NORMAL;
                (*projectionPoint)  = levelSetDefaults::POINT;
            }
 
            resultObject->evalProjection(point, signedLevelSet, projectionPoint, projectionNormal);
            if (signedLevelSet) {
                (*projectionNormal) *= static_cast<double>(result.getObjectSign());
            }
        });
}

int LevelSetBooleanObject<LevelSetSegmentationBaseObject>::_evalPart(const std::array<double,3> &point) const
{
    return getReferenceObject(point)->evalPart(point);
}

LevelSetComplementObject<LevelSetSegmentationBaseObject>::LevelSetComplementObject(int id, const LevelSetSegmentationBaseObject *source)
    : LevelSetComplementBaseObject<LevelSetSegmentationBaseObject>(id, source)
{
}

LevelSetComplementObject<LevelSetSegmentationBaseObject> * LevelSetComplementObject<LevelSetSegmentationBaseObject>::clone() const {
    return new LevelSetComplementObject<LevelSetSegmentationBaseObject>(*this );
}

const SurfUnstructured & LevelSetComplementObject<LevelSetSegmentationBaseObject>::_evalCellSurface(long id) const
{
    return getSourceObject()->evalCellSurface(id);
}

long LevelSetComplementObject<LevelSetSegmentationBaseObject>::_evalCellSupport(long id, double searchRadius) const
{
    return getSourceObject()->evalCellSupport(id, searchRadius);
}

int LevelSetComplementObject<LevelSetSegmentationBaseObject>::_evalCellPart(long id) const
{
    return getCellReferenceObject(id)->evalCellPart(id);
}

std::array<double,3> LevelSetComplementObject<LevelSetSegmentationBaseObject>::_evalCellNormal(long id, bool signedLevelSet) const
{
    std::array<double,3> normal = getSourceObject()->evalCellNormal(id, signedLevelSet);
    if (signedLevelSet) {
        normal *= -1.;
    }
 
    return normal;
}

const SurfUnstructured & LevelSetComplementObject<LevelSetSegmentationBaseObject>::_evalSurface(const std::array<double,3> &point) const
{
    return getSourceObject()->evalSurface(point);
}

long LevelSetComplementObject<LevelSetSegmentationBaseObject>::_evalSupport(const std::array<double,3> &point) const
{
    return getSourceObject()->evalSupport(point);
}

long LevelSetComplementObject<LevelSetSegmentationBaseObject>::_evalSupport(const std::array<double,3> &point, double searchRadius) const
{
    return getSourceObject()->evalSupport(point, searchRadius);
}

int LevelSetComplementObject<LevelSetSegmentationBaseObject>::_evalPart(const std::array<double,3> &point) const
{
    return getSourceObject()->evalPart(point);
}

void LevelSetComplementObject<LevelSetSegmentationBaseObject>::_evalProjection(const std::array<double,3> &point,
                                                                               bool signedLevelSet,
                                                                               std::array<double, 3> *projectionPoint,
                                                                               std::array<double, 3> *projectionNormal) const
{
    getSourceObject()->evalProjection(point, signedLevelSet, projectionPoint, projectionNormal);
    if (signedLevelSet) {
        (*projectionNormal) *= -1.;
    }
}
 
}