FreeCAD / src /Mod /CAM /App /VoronoiEdgePyImp.cpp
Upload folder using huggingface_hub
985c397 verified 9 months ago
26.4 kB
// SPDX-License-Identifier: LGPL-2.1-or-later
/***************************************************************************
 *   Copyright (c) 2020 sliptonic <shopinthewoods@gmail.com>               *
 *                                                                         *
 *   This file is part of the FreeCAD CAx development system.              *
 *                                                                         *
 *   This library is free software; you can redistribute it and/or         *
 *   modify it under the terms of the GNU Library General Public           *
 *   License as published by the Free Software Foundation; either          *
 *   version 2 of the License, or (at your option) any later version.      *
 *                                                                         *
 *   This library  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 Library General Public License for more details.                  *
 *                                                                         *
 *   You should have received a copy of the GNU Library General Public     *
 *   License along with this library; see the file COPYING.LIB. If not,    *
 *   write to the Free Software Foundation, Inc., 59 Temple Place,         *
 *   Suite 330, Boston, MA  02111-1307, USA                                *
 *                                                                         *
 ***************************************************************************/

#include <limits>

#include <BRepBuilderAPI_MakeEdge.hxx>
#include <Geom_Parabola.hxx>


#include "Mod/Part/App/Geometry.h"
#include "Mod/Part/App/TopoShapeEdgePy.h"

#include "VoronoiEdgePy.h"
#include "VoronoiEdgePy.cpp"
#include "VoronoiCellPy.h"
#include "VoronoiVertexPy.h"


using namespace Path;

namespace
{

Voronoi::point_type pointFromVertex(const Voronoi::vertex_type v)
{
    Voronoi::point_type pt;
    pt.x(v.x());
    pt.y(v.y());
    return pt;
}

Voronoi::point_type orthognalProjection(
    const Voronoi::point_type& point,
    const Voronoi::segment_type& segment
)
{
    // move segment so it goes through the origin (s)
    Voronoi::point_type offset;
    {
        offset.x(low(segment).x());
        offset.y(low(segment).y());
    }
    Voronoi::point_type s;
    {
        s.x(high(segment).x() - offset.x());
        s.y(high(segment).y() - offset.y());
    }
    // move point accordingly so it maintains it's relation to s (p)
    Voronoi::point_type p;
    {
        p.x(point.x() - offset.x());
        p.y(point.y() - offset.y());
    }
    // calculate the orthogonal projection of p onto s
    // ((p dot s) / (s dot s)) * s
    // (https://en.wikibooks.org/wiki/Linear_Algebra/Orthogonal_Projection_Onto_a_Line) and it back
    // by original offset to get the projected point
    const double proj = (p.x() * s.x() + p.y() * s.y())
        / (s.x() * s.x() + s.y() * s.y() + std::numeric_limits<double>::epsilon());
    Voronoi::point_type pt;
    {
        pt.x(offset.x() + proj * s.x());
        pt.y(offset.y() + proj * s.y());
    }
    return pt;
}

double length(const Voronoi::point_type& p)
{
    return sqrt(p.x() * p.x() + p.y() * p.y());
}

int sideOf(const Voronoi::point_type& p, const Voronoi::segment_type& s)
{
    Voronoi::coordinate_type dxp = p.x() - low(s).x();
    Voronoi::coordinate_type dyp = p.y() - low(s).y();
    Voronoi::coordinate_type dxs = high(s).x() - low(s).x();
    Voronoi::coordinate_type dys = high(s).y() - low(s).y();

    double d = -dxs * dyp + dys * dxp;
    if (d < 0) {
        return -1;
    }
    if (d > 0) {
        return +1;
    }
    return 0;
}

template<typename pt0_type, typename pt1_type>
double distanceBetween(const pt0_type& p0, const pt1_type& p1, double scale)
{
    Voronoi::point_type dist;
    dist.x(p0.x() - p1.x());
    dist.y(p0.y() - p1.y());
    return length(dist) / scale;
}

template<typename pt0_type, typename pt1_type>
double signedDistanceBetween(const pt0_type& p0, const pt1_type& p1, double scale)
{
    if (length(p0) > length(p1)) {
        return -distanceBetween(p0, p1, scale);
    }
    return distanceBetween(p0, p1, scale);
}


void addDistanceBetween(
    const Voronoi::diagram_type::vertex_type* v0,
    const Voronoi::point_type& p1,
    Py::List* list,
    double scale
)
{
    if (v0) {
        list->append(Py::Float(distanceBetween(*v0, p1, scale)));
    }
    else {
        Py_INCREF(Py_None);
        list->append(Py::asObject(Py_None));
    }
}

void addProjectedDistanceBetween(
    const Voronoi::diagram_type::vertex_type* v0,
    const Voronoi::segment_type& segment,
    Py::List* list,
    double scale
)
{
    if (v0) {
        Voronoi::point_type p0;
        {
            p0.x(v0->x());
            p0.y(v0->y());
        }
        Voronoi::point_type p1 = orthognalProjection(p0, segment);
        list->append(Py::Float(distanceBetween(*v0, p1, scale)));
    }
    else {
        Py_INCREF(Py_None);
        list->append(Py::asObject(Py_None));
    }
}

bool addDistancesToPoint(const VoronoiEdge* edge, Voronoi::point_type p, Py::List* list, double scale)
{
    addDistanceBetween(edge->ptr->vertex0(), p, list, scale);
    addDistanceBetween(edge->ptr->vertex1(), p, list, scale);
    return true;
}

bool retrieveDistances(const VoronoiEdge* edge, Py::List* list)
{
    const Voronoi::diagram_type::cell_type* c0 = edge->ptr->cell();
    if (c0->contains_point()) {
        return addDistancesToPoint(edge, edge->dia->retrievePoint(c0), list, edge->dia->getScale());
    }
    const Voronoi::diagram_type::cell_type* c1 = edge->ptr->twin()->cell();
    if (c1->contains_point()) {
        return addDistancesToPoint(edge, edge->dia->retrievePoint(c1), list, edge->dia->getScale());
    }
    // at this point both cells are sourced from segments and it does not matter which one we use
    Voronoi::segment_type segment = edge->dia->retrieveSegment(c0);
    addProjectedDistanceBetween(edge->ptr->vertex0(), segment, list, edge->dia->getScale());
    addProjectedDistanceBetween(edge->ptr->vertex1(), segment, list, edge->dia->getScale());
    return false;
}

bool pointsMatch(const Voronoi::point_type& p0, const Voronoi::point_type& p1, double scale)
{
    return 1e-6 > distanceBetween(p0, p1, scale);
}

bool isPointOnSegment(const Voronoi::point_type& point, const Voronoi::segment_type& segment, double scale)
{
    return pointsMatch(point, low(segment), scale) || pointsMatch(point, high(segment), scale);
}

template<typename T>
PyObject* makeLineSegment(const VoronoiEdge* e, const T& p0, double z0, const T& p1, double z1)
{
    Part::GeomLineSegment p;
    p.setPoints(e->dia->scaledVector(p0, z0), e->dia->scaledVector(p1, z1));
    Handle(Geom_Curve) h = Handle(Geom_Curve)::DownCast(p.handle());
    BRepBuilderAPI_MakeEdge mkBuilder(h, h->FirstParameter(), h->LastParameter());
    return new Part::TopoShapeEdgePy(new Part::TopoShape(mkBuilder.Shape()));
}
}  // namespace

std::ostream& operator<<(std::ostream& os, const Voronoi::vertex_type& v)
{
    return os << '(' << v.x() << ", " << v.y() << ')';
}

std::ostream& operator<<(std::ostream& os, const Voronoi::point_type& p)
{
    return os << '(' << p.x() << ", " << p.y() << ')';
}

std::ostream& operator<<(std::ostream& os, const Voronoi::segment_type& s)
{
    return os << '<' << low(s) << ", " << high(s) << '>';
}


// returns a string which represents the object e.g. when printed in python
std::string VoronoiEdgePy::representation() const
{
    std::stringstream ss;
    ss.precision(5);
    ss << "VoronoiEdge(";
    VoronoiEdge* e = getVoronoiEdgePtr();
    if (e->isBound()) {
        const Voronoi::diagram_type::vertex_type* v0 = e->ptr->vertex0();
        const Voronoi::diagram_type::vertex_type* v1 = e->ptr->vertex1();
        if (v0) {
            ss << "[" << (v0->x() / e->dia->getScale()) << ", " << (v0->y() / e->dia->getScale())
               << "]";
        }
        else {
            ss << "[~]";
        }
        ss << ", ";
        if (v1) {
            ss << "[" << (v1->x() / e->dia->getScale()) << ", " << (v1->y() / e->dia->getScale())
               << "]";
        }
        else {
            ss << "[~]";
        }
    }
    ss << ")";
    return ss.str();
}

PyObject* VoronoiEdgePy::PyMake(struct _typeobject*, PyObject*, PyObject*)  // Python wrapper
{
    // create a new instance of VoronoiEdgePy and the Twin object
    return new VoronoiEdgePy(new VoronoiEdge);
}

// constructor method
int VoronoiEdgePy::PyInit(PyObject* args, PyObject* /*kwd*/)
{
    if (!PyArg_ParseTuple(args, "")) {
        PyErr_SetString(PyExc_RuntimeError, "no arguments accepted");
        return -1;
    }
    return 0;
}


PyObject* VoronoiEdgePy::richCompare(PyObject* lhs, PyObject* rhs, int op)
{
    PyObject* cmp = (op == Py_EQ) ? Py_False : Py_True;
    if (PyObject_TypeCheck(lhs, &VoronoiEdgePy::Type)
        && PyObject_TypeCheck(rhs, &VoronoiEdgePy::Type) && (op == Py_EQ || op == Py_NE)) {
        const VoronoiEdge* vl = static_cast<VoronoiEdgePy*>(lhs)->getVoronoiEdgePtr();
        const VoronoiEdge* vr = static_cast<VoronoiEdgePy*>(rhs)->getVoronoiEdgePtr();
        if (vl->dia == vr->dia && vl->index == vr->index) {
            cmp = (op == Py_EQ) ? Py_True : Py_False;
        }
    }
    Py_INCREF(cmp);
    return cmp;
}

const Voronoi::voronoi_diagram_type::edge_type* getEdgeFromPy(
    VoronoiEdgePy* e,
    bool throwIfNotBound = true
)
{
    auto self = e->getVoronoiEdgePtr();
    if (self->isBound()) {
        return self->ptr;
    }
    if (throwIfNotBound) {
        throw Py::TypeError("Edge not bound to voronoi diagram");
    }
    return nullptr;
}

VoronoiEdge* getVoronoiEdgeFromPy(const VoronoiEdgePy* e, PyObject* args = nullptr)
{
    VoronoiEdge* self = e->getVoronoiEdgePtr();
    if (!self->isBound()) {
        throw Py::TypeError("Edge not bound to voronoi diagram");
    }
    if (args && !PyArg_ParseTuple(args, "")) {
        throw Py::RuntimeError("No arguments accepted");
    }
    return self;
}

Py::Long VoronoiEdgePy::getIndex() const
{
    VoronoiEdge* e = getVoronoiEdgePtr();
    if (e->isBound()) {
        return Py::Long(e->dia->index(e->ptr));
    }
    return Py::Long(-1);
}

Py::Long VoronoiEdgePy::getColor() const
{
    VoronoiEdge* e = getVoronoiEdgePtr();
    if (e->isBound()) {
        Voronoi::color_type color = e->ptr->color() & Voronoi::ColorMask;
        return Py::Long(PyLong_FromSize_t(color));
    }
    return Py::Long(0);
}

void VoronoiEdgePy::setColor(Py::Long color)
{
    getEdgeFromPy(this)->color(long(color) & Voronoi::ColorMask);
}

Py::List VoronoiEdgePy::getVertices() const
{
    Py::List list;
    VoronoiEdge* e = getVoronoiEdgePtr();
    if (e->isBound()) {
        auto v0 = e->ptr->vertex0();
        auto v1 = e->ptr->vertex1();
        if (v0) {
            list.append(Py::asObject(new VoronoiVertexPy(new VoronoiVertex(e->dia, v0))));
        }
        else {
            Py_INCREF(Py_None);
            list.append(Py::asObject(Py_None));
        }
        if (v1) {
            list.append(Py::asObject(new VoronoiVertexPy(new VoronoiVertex(e->dia, v1))));
        }
        else {
            Py_INCREF(Py_None);
            list.append(Py::asObject(Py_None));
        }
    }
    return list;
}

Py::Object VoronoiEdgePy::getTwin() const
{
    VoronoiEdge* e = getVoronoiEdgeFromPy(this);
    return Py::asObject(new VoronoiEdgePy(new VoronoiEdge(e->dia, e->ptr->twin())));
}

Py::Object VoronoiEdgePy::getNext() const
{
    VoronoiEdge* e = getVoronoiEdgeFromPy(this);
    return Py::asObject(new VoronoiEdgePy(new VoronoiEdge(e->dia, e->ptr->next())));
}

Py::Object VoronoiEdgePy::getPrev() const
{
    VoronoiEdge* e = getVoronoiEdgeFromPy(this);
    return Py::asObject(new VoronoiEdgePy(new VoronoiEdge(e->dia, e->ptr->prev())));
}

Py::Object VoronoiEdgePy::getRotNext() const
{
    VoronoiEdge* e = getVoronoiEdgeFromPy(this);
    return Py::asObject(new VoronoiEdgePy(new VoronoiEdge(e->dia, e->ptr->rot_next())));
}

Py::Object VoronoiEdgePy::getRotPrev() const
{
    VoronoiEdge* e = getVoronoiEdgeFromPy(this);
    return Py::asObject(new VoronoiEdgePy(new VoronoiEdge(e->dia, e->ptr->rot_prev())));
}

Py::Object VoronoiEdgePy::getCell() const
{
    VoronoiEdge* e = getVoronoiEdgeFromPy(this);
    return Py::asObject(new VoronoiCellPy(new VoronoiCell(e->dia, e->ptr->cell())));
}


PyObject* VoronoiEdgePy::isFinite(PyObject* args) const
{
    VoronoiEdge* e = getVoronoiEdgeFromPy(this, args);
    PyObject* chk = e->ptr->is_finite() ? Py_True : Py_False;
    Py_INCREF(chk);
    return chk;
}

PyObject* VoronoiEdgePy::isInfinite(PyObject* args) const
{
    VoronoiEdge* e = getVoronoiEdgeFromPy(this, args);
    PyObject* chk = e->ptr->is_infinite() ? Py_True : Py_False;
    Py_INCREF(chk);
    return chk;
}

PyObject* VoronoiEdgePy::isLinear(PyObject* args) const
{
    VoronoiEdge* e = getVoronoiEdgeFromPy(this, args);
    PyObject* chk = e->ptr->is_linear() ? Py_True : Py_False;
    Py_INCREF(chk);
    return chk;
}

PyObject* VoronoiEdgePy::isCurved(PyObject* args) const
{
    VoronoiEdge* e = getVoronoiEdgeFromPy(this, args);
    PyObject* chk = e->ptr->is_curved() ? Py_True : Py_False;
    Py_INCREF(chk);
    return chk;
}

PyObject* VoronoiEdgePy::isPrimary(PyObject* args) const
{
    VoronoiEdge* e = getVoronoiEdgeFromPy(this, args);
    PyObject* chk = e->ptr->is_primary() ? Py_True : Py_False;
    Py_INCREF(chk);
    return chk;
}

PyObject* VoronoiEdgePy::isSecondary(PyObject* args) const
{
    VoronoiEdge* e = getVoronoiEdgeFromPy(this, args);
    PyObject* chk = e->ptr->is_secondary() ? Py_True : Py_False;
    Py_INCREF(chk);
    return chk;
}

PyObject* VoronoiEdgePy::isBorderline(PyObject* args) const
{
    VoronoiEdge* e = getVoronoiEdgeFromPy(this, args);
    PyObject* chk = Py_False;
    if (e->isBound() && !e->ptr->is_linear()) {
        Voronoi::point_type point = e->ptr->cell()->contains_point()
            ? e->dia->retrievePoint(e->ptr->cell())
            : e->dia->retrievePoint(e->ptr->twin()->cell());
        Voronoi::segment_type segment = e->ptr->cell()->contains_point()
            ? e->dia->retrieveSegment(e->ptr->twin()->cell())
            : e->dia->retrieveSegment(e->ptr->cell());
        if (isPointOnSegment(point, segment, e->dia->getScale())) {
            chk = Py_True;
        }
    }
    Py_INCREF(chk);
    return chk;
}

PyObject* VoronoiEdgePy::toShape(PyObject* args) const
{
    double z0 = 0.0;
    double z1 = std::numeric_limits<double>::max();
    int dbg = 0;
    if (!PyArg_ParseTuple(args, "|ddp", &z0, &z1, &dbg)) {
        throw Py::RuntimeError("no, one or two arguments of type double accepted");
    }
    if (z1 == std::numeric_limits<double>::max()) {
        z1 = z0;
    }
    VoronoiEdge* e = getVoronoiEdgePtr();
    if (e->isBound()) {
        if (e->ptr->is_linear()) {
            if (e->ptr->is_finite()) {
                auto v0 = e->ptr->vertex0();
                auto v1 = e->ptr->vertex1();
                if (v0 && v1) {
                    return makeLineSegment(e, *v0, z0, *v1, z1);
                }
            }
            else {
                // infinite linear, need to clip somehow
                const Voronoi::diagram_type::cell_type* c0 = e->ptr->cell();
                const Voronoi::diagram_type::cell_type* c1 = e->ptr->twin()->cell();
                Voronoi::point_type origin;
                Voronoi::point_type direction;
                if (c0->contains_point() && c1->contains_point()) {
                    Voronoi::point_type p0 = e->dia->retrievePoint(c0);
                    Voronoi::point_type p1 = e->dia->retrievePoint(c1);
                    origin.x((p0.x() + p1.x()) / 2.);
                    origin.y((p0.y() + p1.y()) / 2.);
                    direction.x(p0.y() - p1.y());
                    direction.y(p1.x() - p0.x());
                }
                else {
                    origin = c0->contains_segment() ? e->dia->retrievePoint(c1)
                                                    : e->dia->retrievePoint(c0);
                    Voronoi::segment_type segment = c0->contains_segment()
                        ? e->dia->retrieveSegment(c0)
                        : e->dia->retrieveSegment(c1);
                    Voronoi::coordinate_type dx = high(segment).x() - low(segment).x();
                    Voronoi::coordinate_type dy = high(segment).y() - low(segment).y();
                    if ((low(segment) == origin) ^ c0->contains_point()) {
                        direction.x(dy);
                        direction.y(-dx);
                    }
                    else {
                        direction.x(-dy);
                        direction.y(dx);
                    }
                }
                double k = 2.5;  // <-- need something smarter here
                Voronoi::point_type begin;
                Voronoi::point_type end;
                if (e->ptr->vertex0()) {
                    begin.x(e->ptr->vertex0()->x());
                    begin.y(e->ptr->vertex0()->y());
                }
                else {
                    begin.x(origin.x() - direction.x() * k);
                    begin.y(origin.y() - direction.y() * k);
                }
                if (e->ptr->vertex1()) {
                    end.x(e->ptr->vertex1()->x());
                    end.y(e->ptr->vertex1()->y());
                }
                else {
                    end.x(origin.x() + direction.x() * k);
                    end.y(origin.y() + direction.y() * k);
                }
                return makeLineSegment(e, begin, z0, end, z1);
            }
        }
        else {
            // parabolic curve, which is always formed by a point and an edge
            Voronoi::point_type point = e->ptr->cell()->contains_point()
                ? e->dia->retrievePoint(e->ptr->cell())
                : e->dia->retrievePoint(e->ptr->twin()->cell());
            Voronoi::segment_type segment = e->ptr->cell()->contains_point()
                ? e->dia->retrieveSegment(e->ptr->twin()->cell())
                : e->dia->retrieveSegment(e->ptr->cell());
            // the location is the mid point between the normal on the segment through point
            // this is only the mid point of the segment if the parabola is symmetric

            if (isPointOnSegment(point, segment, e->dia->getScale())) {
                return makeLineSegment(e, low(segment), z0, high(segment), z1);
            }

            Voronoi::point_type loc;
            {
                Voronoi::point_type proj = orthognalProjection(point, segment);
                // the location is the mid point between the projection on the segment and the point
                loc.x((proj.x() + point.x()) / 2);
                loc.y((proj.y() + point.y()) / 2);
            }
            Voronoi::point_type axis;
            {
                axis.x(point.x() - loc.x());
                axis.y(point.y() - loc.y());
            }
            Voronoi::segment_type xaxis;
            {
                xaxis.low(point);
                xaxis.high(loc);
            }

            // determine distances of the end points from the x-axis, those are the parameters for
            // the arc of the parabola in the horizontal plane
            auto pt0 = pointFromVertex(*e->ptr->vertex0());
            auto pt1 = pointFromVertex(*e->ptr->vertex1());
            Voronoi::point_type pt0x = orthognalProjection(pt0, xaxis);
            Voronoi::point_type pt1x = orthognalProjection(pt1, xaxis);
            double dist0 = distanceBetween(pt0, pt0x, e->dia->getScale()) * sideOf(pt0, xaxis);
            double dist1 = distanceBetween(pt1, pt1x, e->dia->getScale()) * sideOf(pt1, xaxis);
            if (dist1 < dist0) {
                // if the parabola is traversed in the revere direction we need to use the points
                // on the other side of the parabola - 'beauty of symmetric geometries
                dist0 = -dist0;
                dist1 = -dist1;
            }

            // at this point we have the direction of the x-axis and the two end points p0 and p1
            // which means we know the plane of the parabola
            auto p0 = e->dia->scaledVector(pt0, z0);
            auto p1 = e->dia->scaledVector(pt1, z1);
            // we get a third point by moving p0 along the axis of the parabola
            auto p_ = p0 + e->dia->scaledVector(axis, 0);

            // normal of the plane defined by those 3 points
            auto norm = ((p_ - p0).Cross(p1 - p0)).Normalize();

            // the next thing to figure out is the z level of the x-axis,
            double zx = z0 - (dist0 / (dist0 - dist1)) * (z0 - z1);

            auto locn = e->dia->scaledVector(loc, zx);
            auto xdir = e->dia->scaledVector(axis, zx);

            double focal;
            if (z0 == z1) {
                // focal length if parabola in the xy-plane is simply half the distance between the
                // point and segment - aka the distance between point and location, aka the length
                // of axis
                focal = length(axis) / e->dia->getScale();
                if (dbg) {
                    std::cerr << "focal = " << length(axis) << "/" << e->dia->getScale() << "\n";
                }
            }
            else {
                // if the parabola is not in the xy-plane we need to find the
                // (x,y) coordinates of a point on the parabola in the parabola's
                // coordinate system.
                // see: http://amsi.org.au/ESA_Senior_Years/SeniorTopic2/2a/2a_2content_10.html
                // note that above website uses Y as the symmetry axis of the parabola whereas
                // OCC uses X as the symmetry axis. The math below is in the website's system.
                // We already know 2 points on the parabola (p0 and p1), we know their X values
                // (dist0 and dist1) if the parabola is in the xy-plane, and we know their
                // orthogonal projection onto the parabola's symmetry axis Y (pt0x and pt1x). The
                // resulting Y values are the distance between the parabola's location (loc) and the
                // respective orthogonal projection. Pythagoras gives us the X values, using the X
                // from the xy-plane and the difference in z. Note that this calculation also gives
                // correct results if the parabola is in the xy-plane (z0 == z1), it's just that
                // above calculation is so much simpler.
                double flenX0 = sqrt(dist0 * dist0 + (z0 - zx) * (z0 - zx));
                double flenX1 = sqrt(dist1 * dist1 + (zx - z1) * (zx - z1));
                double flenX;
                double flenY;
                // if one of the points is the location, we have to use the other to get sensible
                // values
                if (fabs(dist0) > fabs(dist1)) {
                    flenX = flenX0;
                    flenY = distanceBetween(loc, pt0x, e->dia->getScale());
                }
                else {
                    flenX = flenX1;
                    flenY = distanceBetween(loc, pt1x, e->dia->getScale());
                }
                // parabola: (x - p)^2 = 4*focal*(y - q)   |  (p,q) ... location of parabola
                focal = (flenX * flenX) / (4 * fabs(flenY));
                if (dbg) {
                    std::cerr << "segment" << segment << ", point" << point << std::endl;
                    std::cerr << "  loc" << loc << ", axis" << axis << std::endl;
                    std::cerr << "  dist0(" << dist0 << " : " << flenX0 << ", dist1(" << dist1
                              << " : " << flenX1 << ")" << std::endl;
                    std::cerr << "  z(" << z0 << ", " << zx << ", " << z1 << ")" << std::endl;
                    std::cerr << "  focal = (" << flenX << " * " << flenX << ") / (4 * fabs("
                              << flenY << "))\n";
                }
                // use new X values to set the parameters
                dist0 = dist0 >= 0 ? flenX0 : -flenX0;
                dist1 = dist1 >= 0 ? flenX1 : -flenX1;
            }

            gp_Pnt pbLocn(locn.x, locn.y, locn.z);
            gp_Dir pbNorm(norm.x, norm.y, norm.z);
            gp_Dir pbXdir(xdir.x, xdir.y, 0);

            gp_Ax2 pb(pbLocn, pbNorm, pbXdir);
            Handle(Geom_Parabola) parabola = new Geom_Parabola(pb, focal);

            Part::GeomArcOfParabola arc;
            arc.setHandle(parabola);
            arc.setRange(dist0, dist1, false);

            // get a shape for the parabola arc
            Handle(Geom_Curve) h = Handle(Geom_Curve)::DownCast(arc.handle());
            BRepBuilderAPI_MakeEdge mkBuilder(h, h->FirstParameter(), h->LastParameter());
            return new Part::TopoShapeEdgePy(new Part::TopoShape(mkBuilder.Shape()));
        }
    }
    Py_INCREF(Py_None);
    return Py_None;
}


PyObject* VoronoiEdgePy::getDistances(PyObject* args) const
{
    VoronoiEdge* e = getVoronoiEdgeFromPy(this, args);
    Py::List list;
    retrieveDistances(e, &list);
    return Py::new_reference_to(list);
}

PyObject* VoronoiEdgePy::getSegmentAngle(PyObject* args) const
{
    using std::numbers::pi;

    VoronoiEdge* e = getVoronoiEdgeFromPy(this, args);

    if (e->ptr->cell()->contains_segment() && e->ptr->twin()->cell()->contains_segment()) {
        int i0 = e->ptr->cell()->source_index() - e->dia->points.size();
        int i1 = e->ptr->twin()->cell()->source_index() - e->dia->points.size();
        if (e->dia->segmentsAreConnected(i0, i1)) {
            double a0 = e->dia->angleOfSegment(i0);
            double a1 = e->dia->angleOfSegment(i1);
            double a = a0 - a1;
            if (a > pi / 2) {
                a -= pi;
            }
            else if (a < -pi / 2) {
                a += pi;
            }
            return Py::new_reference_to(Py::Float(a));
        }
    }
    Py_INCREF(Py_None);
    return Py_None;
}

// custom attributes get/set

PyObject* VoronoiEdgePy::getCustomAttributes(const char* /*attr*/) const
{
    return nullptr;
}

int VoronoiEdgePy::setCustomAttributes(const char* /*attr*/, PyObject* /*obj*/)
{
    return 0;
}