00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024 #include "PreCompiled.h"
00025 #ifndef _PreComp_
00026 # ifdef FC_OS_LINUX
00027 # include <unistd.h>
00028 # endif
00029 #endif
00030
00031
00032 #include "MeshAlgos.h"
00033 #include "CurveProjector.h"
00034
00035 #include <Mod/Mesh/App/Core/MeshIO.h>
00036 #include <Mod/Mesh/App/Core/MeshKernel.h>
00037 #include <Mod/Mesh/App/Core/Iterator.h>
00038 #include <Mod/Mesh/App/Core/Algorithm.h>
00039 #include <Mod/Mesh/App/Mesh.h>
00040
00041 #include <Base/Exception.h>
00042 #include <Base/Console.h>
00043 #include <Base/Sequencer.h>
00044
00045 #include <TopExp_Explorer.hxx>
00046 #include <TopoDS.hxx>
00047 #include <Geom_Curve.hxx>
00048 #include <Geom_Plane.hxx>
00049 #include <BRep_Tool.hxx>
00050 #include <GeomAPI_IntCS.hxx>
00051
00052 using namespace MeshPart;
00053 using namespace MeshCore;
00054
00055
00056
00057
00058 CurveProjector::CurveProjector(const TopoDS_Shape &aShape, const MeshKernel &pMesh)
00059 : _Shape(aShape), _Mesh(pMesh)
00060 {
00061 }
00062
00063 void CurveProjector::writeIntersectionPointsToFile(const char *name)
00064 {
00065
00066 std::ofstream str(name, std::ios::out | std::ios::binary);
00067 str.precision(4);
00068 str.setf(std::ios::fixed | std::ios::showpoint);
00069 for (result_type::const_iterator it1 = mvEdgeSplitPoints.begin();it1!=mvEdgeSplitPoints.end();++it1) {
00070 for (std::vector<FaceSplitEdge>::const_iterator it2 = it1->second.begin();it2!=it1->second.end();++it2) {
00071 str << it2->p1.x << " " << it2->p1.y << " " << it2->p1.z << std::endl;
00072 }
00073 }
00074 str.close();
00075 }
00076
00077
00078
00079
00080
00081
00082
00083 CurveProjectorShape::CurveProjectorShape(const TopoDS_Shape &aShape, const MeshKernel &pMesh)
00084 : CurveProjector(aShape,pMesh)
00085 {
00086 Do();
00087 }
00088
00089 void CurveProjectorShape::Do(void)
00090 {
00091 TopExp_Explorer Ex;
00092 TopoDS_Shape Edge;
00093
00094 for (Ex.Init(_Shape, TopAbs_EDGE); Ex.More(); Ex.Next())
00095 {
00096 const TopoDS_Edge& aEdge = TopoDS::Edge(Ex.Current());
00097
00098
00099 projectCurve(aEdge, mvEdgeSplitPoints[aEdge]);
00100
00101 }
00102
00103 }
00104
00105
00106 void CurveProjectorShape::projectCurve( const TopoDS_Edge& aEdge,
00107 std::vector<FaceSplitEdge> &vSplitEdges)
00108 {
00109 Standard_Real fFirst, fLast;
00110 Handle(Geom_Curve) hCurve = BRep_Tool::Curve( aEdge,fFirst,fLast );
00111
00112
00113 gp_Pnt gpPt = hCurve->Value(fFirst);
00114
00115
00116 Base::Vector3f cStartPoint = Base::Vector3f((float)gpPt.X(),
00117 (float)gpPt.Y(),
00118 (float)gpPt.Z());
00119 Base::Vector3f cResultPoint, cSplitPoint, cPlanePnt, cPlaneNormal;
00120 unsigned long uStartFacetIdx,uCurFacetIdx;
00121 unsigned long uLastFacetIdx=ULONG_MAX-1;
00122 unsigned long auNeighboursIdx[3];
00123 bool GoOn;
00124
00125 if( !findStartPoint(_Mesh,cStartPoint,cResultPoint,uStartFacetIdx) )
00126 return;
00127
00128 uCurFacetIdx = uStartFacetIdx;
00129 do{
00130 MeshGeomFacet cCurFacet= _Mesh.GetFacet(uCurFacetIdx);
00131 _Mesh.GetFacetNeighbours ( uCurFacetIdx, auNeighboursIdx[0], auNeighboursIdx[1], auNeighboursIdx[2]);
00132 Base::Vector3f PointOnEdge[3];
00133
00134 GoOn = false;
00135 int NbrOfHits = 0,HitIdx=0;
00136
00137 for(int i=0; i<3; i++)
00138 {
00139
00140 if ( auNeighboursIdx[i] == uLastFacetIdx )
00141 continue;
00142
00143
00144 const Base::Vector3f& cP0 = cCurFacet._aclPoints[i];
00145 const Base::Vector3f& cP1 = cCurFacet._aclPoints[(i+1)%3];
00146
00147 if ( auNeighboursIdx[i] != ULONG_MAX )
00148 {
00149
00150 MeshGeomFacet N = _Mesh.GetFacet( auNeighboursIdx[i] );
00151 cPlaneNormal = ( N.GetNormal() + cCurFacet.GetNormal() ) % ( cP1 - cP0 );
00152 cPlanePnt = cP0;
00153 }else{
00154
00155 cPlaneNormal = cCurFacet.GetNormal() % ( cP1 - cP0 );
00156 cPlanePnt = cP0;
00157 }
00158
00159 Handle(Geom_Plane) hPlane = new Geom_Plane(gp_Pln(gp_Pnt(cPlanePnt.x,cPlanePnt.y,cPlanePnt.z),
00160 gp_Dir(cPlaneNormal.x,cPlaneNormal.y,cPlaneNormal.z)));
00161
00162 GeomAPI_IntCS Alg(hCurve,hPlane);
00163
00164 if ( Alg.IsDone() )
00165 {
00166
00167 if( Alg.NbPoints() == 1)
00168 {
00169 gp_Pnt P = Alg.Point(1);
00170 float l = ((Base::Vector3f((float)P.X(),(float)P.Y(),(float)P.Z()) - cP0)
00171 * (cP1 - cP0) ) / ((cP1 - cP0) * (cP1 - cP0));
00172
00173 if(l<0.0 || l>1.0)
00174 PointOnEdge[i] = Base::Vector3f(FLOAT_MAX,0,0);
00175 else{
00176 cSplitPoint = (1-l) * cP0 + l * cP1;
00177 PointOnEdge[i] = (1-l)*cP0 + l * cP1;
00178 NbrOfHits ++;
00179 HitIdx = i;
00180 }
00181
00182 }else if(Alg.NbPoints() == 0){
00183 PointOnEdge[i] = Base::Vector3f(FLOAT_MAX,0,0);
00184
00185 }else if(Alg.NbPoints() > 1){
00186 PointOnEdge[i] = Base::Vector3f(FLOAT_MAX,0,0);
00187 Base::Console().Log("MeshAlgos::projectCurve(): More then one intersection in Facet %ld, Edge %d\n",uCurFacetIdx,i);
00188 }
00189 }
00190 }
00191
00192 uLastFacetIdx = uCurFacetIdx;
00193
00194 if(NbrOfHits == 1)
00195 {
00196 uCurFacetIdx = auNeighboursIdx[HitIdx];
00197 FaceSplitEdge splitEdge;
00198 splitEdge.ulFaceIndex = uCurFacetIdx;
00199 splitEdge.p1 = cResultPoint;
00200 splitEdge.p2 = cSplitPoint;
00201 vSplitEdges.push_back( splitEdge );
00202 cResultPoint = cSplitPoint;
00203 GoOn = true;
00204 }else{
00205 Base::Console().Log("MeshAlgos::projectCurve(): Posibel reentry in Facet %ld\n", uCurFacetIdx);
00206 }
00207
00208 if( uCurFacetIdx == uStartFacetIdx )
00209 GoOn = false;
00210
00211 }while(GoOn);
00212
00213 }
00214
00215 bool CurveProjectorShape::findStartPoint(const MeshKernel &MeshK,const Base::Vector3f &Pnt,Base::Vector3f &Rslt,unsigned long &FaceIndex)
00216 {
00217 Base::Vector3f TempResultPoint;
00218 float MinLength = FLOAT_MAX;
00219 bool bHit = false;
00220
00221
00222 MeshFacetIterator It(MeshK);
00223 for(It.Init();It.More();It.Next())
00224 {
00225
00226 if(It->Foraminate (Pnt, It->GetNormal(), TempResultPoint) )
00227 {
00228
00229 float Dist = (Pnt-TempResultPoint).Length();
00230 if(Dist < MinLength)
00231 {
00232
00233 bHit = true;
00234 MinLength = Dist;
00235 Rslt = TempResultPoint;
00236 FaceIndex = It.Position();
00237 }
00238 }
00239 }
00240 return bHit;
00241 }
00242
00243
00244
00245
00246
00247
00248
00249
00250 CurveProjectorSimple::CurveProjectorSimple(const TopoDS_Shape &aShape, const MeshKernel &pMesh)
00251 : CurveProjector(aShape,pMesh)
00252 {
00253 Do();
00254 }
00255
00256
00257 void CurveProjectorSimple::Do(void)
00258 {
00259 TopExp_Explorer Ex;
00260 TopoDS_Shape Edge;
00261
00262 std::vector<Base::Vector3f> vEdgePolygon;
00263
00264 for (Ex.Init(_Shape, TopAbs_EDGE); Ex.More(); Ex.Next())
00265 {
00266 const TopoDS_Edge& aEdge = TopoDS::Edge(Ex.Current());
00267
00268
00269
00270 projectCurve(aEdge,vEdgePolygon, mvEdgeSplitPoints[aEdge]);
00271
00272 }
00273
00274 }
00275
00276
00277 void CurveProjectorSimple::GetSampledCurves( const TopoDS_Edge& aEdge, std::vector<Base::Vector3f>& rclPoints, unsigned long ulNbOfPoints)
00278 {
00279 rclPoints.clear();
00280
00281 Standard_Real fBegin, fEnd;
00282
00283 Handle(Geom_Curve) hCurve = BRep_Tool::Curve(aEdge,fBegin,fEnd);
00284 float fLen = float(fEnd - fBegin);
00285
00286 for (unsigned long i = 0; i < ulNbOfPoints; i++)
00287 {
00288 gp_Pnt gpPt = hCurve->Value(fBegin + (fLen * float(i)) / float(ulNbOfPoints-1));
00289 rclPoints.push_back(Base::Vector3f((float)gpPt.X(),
00290 (float)gpPt.Y(),
00291 (float)gpPt.Z()));
00292 }
00293 }
00294
00295
00296
00297
00298 void CurveProjectorSimple::projectCurve( const TopoDS_Edge& aEdge,
00299 const std::vector<Base::Vector3f> &rclPoints,
00300 std::vector<FaceSplitEdge> &vSplitEdges)
00301 {
00302 Base::Vector3f TempResultPoint;
00303 bool bFirst = true;
00304
00305
00306
00307 Standard_Real fBegin, fEnd;
00308 Handle(Geom_Curve) hCurve = BRep_Tool::Curve(aEdge,fBegin,fEnd);
00309 float fLen = float(fEnd - fBegin);
00310
00311 unsigned long ulNbOfPoints = 1000,PointCount=0;
00312
00313 MeshFacetIterator It(_Mesh);
00314
00315 Base::SequencerLauncher seq("Building up projection map...", ulNbOfPoints+1);
00316 std::ofstream str("projected.asc", std::ios::out | std::ios::binary);
00317 str.precision(4);
00318 str.setf(std::ios::fixed | std::ios::showpoint);
00319
00320 std::map<unsigned long,std::vector<Base::Vector3f> > FaceProjctMap;
00321
00322 for (unsigned long i = 0; i <= ulNbOfPoints; i++)
00323 {
00324 seq.next();
00325 gp_Pnt gpPt = hCurve->Value(fBegin + (fLen * float(i)) / float(ulNbOfPoints-1));
00326
00327
00328 for(It.Init();It.More();It.Next())
00329 {
00330
00331 if (It->IntersectWithLine (Base::Vector3f((float)gpPt.X(),(float)gpPt.Y(),(float)gpPt.Z()),
00332 It->GetNormal(), TempResultPoint))
00333 {
00334 FaceProjctMap[It.Position()].push_back(TempResultPoint);
00335 str << TempResultPoint.x << " "
00336 << TempResultPoint.y << " "
00337 << TempResultPoint.z << std::endl;
00338 Base::Console().Log("IDX %d\n",It.Position());
00339
00340 if(bFirst){
00341 bFirst = false;
00342 }
00343
00344 PointCount++;
00345 }
00346 }
00347 }
00348
00349 str.close();
00350 Base::Console().Log("Projection map [%d facets with %d points]\n",FaceProjctMap.size(),PointCount);
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395 }
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440
00441
00442
00443
00444
00445
00446
00447
00448
00449
00450
00451
00452
00453
00454
00455
00456
00457
00458
00459
00460
00461
00462
00463
00464
00465
00466
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504
00505
00506
00507
00508
00509
00510
00511
00512
00513
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523
00524
00525
00526
00527
00528
00529
00530 bool CurveProjectorSimple::findStartPoint(const MeshKernel &MeshK,const Base::Vector3f &Pnt,Base::Vector3f &Rslt,unsigned long &FaceIndex)
00531 {
00532 Base::Vector3f TempResultPoint;
00533 float MinLength = FLOAT_MAX;
00534 bool bHit = false;
00535
00536
00537 MeshFacetIterator It(MeshK);
00538 for(It.Init();It.More();It.Next())
00539 {
00540
00541 if(It->Foraminate (Pnt, It->GetNormal(), TempResultPoint) )
00542 {
00543
00544 float Dist = (Pnt-TempResultPoint).Length();
00545 if(Dist < MinLength)
00546 {
00547
00548 bHit = true;
00549 MinLength = Dist;
00550 Rslt = TempResultPoint;
00551 FaceIndex = It.Position();
00552 }
00553 }
00554 }
00555 return bHit;
00556 }
00557
00558
00559
00560
00561
00562
00563
00564 CurveProjectorWithToolMesh::CurveProjectorWithToolMesh(const TopoDS_Shape &aShape, const MeshKernel &pMesh,MeshKernel &rToolMesh)
00565 : CurveProjector(aShape,pMesh),ToolMesh(rToolMesh)
00566 {
00567 Do();
00568 }
00569
00570
00571 void CurveProjectorWithToolMesh::Do(void)
00572 {
00573 TopExp_Explorer Ex;
00574 TopoDS_Shape Edge;
00575 std::vector<MeshGeomFacet> cVAry;
00576
00577 std::vector<Base::Vector3f> vEdgePolygon;
00578
00579 for (Ex.Init(_Shape, TopAbs_EDGE); Ex.More(); Ex.Next())
00580 {
00581 const TopoDS_Edge& aEdge = TopoDS::Edge(Ex.Current());
00582
00583 makeToolMesh(aEdge,cVAry);
00584
00585 }
00586
00587 ToolMesh.AddFacets(cVAry);
00588
00589 }
00590
00591
00592
00593
00594 void CurveProjectorWithToolMesh::makeToolMesh( const TopoDS_Edge& aEdge,std::vector<MeshGeomFacet> &cVAry )
00595 {
00596 Standard_Real fBegin, fEnd;
00597 Handle(Geom_Curve) hCurve = BRep_Tool::Curve(aEdge,fBegin,fEnd);
00598 float fLen = float(fEnd - fBegin);
00599 Base::Vector3f cResultPoint;
00600
00601 unsigned long ulNbOfPoints = 15,PointCount=0;
00602
00603 std::vector<LineSeg> LineSegs;
00604
00605 MeshFacetIterator It(_Mesh);
00606
00607 Base::SequencerLauncher seq("Building up tool mesh...", ulNbOfPoints+1);
00608
00609 std::map<unsigned long,std::vector<Base::Vector3f> > FaceProjctMap;
00610
00611 for (unsigned long i = 0; i < ulNbOfPoints; i++)
00612 {
00613 seq.next();
00614 gp_Pnt gpPt = hCurve->Value(fBegin + (fLen * float(i)) / float(ulNbOfPoints-1));
00615 Base::Vector3f LinePoint((float)gpPt.X(),
00616 (float)gpPt.Y(),
00617 (float)gpPt.Z());
00618
00619 Base::Vector3f ResultNormal;
00620
00621
00622 for(It.Init();It.More();It.Next())
00623 {
00624
00625 if (It->IntersectWithLine (Base::Vector3f((float)gpPt.X(),(float)gpPt.Y(),(float)gpPt.Z()),
00626 It->GetNormal(), cResultPoint) )
00627 {
00628 if(Base::Distance(LinePoint,cResultPoint) < 0.5)
00629 ResultNormal += It->GetNormal();
00630 }
00631 }
00632 LineSeg s;
00633 s.p = Base::Vector3f((float)gpPt.X(),
00634 (float)gpPt.Y(),
00635 (float)gpPt.Z());
00636 s.n = ResultNormal.Normalize();
00637 LineSegs.push_back(s);
00638 }
00639
00640 Base::Console().Log("Projection map [%d facets with %d points]\n",FaceProjctMap.size(),PointCount);
00641
00642
00643
00644 Base::Vector3f lp(FLOAT_MAX,0,0), ln, p1, p2, p3, p4,p5,p6;
00645 float ToolSize = 0.2f;
00646
00647 for (std::vector<LineSeg>::iterator It2=LineSegs.begin(); It2!=LineSegs.end();++It2)
00648 {
00649 if(lp.x != FLOAT_MAX)
00650 {
00651 p1 = lp + (ln * (-ToolSize));
00652 p2 = lp + (ln * ToolSize);
00653 p3 = lp;
00654 p4 = (*It2).p;
00655 p5 = (*It2).p + ((*It2).n * (-ToolSize));
00656 p6 = (*It2).p + ((*It2).n * ToolSize);
00657
00658 cVAry.push_back(MeshGeomFacet(p3,p2,p6));
00659 cVAry.push_back(MeshGeomFacet(p3,p6,p4));
00660 cVAry.push_back(MeshGeomFacet(p1,p3,p4));
00661 cVAry.push_back(MeshGeomFacet(p1,p4,p5));
00662
00663 }
00664
00665 lp = (*It2).p;
00666 ln = (*It2).n;
00667 }
00668
00669
00670
00671 }