1005 lines
38 KiB
C++
1005 lines
38 KiB
C++
// SPDX-License-Identifier: LGPL-2.1-or-later
|
|
|
|
/***************************************************************************
|
|
* Copyright (c) 2008 Juergen Riegel <juergen.riegel@web.de> *
|
|
* *
|
|
* 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 <FCConfig.h>
|
|
|
|
#ifdef FC_OS_LINUX
|
|
#include <unistd.h>
|
|
#endif
|
|
#include <BRepAdaptor_Curve.hxx>
|
|
#include <BRepBuilderAPI_MakeVertex.hxx>
|
|
#include <BRepExtrema_DistShapeShape.hxx>
|
|
#include <BRep_Tool.hxx>
|
|
#include <BndLib_Add3dCurve.hxx>
|
|
#include <Bnd_Box.hxx>
|
|
#include <GCPnts_AbscissaPoint.hxx>
|
|
#include <GCPnts_UniformAbscissa.hxx>
|
|
#include <GCPnts_UniformDeflection.hxx>
|
|
#include <GeomAPI_IntCS.hxx>
|
|
#include <Geom_Curve.hxx>
|
|
#include <Geom_Plane.hxx>
|
|
#include <Standard_Failure.hxx>
|
|
#include <TopExp_Explorer.hxx>
|
|
#include <TopoDS.hxx>
|
|
#include <TopoDS_Edge.hxx>
|
|
#include <gp_Pln.hxx>
|
|
|
|
#include <Base/Console.h>
|
|
#include <Base/FileInfo.h>
|
|
#include <Base/Sequencer.h>
|
|
#include <Base/Stream.h>
|
|
|
|
#include <Mod/Mesh/App/Core/Algorithm.h>
|
|
#include <Mod/Mesh/App/Core/Grid.h>
|
|
#include <Mod/Mesh/App/Core/Iterator.h>
|
|
#include <Mod/Mesh/App/Core/MeshKernel.h>
|
|
#include <Mod/Mesh/App/Core/Projection.h>
|
|
|
|
#include "MeshAlgos.h"
|
|
|
|
|
|
using namespace MeshPart;
|
|
using MeshCore::MeshAlgorithm;
|
|
using MeshCore::MeshFacet;
|
|
using MeshCore::MeshFacetGrid;
|
|
using MeshCore::MeshFacetIterator;
|
|
using MeshCore::MeshKernel;
|
|
using MeshCore::MeshPointIterator;
|
|
|
|
CurveProjector::CurveProjector(const TopoDS_Shape& aShape, const MeshKernel& pMesh)
|
|
: _Shape(aShape)
|
|
, _Mesh(pMesh)
|
|
{}
|
|
|
|
void CurveProjector::writeIntersectionPointsToFile(const char* name)
|
|
{
|
|
// export points
|
|
Base::FileInfo fi(name);
|
|
Base::ofstream str(fi, std::ios::out | std::ios::binary);
|
|
str.precision(4);
|
|
str.setf(std::ios::fixed | std::ios::showpoint);
|
|
for (const auto& it1 : mvEdgeSplitPoints) {
|
|
for (const auto& it2 : it1.second) {
|
|
str << it2.p1.x << " " << it2.p1.y << " " << it2.p1.z << std::endl;
|
|
}
|
|
}
|
|
str.close();
|
|
}
|
|
|
|
|
|
//**************************************************************************
|
|
//**************************************************************************
|
|
// Separator for additional classes
|
|
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|
|
|
|
CurveProjectorShape::CurveProjectorShape(const TopoDS_Shape& aShape, const MeshKernel& pMesh)
|
|
: CurveProjector(aShape, pMesh)
|
|
{
|
|
CurveProjectorShape::Do();
|
|
}
|
|
|
|
void CurveProjectorShape::Do()
|
|
{
|
|
TopExp_Explorer Ex;
|
|
for (Ex.Init(_Shape, TopAbs_EDGE); Ex.More(); Ex.Next()) {
|
|
const TopoDS_Edge& aEdge = TopoDS::Edge(Ex.Current());
|
|
projectCurve(aEdge, mvEdgeSplitPoints[aEdge]);
|
|
}
|
|
}
|
|
|
|
|
|
void CurveProjectorShape::projectCurve(const TopoDS_Edge& aEdge,
|
|
std::vector<FaceSplitEdge>& vSplitEdges)
|
|
{
|
|
Standard_Real fFirst, fLast;
|
|
Handle(Geom_Curve) hCurve = BRep_Tool::Curve(aEdge, fFirst, fLast);
|
|
|
|
// getting start point
|
|
gp_Pnt gpPt = hCurve->Value(fFirst);
|
|
|
|
// projection of the first point
|
|
Base::Vector3f cStartPoint = Base::Vector3f((float)gpPt.X(), (float)gpPt.Y(), (float)gpPt.Z());
|
|
Base::Vector3f cResultPoint, cSplitPoint, cPlanePnt, cPlaneNormal;
|
|
MeshCore::FacetIndex uStartFacetIdx, uCurFacetIdx;
|
|
MeshCore::FacetIndex uLastFacetIdx =
|
|
MeshCore::FACET_INDEX_MAX - 1; // use another value as FACET_INDEX_MAX
|
|
MeshCore::FacetIndex auNeighboursIdx[3];
|
|
bool GoOn;
|
|
|
|
if (!findStartPoint(_Mesh, cStartPoint, cResultPoint, uStartFacetIdx)) {
|
|
return;
|
|
}
|
|
|
|
uCurFacetIdx = uStartFacetIdx;
|
|
do {
|
|
MeshGeomFacet cCurFacet = _Mesh.GetFacet(uCurFacetIdx);
|
|
_Mesh.GetFacetNeighbours(uCurFacetIdx,
|
|
auNeighboursIdx[0],
|
|
auNeighboursIdx[1],
|
|
auNeighboursIdx[2]);
|
|
Base::Vector3f PointOnEdge[3];
|
|
|
|
GoOn = false;
|
|
int NbrOfHits = 0, HitIdx = 0;
|
|
|
|
for (int i = 0; i < 3; i++) {
|
|
// ignore last visited facet
|
|
if (auNeighboursIdx[i] == uLastFacetIdx) {
|
|
continue;
|
|
}
|
|
|
|
// get points of the edge i
|
|
const Base::Vector3f& cP0 = cCurFacet._aclPoints[i];
|
|
const Base::Vector3f& cP1 = cCurFacet._aclPoints[(i + 1) % 3];
|
|
|
|
if (auNeighboursIdx[i] != MeshCore::FACET_INDEX_MAX) {
|
|
// calculate the normal by the edge vector and the middle between the two face
|
|
// normals
|
|
MeshGeomFacet N = _Mesh.GetFacet(auNeighboursIdx[i]);
|
|
cPlaneNormal = (N.GetNormal() + cCurFacet.GetNormal()) % (cP1 - cP0);
|
|
cPlanePnt = cP0;
|
|
}
|
|
else {
|
|
// with no neighbours the face normal is used
|
|
cPlaneNormal = cCurFacet.GetNormal() % (cP1 - cP0);
|
|
cPlanePnt = cP0;
|
|
}
|
|
|
|
Handle(Geom_Plane) hPlane =
|
|
new Geom_Plane(gp_Pln(gp_Pnt(cPlanePnt.x, cPlanePnt.y, cPlanePnt.z),
|
|
gp_Dir(cPlaneNormal.x, cPlaneNormal.y, cPlaneNormal.z)));
|
|
|
|
GeomAPI_IntCS Alg(hCurve, hPlane);
|
|
|
|
if (Alg.IsDone()) {
|
|
// deciding by the number of result points (intersections)
|
|
if (Alg.NbPoints() == 1) {
|
|
gp_Pnt P = Alg.Point(1);
|
|
float l = ((Base::Vector3f((float)P.X(), (float)P.Y(), (float)P.Z()) - cP0)
|
|
* (cP1 - cP0))
|
|
/ ((cP1 - cP0) * (cP1 - cP0));
|
|
// is the Point on the Edge of the facet?
|
|
if (l < 0.0 || l > 1.0) {
|
|
PointOnEdge[i] = Base::Vector3f(std::numeric_limits<float>::max(), 0, 0);
|
|
}
|
|
else {
|
|
cSplitPoint = (1 - l) * cP0 + l * cP1;
|
|
PointOnEdge[i] = (1 - l) * cP0 + l * cP1;
|
|
NbrOfHits++;
|
|
HitIdx = i;
|
|
}
|
|
// no intersection
|
|
}
|
|
else if (Alg.NbPoints() == 0) {
|
|
PointOnEdge[i] = Base::Vector3f(std::numeric_limits<float>::max(), 0, 0);
|
|
// more the one intersection (@ToDo)
|
|
}
|
|
else if (Alg.NbPoints() > 1) {
|
|
PointOnEdge[i] = Base::Vector3f(std::numeric_limits<float>::max(), 0, 0);
|
|
Base::Console().log("MeshAlgos::projectCurve(): More then one intersection in "
|
|
"Facet %lu, Edge %d\n",
|
|
uCurFacetIdx,
|
|
i);
|
|
}
|
|
}
|
|
}
|
|
|
|
uLastFacetIdx = uCurFacetIdx;
|
|
|
|
if (NbrOfHits == 1) {
|
|
uCurFacetIdx = auNeighboursIdx[HitIdx];
|
|
FaceSplitEdge splitEdge;
|
|
splitEdge.ulFaceIndex = uCurFacetIdx;
|
|
splitEdge.p1 = cResultPoint;
|
|
splitEdge.p2 = cSplitPoint;
|
|
vSplitEdges.push_back(splitEdge);
|
|
cResultPoint = cSplitPoint;
|
|
GoOn = true;
|
|
}
|
|
else {
|
|
Base::Console().log("MeshAlgos::projectCurve(): Possible reentry in Facet %lu\n",
|
|
uCurFacetIdx);
|
|
}
|
|
|
|
if (uCurFacetIdx == uStartFacetIdx) {
|
|
GoOn = false;
|
|
}
|
|
|
|
} while (GoOn);
|
|
}
|
|
|
|
bool CurveProjectorShape::findStartPoint(const MeshKernel& MeshK,
|
|
const Base::Vector3f& Pnt,
|
|
Base::Vector3f& Rslt,
|
|
MeshCore::FacetIndex& FaceIndex)
|
|
{
|
|
Base::Vector3f TempResultPoint;
|
|
float MinLength = std::numeric_limits<float>::max();
|
|
bool bHit = false;
|
|
|
|
// go through the whole Mesh
|
|
MeshFacetIterator It(MeshK);
|
|
for (It.Init(); It.More(); It.Next()) {
|
|
// try to project (with angle) to the face
|
|
if (It->Foraminate(Pnt, It->GetNormal(), TempResultPoint)) {
|
|
// distance to the projected point
|
|
float Dist = (Pnt - TempResultPoint).Length();
|
|
if (Dist < MinLength) {
|
|
// remember the point with the closest distance
|
|
bHit = true;
|
|
MinLength = Dist;
|
|
Rslt = TempResultPoint;
|
|
FaceIndex = It.Position();
|
|
}
|
|
}
|
|
}
|
|
return bHit;
|
|
}
|
|
|
|
|
|
//**************************************************************************
|
|
//**************************************************************************
|
|
// Separator for CurveProjectorSimple classes
|
|
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|
|
|
|
|
|
CurveProjectorSimple::CurveProjectorSimple(const TopoDS_Shape& aShape, const MeshKernel& pMesh)
|
|
: CurveProjector(aShape, pMesh)
|
|
{
|
|
Do();
|
|
}
|
|
|
|
|
|
void CurveProjectorSimple::Do()
|
|
{
|
|
TopExp_Explorer Ex;
|
|
|
|
std::vector<Base::Vector3f> vEdgePolygon;
|
|
|
|
for (Ex.Init(_Shape, TopAbs_EDGE); Ex.More(); Ex.Next()) {
|
|
const TopoDS_Edge& aEdge = TopoDS::Edge(Ex.Current());
|
|
projectCurve(aEdge, vEdgePolygon, mvEdgeSplitPoints[aEdge]);
|
|
}
|
|
}
|
|
|
|
void CurveProjectorSimple::GetSampledCurves(const TopoDS_Edge& aEdge,
|
|
std::vector<Base::Vector3f>& rclPoints,
|
|
unsigned long ulNbOfPoints)
|
|
{
|
|
rclPoints.clear();
|
|
|
|
Standard_Real fBegin, fEnd;
|
|
|
|
Handle(Geom_Curve) hCurve = BRep_Tool::Curve(aEdge, fBegin, fEnd);
|
|
float fLen = float(fEnd - fBegin);
|
|
|
|
for (unsigned long i = 0; i < ulNbOfPoints; i++) {
|
|
gp_Pnt gpPt = hCurve->Value(fBegin + (fLen * float(i)) / float(ulNbOfPoints - 1));
|
|
rclPoints.emplace_back((float)gpPt.X(), (float)gpPt.Y(), (float)gpPt.Z());
|
|
}
|
|
}
|
|
|
|
|
|
// projectToNeighbours(Handle(Geom_Curve) hCurve,float pos
|
|
|
|
void CurveProjectorSimple::projectCurve(const TopoDS_Edge& aEdge,
|
|
const std::vector<Base::Vector3f>&,
|
|
std::vector<FaceSplitEdge>&)
|
|
{
|
|
Base::Vector3f TempResultPoint;
|
|
bool bFirst = true;
|
|
|
|
Standard_Real fBegin, fEnd;
|
|
Handle(Geom_Curve) hCurve = BRep_Tool::Curve(aEdge, fBegin, fEnd);
|
|
float fLen = float(fEnd - fBegin);
|
|
|
|
unsigned long ulNbOfPoints = 1000, PointCount = 0;
|
|
|
|
MeshFacetIterator It(_Mesh);
|
|
|
|
Base::SequencerLauncher seq("Building up projection map...", ulNbOfPoints + 1);
|
|
Base::FileInfo fi("projected.asc");
|
|
Base::ofstream str(fi, std::ios::out | std::ios::binary);
|
|
str.precision(4);
|
|
str.setf(std::ios::fixed | std::ios::showpoint);
|
|
|
|
std::map<MeshCore::FacetIndex, std::vector<Base::Vector3f>> FaceProjctMap;
|
|
|
|
for (unsigned long i = 0; i <= ulNbOfPoints; i++) {
|
|
seq.next();
|
|
gp_Pnt gpPt = hCurve->Value(fBegin + (fLen * float(i)) / float(ulNbOfPoints - 1));
|
|
|
|
// go through the whole Mesh
|
|
for (It.Init(); It.More(); It.Next()) {
|
|
// try to project (with angle) to the face
|
|
if (It->IntersectWithLine(
|
|
Base::Vector3f((float)gpPt.X(), (float)gpPt.Y(), (float)gpPt.Z()),
|
|
It->GetNormal(),
|
|
TempResultPoint)) {
|
|
FaceProjctMap[It.Position()].push_back(TempResultPoint);
|
|
str << TempResultPoint.x << " " << TempResultPoint.y << " " << TempResultPoint.z
|
|
<< std::endl;
|
|
Base::Console().log("IDX %d\n", It.Position());
|
|
|
|
if (bFirst) {
|
|
bFirst = false;
|
|
}
|
|
|
|
PointCount++;
|
|
}
|
|
}
|
|
}
|
|
|
|
str.close();
|
|
Base::Console().log("Projection map [%d facets with %d points]\n",
|
|
FaceProjctMap.size(),
|
|
PointCount);
|
|
}
|
|
|
|
|
|
bool CurveProjectorSimple::findStartPoint(const MeshKernel& MeshK,
|
|
const Base::Vector3f& Pnt,
|
|
Base::Vector3f& Rslt,
|
|
MeshCore::FacetIndex& FaceIndex)
|
|
{
|
|
Base::Vector3f TempResultPoint;
|
|
float MinLength = std::numeric_limits<float>::max();
|
|
bool bHit = false;
|
|
|
|
// go through the whole Mesh
|
|
MeshFacetIterator It(MeshK);
|
|
for (It.Init(); It.More(); It.Next()) {
|
|
// try to project (with angle) to the face
|
|
if (It->Foraminate(Pnt, It->GetNormal(), TempResultPoint)) {
|
|
// distance to the projected point
|
|
float Dist = (Pnt - TempResultPoint).Length();
|
|
if (Dist < MinLength) {
|
|
// remember the point with the closest distance
|
|
bHit = true;
|
|
MinLength = Dist;
|
|
Rslt = TempResultPoint;
|
|
FaceIndex = It.Position();
|
|
}
|
|
}
|
|
}
|
|
return bHit;
|
|
}
|
|
|
|
//**************************************************************************
|
|
//**************************************************************************
|
|
// Separator for CurveProjectorSimple classes
|
|
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|
|
|
|
|
|
CurveProjectorWithToolMesh::CurveProjectorWithToolMesh(const TopoDS_Shape& aShape,
|
|
const MeshKernel& pMesh,
|
|
MeshKernel& rToolMesh)
|
|
: CurveProjector(aShape, pMesh)
|
|
, ToolMesh(rToolMesh)
|
|
{
|
|
Do();
|
|
}
|
|
|
|
|
|
void CurveProjectorWithToolMesh::Do()
|
|
{
|
|
TopExp_Explorer Ex;
|
|
std::vector<MeshGeomFacet> cVAry;
|
|
|
|
for (Ex.Init(_Shape, TopAbs_EDGE); Ex.More(); Ex.Next()) {
|
|
const TopoDS_Edge& aEdge = TopoDS::Edge(Ex.Current());
|
|
makeToolMesh(aEdge, cVAry);
|
|
}
|
|
|
|
ToolMesh.AddFacets(cVAry);
|
|
}
|
|
|
|
// projectToNeighbours(Handle(Geom_Curve) hCurve,float pos
|
|
|
|
void CurveProjectorWithToolMesh::makeToolMesh(const TopoDS_Edge& aEdge,
|
|
std::vector<MeshGeomFacet>& cVAry)
|
|
{
|
|
Standard_Real fBegin, fEnd;
|
|
Handle(Geom_Curve) hCurve = BRep_Tool::Curve(aEdge, fBegin, fEnd);
|
|
float fLen = float(fEnd - fBegin);
|
|
Base::Vector3f cResultPoint;
|
|
|
|
unsigned long ulNbOfPoints = 15, PointCount = 0 /*,uCurFacetIdx*/;
|
|
|
|
std::vector<LineSeg> LineSegs;
|
|
|
|
MeshFacetIterator It(_Mesh);
|
|
|
|
Base::SequencerLauncher seq("Building up tool mesh...", ulNbOfPoints + 1);
|
|
|
|
std::map<MeshCore::FacetIndex, std::vector<Base::Vector3f>> FaceProjctMap;
|
|
|
|
for (unsigned long i = 0; i < ulNbOfPoints; i++) {
|
|
seq.next();
|
|
gp_Pnt gpPt = hCurve->Value(fBegin + (fLen * float(i)) / float(ulNbOfPoints - 1));
|
|
Base::Vector3f LinePoint((float)gpPt.X(), (float)gpPt.Y(), (float)gpPt.Z());
|
|
|
|
Base::Vector3f ResultNormal;
|
|
|
|
// go through the whole Mesh
|
|
for (It.Init(); It.More(); It.Next()) {
|
|
// try to project (with angle) to the face
|
|
if (It->IntersectWithLine(
|
|
Base::Vector3f((float)gpPt.X(), (float)gpPt.Y(), (float)gpPt.Z()),
|
|
It->GetNormal(),
|
|
cResultPoint)) {
|
|
if (Base::Distance(LinePoint, cResultPoint) < 0.5) {
|
|
ResultNormal += It->GetNormal();
|
|
}
|
|
}
|
|
}
|
|
LineSeg s;
|
|
s.p = Base::Vector3f((float)gpPt.X(), (float)gpPt.Y(), (float)gpPt.Z());
|
|
s.n = ResultNormal.Normalize();
|
|
LineSegs.push_back(s);
|
|
}
|
|
|
|
Base::Console().log("Projection map [%d facets with %d points]\n",
|
|
FaceProjctMap.size(),
|
|
PointCount);
|
|
|
|
|
|
// build up the new mesh
|
|
Base::Vector3f lp(std::numeric_limits<float>::max(), 0, 0), ln, p1, p2, p3, p4, p5, p6;
|
|
float ToolSize = 0.2f;
|
|
|
|
for (const auto& It2 : LineSegs) {
|
|
if (lp.x != std::numeric_limits<float>::max()) {
|
|
p1 = lp + (ln * (-ToolSize));
|
|
p2 = lp + (ln * ToolSize);
|
|
p3 = lp;
|
|
p4 = It2.p;
|
|
p5 = It2.p + (It2.n * (-ToolSize));
|
|
p6 = It2.p + (It2.n * ToolSize);
|
|
|
|
cVAry.emplace_back(p3, p2, p6);
|
|
cVAry.emplace_back(p3, p6, p4);
|
|
cVAry.emplace_back(p1, p3, p4);
|
|
cVAry.emplace_back(p1, p4, p5);
|
|
}
|
|
|
|
lp = It2.p;
|
|
ln = It2.n;
|
|
}
|
|
}
|
|
|
|
// ----------------------------------------------------------------------------
|
|
|
|
MeshProjection::MeshProjection(const MeshKernel& rMesh)
|
|
: _rcMesh(rMesh)
|
|
{}
|
|
|
|
void MeshProjection::discretize(const TopoDS_Edge& aEdge,
|
|
std::vector<Base::Vector3f>& polyline,
|
|
std::size_t minPoints) const
|
|
{
|
|
BRepAdaptor_Curve clCurve(aEdge);
|
|
|
|
Standard_Real fFirst = clCurve.FirstParameter();
|
|
Standard_Real fLast = clCurve.LastParameter();
|
|
|
|
GCPnts_UniformDeflection clDefl(clCurve, 0.01f, fFirst, fLast);
|
|
if (clDefl.IsDone() == Standard_True) {
|
|
Standard_Integer nNbPoints = clDefl.NbPoints();
|
|
for (Standard_Integer i = 1; i <= nNbPoints; i++) {
|
|
gp_Pnt gpPt = clCurve.Value(clDefl.Parameter(i));
|
|
polyline.emplace_back((float)gpPt.X(), (float)gpPt.Y(), (float)gpPt.Z());
|
|
}
|
|
}
|
|
|
|
if (polyline.size() < minPoints) {
|
|
GCPnts_UniformAbscissa clAbsc(clCurve,
|
|
static_cast<Standard_Integer>(minPoints),
|
|
fFirst,
|
|
fLast);
|
|
if (clAbsc.IsDone() == Standard_True) {
|
|
polyline.clear();
|
|
Standard_Integer nNbPoints = clAbsc.NbPoints();
|
|
for (Standard_Integer i = 1; i <= nNbPoints; i++) {
|
|
gp_Pnt gpPt = clCurve.Value(clAbsc.Parameter(i));
|
|
polyline.emplace_back((float)gpPt.X(), (float)gpPt.Y(), (float)gpPt.Z());
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
void MeshProjection::splitMeshByShape(const TopoDS_Shape& aShape, float fMaxDist) const
|
|
{
|
|
std::vector<PolyLine> rPolyLines;
|
|
projectToMesh(aShape, fMaxDist, rPolyLines);
|
|
|
|
Base::FileInfo fi("output.asc");
|
|
Base::ofstream str(fi, std::ios::out | std::ios::binary);
|
|
str.precision(4);
|
|
str.setf(std::ios::fixed | std::ios::showpoint);
|
|
for (const auto& it : rPolyLines) {
|
|
for (const auto& jt : it.points) {
|
|
str << jt.x << " " << jt.y << " " << jt.z << std::endl;
|
|
}
|
|
}
|
|
str.close();
|
|
}
|
|
|
|
bool MeshProjection::findIntersection(const Edge& edgeSegm,
|
|
const Edge& meshEdge,
|
|
const Base::Vector3f& dir,
|
|
Base::Vector3f& res) const
|
|
{
|
|
Base::Vector3f planeNormal;
|
|
planeNormal = dir.Cross(edgeSegm.cPt2 - edgeSegm.cPt1);
|
|
float dist1 = planeNormal.Dot(meshEdge.cPt1 - edgeSegm.cPt1);
|
|
float dist2 = planeNormal.Dot(meshEdge.cPt2 - edgeSegm.cPt1);
|
|
if (dist1 * dist2 < 0) {
|
|
planeNormal = dir.Cross(meshEdge.cPt2 - meshEdge.cPt1);
|
|
dist1 = planeNormal.Dot(edgeSegm.cPt1 - meshEdge.cPt1);
|
|
dist2 = planeNormal.Dot(edgeSegm.cPt2 - meshEdge.cPt1);
|
|
if (dist1 * dist2 < 0) {
|
|
// intersection detected
|
|
float t = planeNormal.Dot(meshEdge.cPt1 - edgeSegm.cPt1)
|
|
/ planeNormal.Dot(edgeSegm.cPt2 - edgeSegm.cPt1);
|
|
res = edgeSegm.cPt1 * (1 - t) + edgeSegm.cPt2 * t;
|
|
return true;
|
|
}
|
|
}
|
|
return false;
|
|
}
|
|
|
|
void MeshProjection::findSectionParameters(const TopoDS_Edge& edge,
|
|
const Base::Vector3f& dir,
|
|
std::set<double>& parameters) const
|
|
{
|
|
MeshAlgorithm clAlg(_rcMesh);
|
|
float fAvgLen = clAlg.GetAverageEdgeLength();
|
|
BRepAdaptor_Curve adapt(edge);
|
|
double edgeLen = GCPnts_AbscissaPoint::Length(adapt, Precision::Confusion());
|
|
|
|
std::vector<Base::Vector3f> polyline;
|
|
discretize(edge, polyline, std::max<size_t>(10, static_cast<size_t>(edgeLen / fAvgLen)));
|
|
if (polyline.empty()) {
|
|
return;
|
|
}
|
|
|
|
std::vector<Edge> lines;
|
|
Base::Vector3f start = polyline.front();
|
|
for (auto it = polyline.begin() + 1; it != polyline.end(); ++it) {
|
|
Edge line;
|
|
line.cPt1 = start;
|
|
line.cPt2 = *it;
|
|
start = line.cPt2;
|
|
lines.push_back(line);
|
|
}
|
|
|
|
const MeshCore::MeshFacetArray& facets = _rcMesh.GetFacets();
|
|
const MeshCore::MeshPointArray& points = _rcMesh.GetPoints();
|
|
|
|
Base::Vector3f res;
|
|
for (const auto& it : facets) {
|
|
for (int i = 0; i < 3; i++) {
|
|
Base::Vector3f pt1 = points[it._aulPoints[i]];
|
|
Base::Vector3f pt2 = points[it._aulPoints[(i + 1) % 3]];
|
|
Edge line;
|
|
line.cPt1 = pt1;
|
|
line.cPt2 = pt2;
|
|
|
|
for (auto jt : lines) {
|
|
if (findIntersection(jt, line, dir, res)) {
|
|
try {
|
|
BRepBuilderAPI_MakeVertex aBuilder(gp_Pnt(res.x, res.y, res.z));
|
|
BRepExtrema_DistShapeShape extss(aBuilder.Vertex(), edge);
|
|
if (extss.NbSolution() == 1) {
|
|
Standard_Real par;
|
|
extss.ParOnEdgeS2(1, par);
|
|
parameters.insert(par);
|
|
break;
|
|
}
|
|
}
|
|
catch (const Standard_Failure&) {
|
|
// ignore
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
void MeshProjection::projectToMesh(const TopoDS_Shape& aShape,
|
|
float fMaxDist,
|
|
std::vector<PolyLine>& rPolyLines) const
|
|
{
|
|
// calculate the average edge length and create a grid
|
|
MeshAlgorithm clAlg(_rcMesh);
|
|
float fAvgLen = clAlg.GetAverageEdgeLength();
|
|
MeshFacetGrid cGrid(_rcMesh, 5.0f * fAvgLen);
|
|
|
|
TopExp_Explorer Ex;
|
|
|
|
int iCnt = 0;
|
|
for (Ex.Init(aShape, TopAbs_EDGE); Ex.More(); Ex.Next()) {
|
|
iCnt++;
|
|
}
|
|
|
|
Base::SequencerLauncher seq("Project curve on mesh", iCnt);
|
|
|
|
for (Ex.Init(aShape, TopAbs_EDGE); Ex.More(); Ex.Next()) {
|
|
const TopoDS_Edge& aEdge = TopoDS::Edge(Ex.Current());
|
|
std::vector<SplitEdge> rSplitEdges;
|
|
projectEdgeToEdge(aEdge, fMaxDist, cGrid, rSplitEdges);
|
|
PolyLine polyline;
|
|
polyline.points.reserve(rSplitEdges.size());
|
|
for (auto it : rSplitEdges) {
|
|
polyline.points.push_back(it.cPt);
|
|
}
|
|
rPolyLines.push_back(polyline);
|
|
seq.next();
|
|
}
|
|
}
|
|
|
|
void MeshProjection::projectOnMesh(const std::vector<Base::Vector3f>& pointsIn,
|
|
const Base::Vector3f& dir,
|
|
float tolerance,
|
|
std::vector<Base::Vector3f>& pointsOut) const
|
|
{
|
|
// calculate the average edge length and create a grid
|
|
MeshAlgorithm clAlg(_rcMesh);
|
|
float fAvgLen = clAlg.GetAverageEdgeLength();
|
|
MeshFacetGrid cGrid(_rcMesh, 5.0f * fAvgLen);
|
|
|
|
// get all boundary points and edges of the mesh
|
|
std::vector<Base::Vector3f> boundaryPoints;
|
|
std::vector<MeshCore::MeshGeomEdge> boundaryEdges;
|
|
|
|
const MeshCore::MeshFacetArray& facets = _rcMesh.GetFacets();
|
|
const MeshCore::MeshPointArray& points = _rcMesh.GetPoints();
|
|
for (const auto& it : facets) {
|
|
for (int i = 0; i < 3; i++) {
|
|
if (!it.HasNeighbour(i)) {
|
|
boundaryPoints.push_back(points[it._aulPoints[i]]);
|
|
|
|
MeshCore::MeshGeomEdge edge;
|
|
edge._bBorder = true;
|
|
edge._aclPoints[0] = points[it._aulPoints[i]];
|
|
edge._aclPoints[1] = points[it._aulPoints[(i + 1) % 3]];
|
|
boundaryEdges.push_back(edge);
|
|
}
|
|
}
|
|
}
|
|
|
|
Base::SequencerLauncher seq("Project points on mesh", pointsIn.size());
|
|
|
|
for (auto it : pointsIn) {
|
|
Base::Vector3f result;
|
|
MeshCore::FacetIndex index;
|
|
if (clAlg.NearestFacetOnRay(it, dir, cGrid, result, index)) {
|
|
MeshCore::MeshGeomFacet geomFacet = _rcMesh.GetFacet(index);
|
|
if (tolerance > 0 && geomFacet.IntersectPlaneWithLine(it, dir, result)) {
|
|
if (geomFacet.IsPointOfFace(result, tolerance)) {
|
|
pointsOut.push_back(result);
|
|
}
|
|
}
|
|
else {
|
|
pointsOut.push_back(result);
|
|
}
|
|
}
|
|
else {
|
|
// go through the boundary points and check if the point can be directly projected
|
|
// onto one of them
|
|
auto boundaryPnt = std::find_if(boundaryPoints.begin(),
|
|
boundaryPoints.end(),
|
|
[&it, &dir](const Base::Vector3f& pnt) -> bool {
|
|
Base::Vector3f vec = pnt - it;
|
|
float angle = vec.GetAngle(dir);
|
|
return angle < 1e-6f;
|
|
});
|
|
|
|
if (boundaryPnt != boundaryPoints.end()) {
|
|
pointsOut.push_back(*boundaryPnt);
|
|
}
|
|
else {
|
|
// go through the boundary edges and check if the point can be directly projected
|
|
// onto one of them
|
|
Base::Vector3f result1, result2;
|
|
for (auto jt : boundaryEdges) {
|
|
jt.ClosestPointsToLine(it, dir, result1, result2);
|
|
float dot = (result1 - jt._aclPoints[0]).Dot(result1 - jt._aclPoints[1]);
|
|
Base::Vector3f vec = result1 - it;
|
|
float angle = vec.GetAngle(dir);
|
|
if (dot <= 0 && angle < 1e-6f) {
|
|
pointsOut.push_back(result1);
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
seq.next();
|
|
}
|
|
}
|
|
|
|
void MeshProjection::projectParallelToMesh(const TopoDS_Shape& aShape,
|
|
const Base::Vector3f& dir,
|
|
std::vector<PolyLine>& rPolyLines) const
|
|
{
|
|
// calculate the average edge length and create a grid
|
|
MeshAlgorithm clAlg(_rcMesh);
|
|
float fAvgLen = clAlg.GetAverageEdgeLength();
|
|
MeshFacetGrid cGrid(_rcMesh, 5.0f * fAvgLen);
|
|
TopExp_Explorer Ex;
|
|
|
|
int iCnt = 0;
|
|
for (Ex.Init(aShape, TopAbs_EDGE); Ex.More(); Ex.Next()) {
|
|
iCnt++;
|
|
}
|
|
|
|
Base::SequencerLauncher seq("Project curve on mesh", iCnt);
|
|
|
|
for (Ex.Init(aShape, TopAbs_EDGE); Ex.More(); Ex.Next()) {
|
|
const TopoDS_Edge& aEdge = TopoDS::Edge(Ex.Current());
|
|
std::vector<Base::Vector3f> points;
|
|
discretize(aEdge, points, 5);
|
|
|
|
using HitPoint = std::pair<Base::Vector3f, MeshCore::FacetIndex>;
|
|
std::vector<HitPoint> hitPoints;
|
|
using HitPoints = std::pair<HitPoint, HitPoint>;
|
|
std::vector<HitPoints> hitPointPairs;
|
|
for (auto it : points) {
|
|
Base::Vector3f result;
|
|
MeshCore::FacetIndex index;
|
|
if (clAlg.NearestFacetOnRay(it, dir, cGrid, result, index)) {
|
|
hitPoints.emplace_back(result, index);
|
|
|
|
if (hitPoints.size() > 1) {
|
|
HitPoint p1 = hitPoints[hitPoints.size() - 2];
|
|
HitPoint p2 = hitPoints[hitPoints.size() - 1];
|
|
hitPointPairs.emplace_back(p1, p2);
|
|
}
|
|
}
|
|
}
|
|
|
|
MeshCore::MeshProjection meshProjection(_rcMesh);
|
|
PolyLine polyline;
|
|
for (auto it : hitPointPairs) {
|
|
points.clear();
|
|
if (meshProjection.projectLineOnMesh(cGrid,
|
|
it.first.first,
|
|
it.first.second,
|
|
it.second.first,
|
|
it.second.second,
|
|
dir,
|
|
points)) {
|
|
polyline.points.insert(polyline.points.end(), points.begin(), points.end());
|
|
}
|
|
}
|
|
rPolyLines.push_back(polyline);
|
|
|
|
seq.next();
|
|
}
|
|
}
|
|
|
|
void MeshProjection::projectParallelToMesh(const std::vector<PolyLine>& aEdges,
|
|
const Base::Vector3f& dir,
|
|
std::vector<PolyLine>& rPolyLines) const
|
|
{
|
|
// calculate the average edge length and create a grid
|
|
MeshAlgorithm clAlg(_rcMesh);
|
|
float fAvgLen = clAlg.GetAverageEdgeLength();
|
|
MeshFacetGrid cGrid(_rcMesh, 5.0f * fAvgLen);
|
|
|
|
Base::SequencerLauncher seq("Project curve on mesh", aEdges.size());
|
|
|
|
for (const auto& it : aEdges) {
|
|
std::vector<Base::Vector3f> points = it.points;
|
|
|
|
using HitPoint = std::pair<Base::Vector3f, MeshCore::FacetIndex>;
|
|
std::vector<HitPoint> hitPoints;
|
|
using HitPoints = std::pair<HitPoint, HitPoint>;
|
|
std::vector<HitPoints> hitPointPairs;
|
|
for (auto it : points) {
|
|
Base::Vector3f result;
|
|
MeshCore::FacetIndex index;
|
|
if (clAlg.NearestFacetOnRay(it, dir, cGrid, result, index)) {
|
|
hitPoints.emplace_back(result, index);
|
|
|
|
if (hitPoints.size() > 1) {
|
|
HitPoint p1 = hitPoints[hitPoints.size() - 2];
|
|
HitPoint p2 = hitPoints[hitPoints.size() - 1];
|
|
hitPointPairs.emplace_back(p1, p2);
|
|
}
|
|
}
|
|
}
|
|
|
|
MeshCore::MeshProjection meshProjection(_rcMesh);
|
|
PolyLine polyline;
|
|
for (auto it : hitPointPairs) {
|
|
points.clear();
|
|
if (meshProjection.projectLineOnMesh(cGrid,
|
|
it.first.first,
|
|
it.first.second,
|
|
it.second.first,
|
|
it.second.second,
|
|
dir,
|
|
points)) {
|
|
polyline.points.insert(polyline.points.end(), points.begin(), points.end());
|
|
}
|
|
}
|
|
rPolyLines.push_back(polyline);
|
|
|
|
seq.next();
|
|
}
|
|
}
|
|
|
|
void MeshProjection::projectEdgeToEdge(const TopoDS_Edge& aEdge,
|
|
float fMaxDist,
|
|
const MeshFacetGrid& rGrid,
|
|
std::vector<SplitEdge>& rSplitEdges) const
|
|
{
|
|
std::vector<MeshCore::FacetIndex> auFInds;
|
|
std::map<std::pair<MeshCore::PointIndex, MeshCore::PointIndex>, std::list<MeshCore::FacetIndex>>
|
|
pEdgeToFace;
|
|
const std::vector<MeshFacet>& rclFAry = _rcMesh.GetFacets();
|
|
|
|
// search the facets in the local area of the curve
|
|
std::vector<Base::Vector3f> acPolyLine;
|
|
discretize(aEdge, acPolyLine);
|
|
|
|
MeshAlgorithm(_rcMesh).SearchFacetsFromPolyline(acPolyLine, fMaxDist, rGrid, auFInds);
|
|
// remove duplicated elements
|
|
std::sort(auFInds.begin(), auFInds.end());
|
|
auFInds.erase(std::unique(auFInds.begin(), auFInds.end()), auFInds.end());
|
|
|
|
// facet to edge
|
|
for (MeshCore::FacetIndex index : auFInds) {
|
|
const MeshFacet& rF = rclFAry[index];
|
|
for (int i = 0; i < 3; i++) {
|
|
MeshCore::PointIndex ulPt0 =
|
|
std::min<MeshCore::PointIndex>(rF._aulPoints[i], rF._aulPoints[(i + 1) % 3]);
|
|
MeshCore::PointIndex ulPt1 =
|
|
std::max<MeshCore::PointIndex>(rF._aulPoints[i], rF._aulPoints[(i + 1) % 3]);
|
|
pEdgeToFace[std::pair<MeshCore::PointIndex, MeshCore::PointIndex>(ulPt0, ulPt1)]
|
|
.push_front(index);
|
|
}
|
|
}
|
|
|
|
// sort intersection points by parameter
|
|
std::map<Standard_Real, SplitEdge> rParamSplitEdges;
|
|
|
|
BRepAdaptor_Curve clCurve(aEdge);
|
|
Standard_Real fFirst = clCurve.FirstParameter();
|
|
Standard_Real fLast = clCurve.LastParameter();
|
|
Handle(Geom_Curve) hCurve = BRep_Tool::Curve(aEdge, fFirst, fLast);
|
|
|
|
MeshPointIterator cPI(_rcMesh);
|
|
MeshFacetIterator cFI(_rcMesh);
|
|
|
|
Base::SequencerLauncher seq("Project curve on mesh", pEdgeToFace.size());
|
|
std::map<std::pair<MeshCore::PointIndex, MeshCore::PointIndex>,
|
|
std::list<MeshCore::FacetIndex>>::iterator it;
|
|
for (it = pEdgeToFace.begin(); it != pEdgeToFace.end(); ++it) {
|
|
seq.next();
|
|
|
|
// edge points
|
|
MeshCore::PointIndex uE0 = it->first.first;
|
|
cPI.Set(uE0);
|
|
Base::Vector3f cE0 = *cPI;
|
|
MeshCore::PointIndex uE1 = it->first.second;
|
|
cPI.Set(uE1);
|
|
Base::Vector3f cE1 = *cPI;
|
|
|
|
const std::list<MeshCore::FacetIndex>& auFaces = it->second;
|
|
if (auFaces.size() > 2) {
|
|
continue; // non-manifold edge -> don't handle this
|
|
}
|
|
|
|
Base::Vector3f cEdgeNormal;
|
|
for (MeshCore::FacetIndex itF : auFaces) {
|
|
cFI.Set(itF);
|
|
cEdgeNormal += cFI->GetNormal();
|
|
}
|
|
|
|
// create a plane from the edge normal and point
|
|
Base::Vector3f cPlaneNormal = cEdgeNormal % (cE1 - cE0);
|
|
Handle(Geom_Plane) hPlane =
|
|
new Geom_Plane(gp_Pln(gp_Pnt(cE0.x, cE0.y, cE0.z),
|
|
gp_Dir(cPlaneNormal.x, cPlaneNormal.y, cPlaneNormal.z)));
|
|
|
|
// get intersection of curve and plane
|
|
GeomAPI_IntCS Alg(hCurve, hPlane);
|
|
if (Alg.IsDone()) {
|
|
Standard_Integer nNbPoints = Alg.NbPoints();
|
|
if (nNbPoints == 1) {
|
|
Standard_Real fU, fV, fW;
|
|
Alg.Parameters(1, fU, fV, fW);
|
|
|
|
gp_Pnt P = Alg.Point(1);
|
|
Base::Vector3f cP0((float)P.X(), (float)P.Y(), (float)P.Z());
|
|
|
|
float l = ((cP0 - cE0) * (cE1 - cE0)) / ((cE1 - cE0) * (cE1 - cE0));
|
|
|
|
// lies the point inside the edge?
|
|
if (l >= 0.0f && l <= 1.0f) {
|
|
Base::Vector3f cSplitPoint = (1 - l) * cE0 + l * cE1;
|
|
float fDist = Base::Distance(cP0, cSplitPoint);
|
|
|
|
if (fDist <= fMaxDist) {
|
|
SplitEdge splitEdge;
|
|
splitEdge.uE0 = uE0;
|
|
splitEdge.uE1 = uE1;
|
|
splitEdge.cPt = cSplitPoint;
|
|
rParamSplitEdges[fW] = splitEdge;
|
|
}
|
|
}
|
|
}
|
|
// search for the right solution
|
|
else if (nNbPoints > 1) {
|
|
int nCntSol = 0;
|
|
Base::Vector3f cSplitPoint;
|
|
Standard_Real fSol;
|
|
Base::Vector3f cP0;
|
|
for (int j = 1; j <= nNbPoints; j++) {
|
|
Standard_Real fU, fV, fW;
|
|
Alg.Parameters(j, fU, fV, fW);
|
|
gp_Pnt P = Alg.Point(j);
|
|
cP0.Set((float)P.X(), (float)P.Y(), (float)P.Z());
|
|
|
|
float l = ((cP0 - cE0) * (cE1 - cE0)) / ((cE1 - cE0) * (cE1 - cE0));
|
|
|
|
// lies the point inside the edge?
|
|
if (l >= 0.0 && l <= 1.0) {
|
|
cSplitPoint = (1 - l) * cE0 + l * cE1;
|
|
float fDist = Base::Distance(cP0, cSplitPoint);
|
|
|
|
if (fDist <= fMaxDist) {
|
|
nCntSol++;
|
|
fSol = fW;
|
|
}
|
|
}
|
|
}
|
|
|
|
// ok, only one sensible solution
|
|
if (nCntSol == 1) {
|
|
SplitEdge splitEdge;
|
|
splitEdge.uE0 = uE0;
|
|
splitEdge.uE1 = uE1;
|
|
splitEdge.cPt = cSplitPoint;
|
|
rParamSplitEdges[fSol] = splitEdge;
|
|
}
|
|
else if (nCntSol > 1) {
|
|
Base::Console().log("More than one possible intersection points\n");
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// sorted by parameter
|
|
for (const auto& itS : rParamSplitEdges) {
|
|
rSplitEdges.push_back(itS.second);
|
|
}
|
|
}
|