Algorithm.cpp

Go to the documentation of this file.
00001 /***************************************************************************
00002  *   Copyright (c) 2005 Imetric 3D GmbH                                    *
00003  *                                                                         *
00004  *   This file is part of the FreeCAD CAx development system.              *
00005  *                                                                         *
00006  *   This library is free software; you can redistribute it and/or         *
00007  *   modify it under the terms of the GNU Library General Public           *
00008  *   License as published by the Free Software Foundation; either          *
00009  *   version 2 of the License, or (at your option) any later version.      *
00010  *                                                                         *
00011  *   This library  is distributed in the hope that it will be useful,      *
00012  *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
00013  *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
00014  *   GNU Library General Public License for more details.                  *
00015  *                                                                         *
00016  *   You should have received a copy of the GNU Library General Public     *
00017  *   License along with this library; see the file COPYING.LIB. If not,    *
00018  *   write to the Free Software Foundation, Inc., 59 Temple Place,         *
00019  *   Suite 330, Boston, MA  02111-1307, USA                                *
00020  *                                                                         *
00021  ***************************************************************************/
00022 
00023 
00024 #include "PreCompiled.h"
00025 
00026 #ifndef _PreComp_
00027 # include <algorithm>
00028 #endif
00029 
00030 #include "Algorithm.h"
00031 #include "Approximation.h"
00032 #include "Elements.h"
00033 #include "Iterator.h"
00034 #include "Grid.h"
00035 #include "Triangulation.h"
00036 
00037 #include <Base/Console.h>
00038 #include <Base/Sequencer.h>
00039 
00040 using namespace MeshCore;
00041 using Base::BoundBox3f;
00042 using Base::BoundBox2D;
00043 using Base::Polygon2D;
00044 
00045 
00046 bool MeshAlgorithm::IsVertexVisible (const Base::Vector3f &rcVertex, const Base::Vector3f &rcView, const MeshFacetGrid &rclGrid ) const
00047 {
00048   Base::Vector3f cDirection = rcVertex-rcView;
00049   float fDistance = cDirection.Length();
00050   Base::Vector3f cIntsct; unsigned long uInd;
00051 
00052   // search for the nearest facet to rcView in direction to rcVertex
00053   if ( NearestFacetOnRay( rcView, cDirection, /*1.2f*fDistance,*/ rclGrid, cIntsct, uInd) )
00054   {
00055     // now check if the facet overlays the point
00056     float fLen = Base::Distance( rcView, cIntsct );
00057     if ( fLen < fDistance )
00058     {
00059       // is it the same point?
00060       if ( Base::Distance(rcVertex, cIntsct) > 0.001f )
00061       {
00062         // ok facet overlays the vertex
00063         return false;
00064       }
00065     }
00066   }
00067 
00068   return true; // no facet between the two points
00069 }
00070 
00071 bool MeshAlgorithm::NearestFacetOnRay (const Base::Vector3f &rclPt, const Base::Vector3f &rclDir, Base::Vector3f &rclRes,
00072                                        unsigned long &rulFacet) const
00073 {
00074     Base::Vector3f clProj, clRes;
00075     bool bSol = false;
00076     unsigned long ulInd = 0;
00077 
00078     // langsame Ausfuehrung ohne Grid
00079     MeshFacetIterator  clFIter(_rclMesh);
00080     for (clFIter.Init(); clFIter.More(); clFIter.Next()) {
00081         if (clFIter->Foraminate( rclPt, rclDir, clRes ) == true) {
00082             if (bSol == false) { // erste Loesung
00083                 bSol   = true;
00084                 clProj = clRes;
00085                 ulInd  = clFIter.Position();
00086             }
00087             else {  // liegt Punkt naeher
00088                 if ((clRes - rclPt).Length() < (clProj - rclPt).Length()) {
00089                     clProj = clRes;
00090                     ulInd  = clFIter.Position();
00091                 }
00092             }
00093         }
00094     }
00095 
00096     if (bSol) {
00097         rclRes   = clProj;
00098         rulFacet = ulInd;
00099     }
00100 
00101     return bSol;
00102 }
00103 
00104 bool MeshAlgorithm::NearestFacetOnRay (const Base::Vector3f &rclPt, const Base::Vector3f &rclDir, const MeshFacetGrid &rclGrid,
00105                                        Base::Vector3f &rclRes, unsigned long &rulFacet) const
00106 {
00107     std::vector<unsigned long> aulFacets;
00108     MeshGridIterator clGridIter(rclGrid);
00109 
00110     if (clGridIter.InitOnRay(rclPt, rclDir, aulFacets) == true) {
00111         if (RayNearestField(rclPt, rclDir, aulFacets, rclRes, rulFacet) == false) {
00112             aulFacets.clear();
00113             while (clGridIter.NextOnRay(aulFacets) == true) {
00114                 if (RayNearestField(rclPt, rclDir, aulFacets, rclRes, rulFacet) == true)
00115                     return true;
00116             }
00117         }
00118         else
00119             return true;
00120     }
00121 
00122     return false;
00123 }
00124 
00125 bool MeshAlgorithm::NearestFacetOnRay (const Base::Vector3f &rclPt, const Base::Vector3f &rclDir, float fMaxSearchArea,
00126                                        const MeshFacetGrid &rclGrid, Base::Vector3f &rclRes, unsigned long &rulFacet) const
00127 {
00128     std::vector<unsigned long> aulFacets;
00129     MeshGridIterator  clGridIter(rclGrid);
00130 
00131     if (clGridIter.InitOnRay(rclPt, rclDir, fMaxSearchArea, aulFacets) == true) {
00132         if (RayNearestField(rclPt, rclDir, aulFacets, rclRes, rulFacet, 1.75f) == false) {
00133             aulFacets.clear();
00134             while (clGridIter.NextOnRay(aulFacets) == true) {
00135                 if (RayNearestField(rclPt, rclDir, aulFacets, rclRes, rulFacet, 1.75f) == true)
00136                     return true;
00137             }
00138         }
00139         else
00140             return true;
00141     }
00142 
00143     return false;
00144 }
00145 
00146 bool MeshAlgorithm::NearestFacetOnRay (const Base::Vector3f &rclPt, const Base::Vector3f &rclDir, const std::vector<unsigned long> &raulFacets,
00147                                        Base::Vector3f &rclRes, unsigned long &rulFacet) const
00148 {
00149     Base::Vector3f  clProj, clRes;
00150     bool bSol = false;
00151     unsigned long ulInd = 0;
00152 
00153     for (std::vector<unsigned long>::const_iterator pI = raulFacets.begin(); pI != raulFacets.end(); pI++) {
00154         MeshGeomFacet rclSFacet = _rclMesh.GetFacet(*pI);
00155         if (rclSFacet.Foraminate(rclPt, rclDir, clRes) == true) {
00156             if (bSol == false) {// erste Loesung
00157                 bSol   = true;
00158                 clProj = clRes;
00159                 ulInd  = *pI;
00160             }
00161             else {  // liegt Punkt naeher
00162                 if ((clRes - rclPt).Length() < (clProj - rclPt).Length()) {
00163                     clProj = clRes;
00164                     ulInd  = *pI;
00165                 }
00166             }
00167         }
00168     }
00169 
00170     if (bSol) {
00171         rclRes   = clProj;
00172         rulFacet = ulInd;
00173     }
00174 
00175     return bSol;
00176 }
00177 
00178 bool MeshAlgorithm::RayNearestField (const Base::Vector3f &rclPt, const Base::Vector3f &rclDir, const std::vector<unsigned long> &raulFacets,
00179                                      Base::Vector3f &rclRes, unsigned long &rulFacet, float fMaxAngle) const
00180 {
00181     Base::Vector3f  clProj, clRes;
00182     bool bSol = false;
00183     unsigned long ulInd = 0;
00184 
00185     for (std::vector<unsigned long>::const_iterator pF = raulFacets.begin(); pF != raulFacets.end(); pF++) {
00186         if (_rclMesh.GetFacet(*pF).Foraminate(rclPt, rclDir, clRes/*, fMaxAngle*/) == true) {
00187             if (bSol == false) { // erste Loesung
00188                 bSol   = true;
00189                 clProj = clRes;
00190                 ulInd  = *pF;
00191             }
00192             else {  // liegt Punkt naeher
00193                 if ((clRes - rclPt).Length() < (clProj - rclPt).Length()) {
00194                     clProj = clRes;
00195                     ulInd  = *pF;
00196                 }
00197             }
00198         }
00199     }
00200 
00201     if (bSol) {
00202         rclRes   = clProj;
00203         rulFacet = ulInd;
00204     }
00205 
00206     return bSol;
00207 }
00208 
00209 bool MeshAlgorithm::FirstFacetToVertex(const Base::Vector3f &rPt, float fMaxDistance, const MeshFacetGrid &rGrid, unsigned long &uIndex) const
00210 {
00211     const float fEps = 0.001f;
00212 
00213     bool found = false;
00214     std::vector<unsigned long> facets;
00215 
00216     // get the facets of the grid the point lies into
00217     rGrid.GetElements(rPt, facets);
00218 
00219     // Check all facets inside the grid if the point is part of it
00220     for (std::vector<unsigned long>::iterator it = facets.begin(); it != facets.end(); ++it) {
00221         MeshGeomFacet cFacet = this->_rclMesh.GetFacet(*it);
00222         if (cFacet.IsPointOfFace(rPt, fMaxDistance)) {
00223             found = true;
00224             uIndex = *it;
00225             break;
00226         }
00227         else {
00228             // if not then check the distance to the border of the triangle
00229             Base::Vector3f res = rPt;
00230             float fDist;
00231             unsigned short uSide;
00232             cFacet.ProjectPointToPlane(res);
00233             cFacet.NearestEdgeToPoint(res, fDist, uSide);
00234             if (fDist < fEps) {
00235                 found = true;
00236                 uIndex = *it;
00237                 break;
00238             }
00239         }
00240     }
00241 
00242     return found;
00243 }
00244 
00245 float MeshAlgorithm::GetAverageEdgeLength() const
00246 {
00247     float fLen = 0.0f;
00248     MeshFacetIterator cF(_rclMesh);
00249     for (cF.Init(); cF.More(); cF.Next()) {
00250         for (int i=0; i<3; i++)
00251             fLen += Base::Distance(cF->_aclPoints[i], cF->_aclPoints[(i+1)%3]);
00252     }
00253 
00254     fLen = fLen / (3.0f * _rclMesh.CountFacets() );
00255     return fLen;
00256 }
00257 
00258 void MeshAlgorithm::GetMeshBorders (std::list<std::vector<Base::Vector3f> > &rclBorders) const
00259 {
00260     std::vector<unsigned long> aulAllFacets(_rclMesh.CountFacets());
00261     unsigned long k = 0;
00262     for (std::vector<unsigned long>::iterator pI = aulAllFacets.begin(); pI != aulAllFacets.end(); pI++)
00263         *pI = k++;
00264 
00265     GetFacetBorders(aulAllFacets, rclBorders);
00266 }
00267 
00268 void MeshAlgorithm::GetMeshBorders (std::list<std::vector<unsigned long> > &rclBorders) const
00269 {
00270     std::vector<unsigned long> aulAllFacets(_rclMesh.CountFacets());
00271     unsigned long k = 0;
00272     for (std::vector<unsigned long>::iterator pI = aulAllFacets.begin(); pI != aulAllFacets.end(); pI++)
00273         *pI = k++;
00274 
00275     GetFacetBorders(aulAllFacets, rclBorders, true);
00276 }
00277 
00278 void MeshAlgorithm::GetFacetBorders (const std::vector<unsigned long> &raulInd, std::list<std::vector<Base::Vector3f> > &rclBorders) const
00279 {
00280 #if 1
00281   const MeshPointArray &rclPAry = _rclMesh._aclPointArray;
00282   std::list<std::vector<unsigned long> > aulBorders;
00283 
00284   GetFacetBorders (raulInd, aulBorders, true);
00285   for ( std::list<std::vector<unsigned long> >::iterator it = aulBorders.begin(); it != aulBorders.end(); ++it )
00286   {
00287     std::vector<Base::Vector3f> boundary;
00288     boundary.reserve( it->size() );
00289 
00290     for ( std::vector<unsigned long>::iterator jt = it->begin(); jt != it->end(); ++jt )
00291       boundary.push_back(rclPAry[*jt]);
00292 
00293     rclBorders.push_back( boundary );
00294   }
00295 
00296 #else
00297   const MeshFacetArray &rclFAry = _rclMesh._aclFacetArray;
00298 
00299   // alle Facets markieren die in der Indizie-Liste vorkommen
00300   ResetFacetFlag(MeshFacet::VISIT);
00301   for (std::vector<unsigned long>::const_iterator pIter = raulInd.begin(); pIter != raulInd.end(); pIter++)
00302     rclFAry[*pIter].SetFlag(MeshFacet::VISIT);
00303 
00304   std::list<std::pair<unsigned long, unsigned long> >  aclEdges;
00305   // alle Randkanten suchen und ablegen (unsortiert)
00306   for (std::vector<unsigned long>::const_iterator pIter2 = raulInd.begin(); pIter2 != raulInd.end(); pIter2++)
00307   {
00308     const MeshFacet  &rclFacet = rclFAry[*pIter2];
00309     for (int i = 0; i < 3; i++)
00310     {
00311       unsigned long ulNB = rclFacet._aulNeighbours[i];
00312       if (ulNB != ULONG_MAX)
00313       {
00314         if (rclFAry[ulNB].IsFlag(MeshFacet::VISIT) == true)
00315           continue;
00316       }
00317 
00318       aclEdges.push_back(rclFacet.GetEdge(i));
00319     }
00320   }
00321 
00322   if (aclEdges.size() == 0)
00323     return; // no borders found (=> solid)
00324 
00325   // Kanten aus der unsortieren Kantenliste suchen
00326   const MeshPointArray &rclPAry = _rclMesh._aclPointArray;
00327   unsigned long              ulFirst, ulLast;
00328   std::list<Base::Vector3f>   clBorder;
00329   ulFirst = aclEdges.begin()->first;
00330   ulLast  = aclEdges.begin()->second;
00331 
00332   aclEdges.erase(aclEdges.begin());
00333   clBorder.push_back(rclPAry[ulFirst]);
00334   clBorder.push_back(rclPAry[ulLast]);
00335 
00336   while (aclEdges.size() > 0)
00337   {
00338     // naechste anliegende Kante suchen
00339     std::list<std::pair<unsigned long, unsigned long> >::iterator pEI;
00340     for (pEI = aclEdges.begin(); pEI != aclEdges.end(); pEI++)
00341     {
00342       if (pEI->first == ulLast)
00343       {
00344         ulLast = pEI->second;
00345         clBorder.push_back(rclPAry[ulLast]);
00346         aclEdges.erase(pEI);
00347         break;
00348       }
00349       else if (pEI->second == ulLast)
00350       {
00351         ulLast = pEI->first;
00352         clBorder.push_back(rclPAry[ulLast]);
00353         aclEdges.erase(pEI);
00354         break;
00355       }
00356       else if (pEI->first == ulFirst)
00357       {
00358         ulFirst = pEI->second;
00359         clBorder.push_front(rclPAry[ulFirst]);
00360         aclEdges.erase(pEI);
00361         break;
00362       }
00363       else if (pEI->second == ulFirst)
00364       {
00365         ulFirst = pEI->first;
00366         clBorder.push_front(rclPAry[ulFirst]);
00367         aclEdges.erase(pEI);
00368         break;
00369       }
00370     }
00371     if ((pEI == aclEdges.end()) || (ulLast == ulFirst))
00372     {  // keine weitere Kante gefunden bzw. Polylinie geschlossen
00373       rclBorders.push_back(std::vector<Base::Vector3f>(clBorder.begin(), clBorder.end()));
00374       clBorder.clear();
00375 
00376       if (aclEdges.size() > 0)
00377       {  // neue Border anfangen
00378         ulFirst = aclEdges.begin()->first;
00379         ulLast  = aclEdges.begin()->second;
00380         aclEdges.erase(aclEdges.begin());
00381         clBorder.push_back(rclPAry[ulFirst]);
00382         clBorder.push_back(rclPAry[ulLast]);
00383       }
00384     }
00385   }
00386 #endif
00387 }
00388 
00389 void MeshAlgorithm::GetFacetBorders (const std::vector<unsigned long> &raulInd, 
00390                                      std::list<std::vector<unsigned long> > &rclBorders,
00391                                      bool ignoreOrientation) const
00392 {
00393     const MeshFacetArray &rclFAry = _rclMesh._aclFacetArray;
00394 
00395     // mark all facets that are in the indices list
00396     ResetFacetFlag(MeshFacet::VISIT);
00397     for (std::vector<unsigned long>::const_iterator it = raulInd.begin(); it != raulInd.end(); ++it)
00398         rclFAry[*it].SetFlag(MeshFacet::VISIT);
00399 
00400     // collect all boundary edges (unsorted)
00401     std::list<std::pair<unsigned long, unsigned long> >  aclEdges;
00402     for (std::vector<unsigned long>::const_iterator it = raulInd.begin(); it != raulInd.end(); ++it) {
00403         const MeshFacet  &rclFacet = rclFAry[*it];
00404         for (int i = 0; i < 3; i++) {
00405             unsigned long ulNB = rclFacet._aulNeighbours[i];
00406             if (ulNB != ULONG_MAX) {
00407                 if (rclFAry[ulNB].IsFlag(MeshFacet::VISIT) == true)
00408                     continue;
00409                 }
00410 
00411                 aclEdges.push_back(rclFacet.GetEdge(i));
00412             }
00413         }
00414 
00415     if (aclEdges.size() == 0)
00416         return; // no borders found (=> solid)
00417 
00418     // search for edges in the unsorted list
00419     unsigned long              ulFirst, ulLast;
00420     std::list<unsigned long>   clBorder;
00421     ulFirst = aclEdges.begin()->first;
00422     ulLast  = aclEdges.begin()->second;
00423 
00424     aclEdges.erase(aclEdges.begin());
00425     clBorder.push_back(ulFirst);
00426     clBorder.push_back(ulLast);
00427 
00428     while (aclEdges.size() > 0) {
00429         // get adjacent edge
00430         std::list<std::pair<unsigned long, unsigned long> >::iterator pEI;
00431         for (pEI = aclEdges.begin(); pEI != aclEdges.end(); ++pEI) {
00432             if (pEI->first == ulLast) {
00433                 ulLast = pEI->second;
00434                 clBorder.push_back(ulLast);
00435                 aclEdges.erase(pEI);
00436                 pEI = aclEdges.begin();
00437                 break;
00438             }
00439             else if (pEI->second == ulFirst) {
00440                 ulFirst = pEI->first;
00441                 clBorder.push_front(ulFirst);
00442                 aclEdges.erase(pEI);
00443                 pEI = aclEdges.begin();
00444                 break;
00445             }
00446             // Note: Using this might result into boundaries with wrong orientation.
00447             // But if the mesh has some facets with wrong orientation we might get
00448             // broken boundary curves.
00449             else if (pEI->second == ulLast && ignoreOrientation) {
00450                 ulLast = pEI->first;
00451                 clBorder.push_back(ulLast);
00452                 aclEdges.erase(pEI);
00453                 pEI = aclEdges.begin();
00454                 break;
00455             }
00456             else if (pEI->first == ulFirst && ignoreOrientation) {
00457                 ulFirst = pEI->second;
00458                 clBorder.push_front(ulFirst);
00459                 aclEdges.erase(pEI);
00460                 pEI = aclEdges.begin();
00461                 break;
00462             }
00463         }
00464 
00465         // Note: Calling erase on list iterators doesn't force a re-allocation and
00466         // thus doesn't invalidate the iterator itself, only the referenced object
00467         if ((pEI == aclEdges.end()) || aclEdges.empty() || (ulLast == ulFirst)) {
00468             // no further edge found or closed polyline, respectively
00469             rclBorders.push_back(std::vector<unsigned long>(clBorder.begin(), clBorder.end()));
00470             clBorder.clear();
00471 
00472             if (aclEdges.size() > 0) {
00473                 // start new boundary
00474                 ulFirst = aclEdges.begin()->first;
00475                 ulLast  = aclEdges.begin()->second;
00476                 aclEdges.erase(aclEdges.begin());
00477                 clBorder.push_back(ulFirst);
00478                 clBorder.push_back(ulLast);
00479             }
00480         }
00481     }
00482 }
00483 
00484 void MeshAlgorithm::GetMeshBorder(unsigned long uFacet, std::list<unsigned long>& rBorder) const
00485 {
00486     const MeshFacetArray &rFAry = _rclMesh._aclFacetArray;
00487     std::list<std::pair<unsigned long, unsigned long> >  openEdges;
00488     if (uFacet >= rFAry.size())
00489         return;
00490     // add the open edge to the beginning of the list
00491     MeshFacetArray::_TConstIterator face = rFAry.begin() + uFacet;
00492     for (int i = 0; i < 3; i++)
00493     {
00494         if (face->_aulNeighbours[i] == ULONG_MAX)
00495             openEdges.push_back(face->GetEdge(i));
00496     }
00497 
00498     if (openEdges.empty())
00499         return; // facet is not a border facet
00500     
00501     for (MeshFacetArray::_TConstIterator it = rFAry.begin(); it != rFAry.end(); ++it)
00502     {
00503         if (it == face)
00504             continue;
00505         for (int i = 0; i < 3; i++)
00506         {
00507             if (it->_aulNeighbours[i] == ULONG_MAX)
00508                 openEdges.push_back(it->GetEdge(i));
00509         }
00510     }
00511 
00512     // Start with the edge that is associated to uFacet
00513     unsigned long ulFirst = openEdges.begin()->first;
00514     unsigned long ulLast  = openEdges.begin()->second;
00515 
00516     openEdges.erase(openEdges.begin());
00517     rBorder.push_back(ulFirst);
00518     rBorder.push_back(ulLast);
00519 
00520     while (ulLast != ulFirst)
00521     {
00522         // find adjacent edge
00523         std::list<std::pair<unsigned long, unsigned long> >::iterator pEI;
00524         for (pEI = openEdges.begin(); pEI != openEdges.end(); pEI++)
00525         {
00526             if (pEI->first == ulLast)
00527             {
00528                 ulLast = pEI->second;
00529                 rBorder.push_back(ulLast);
00530                 openEdges.erase(pEI);
00531                 break;
00532             }
00533             else if (pEI->second == ulFirst)
00534             {
00535                 ulFirst = pEI->first;
00536                 rBorder.push_front(ulFirst);
00537                 openEdges.erase(pEI);
00538                 break;
00539             }
00540         }
00541 
00542         // cannot close the border
00543         if (pEI == openEdges.end())
00544             break;
00545     }
00546 }
00547 
00548 void MeshAlgorithm::SplitBoundaryLoops( std::list<std::vector<unsigned long> >& aBorders )
00549 {
00550     // Count the number of open edges for each point
00551     std::map<unsigned long, int> openPointDegree;
00552     for (MeshFacetArray::_TConstIterator jt = _rclMesh._aclFacetArray.begin();
00553         jt != _rclMesh._aclFacetArray.end(); ++jt) {
00554         for (int i=0; i<3; i++) {
00555             if (jt->_aulNeighbours[i] == ULONG_MAX) {
00556                 openPointDegree[jt->_aulPoints[i]]++;
00557                 openPointDegree[jt->_aulPoints[(i+1)%3]]++;
00558             }
00559         }
00560     }
00561 
00562     // go through all boundaries and split them if needed
00563     std::list<std::vector<unsigned long> > aSplitBorders;
00564     for (std::list<std::vector<unsigned long> >::iterator it = aBorders.begin();
00565         it != aBorders.end(); ++it) {
00566         bool split=false;
00567         for (std::vector<unsigned long>::iterator jt = it->begin(); jt != it->end(); ++jt) {
00568             // two (ore more) boundaries meet in one non-manifold point
00569             if (openPointDegree[*jt] > 2) {
00570                 split = true;
00571                 break;
00572             }
00573         }
00574 
00575         if (!split)
00576             aSplitBorders.push_back( *it );
00577         else
00578             SplitBoundaryLoops( *it, aSplitBorders );
00579     }
00580 
00581     aBorders = aSplitBorders;
00582 }
00583 
00584 void MeshAlgorithm::SplitBoundaryLoops(const std::vector<unsigned long>& rBound,
00585                                        std::list<std::vector<unsigned long> >& aBorders)
00586 {
00587     std::map<unsigned long, int> aPtDegree;
00588     std::vector<unsigned long> cBound;
00589     for (std::vector<unsigned long>::const_iterator it = rBound.begin(); it != rBound.end(); ++it) {
00590         int deg = (aPtDegree[*it]++);
00591         if (deg > 0) {
00592             for (std::vector<unsigned long>::iterator jt = cBound.begin(); jt != cBound.end(); ++jt) {
00593                 if (*jt == *it) {
00594                     std::vector<unsigned long> cBoundLoop;
00595                     cBoundLoop.insert(cBoundLoop.end(), jt, cBound.end());
00596                     cBoundLoop.push_back(*it);
00597                     cBound.erase(jt, cBound.end());
00598                     aBorders.push_back(cBoundLoop);
00599                     (aPtDegree[*it]--);
00600                     break;
00601                 }
00602             }
00603         }
00604 
00605         cBound.push_back(*it);
00606     }
00607 }
00608 
00609 bool MeshAlgorithm::FillupHole(const std::vector<unsigned long>& boundary, 
00610                                AbstractPolygonTriangulator& cTria, 
00611                                MeshFacetArray& rFaces, MeshPointArray& rPoints,
00612                                int level, const MeshRefPointToFacets* pP2FStructure) const
00613 {
00614     if (boundary.front() == boundary.back()) {
00615         // first and last vertex are identical
00616         if (boundary.size() < 4)
00617             return false; // something strange
00618     }
00619     else if (boundary.size() < 3) {
00620         return false; // something strange
00621     }
00622 
00623     // Get a facet as reference coordinate system
00624     MeshGeomFacet rTriangle;
00625     MeshFacet rFace;
00626     unsigned long refPoint0 = *(boundary.begin());
00627     unsigned long refPoint1 = *(boundary.begin()+1);
00628     if (pP2FStructure) {
00629         const std::set<unsigned long>& ring = (*pP2FStructure)[refPoint0];
00630         bool found = false;
00631         for (std::set<unsigned long>::const_iterator it = ring.begin(); it != ring.end() && !found; ++it) {
00632             for (int i=0; i<3; i++) {
00633                 if (pP2FStructure->getFacet(*it)->_aulPoints[i] == refPoint1) {
00634                     rFace = *pP2FStructure->getFacet(*it);
00635                     rTriangle = _rclMesh.GetFacet(rFace);
00636                     found = true;
00637                     break;
00638                 }
00639             }
00640         }
00641     } else {
00642         bool ready = false;
00643         for (MeshFacetArray::_TConstIterator it = _rclMesh._aclFacetArray.begin(); it != _rclMesh._aclFacetArray.end(); ++it) {
00644             for (int i=0; i<3; i++) {
00645                 if ((it->_aulPoints[i] == refPoint0) && (it->_aulPoints[(i+1)%3] == refPoint1)) {
00646                     rFace = *it;
00647                     rTriangle = _rclMesh.GetFacet(*it);
00648                     ready = true;
00649                     break;
00650                 }
00651             }
00652 
00653             if (ready)
00654                 break;
00655         }
00656     }
00657 
00658     // add points to the polygon
00659     std::vector<Base::Vector3f> polygon;
00660     for (std::vector<unsigned long>::const_iterator jt = boundary.begin(); jt != boundary.end(); ++jt) {
00661         polygon.push_back(_rclMesh._aclPointArray[*jt]);
00662         rPoints.push_back(_rclMesh._aclPointArray[*jt]);
00663     }
00664 
00665     // remove the last added point if it is duplicated
00666     if (boundary.front() == boundary.back()) {
00667         polygon.pop_back();
00668         rPoints.pop_back();
00669     }
00670 
00671     // There is no easy way to check whether the boundary is interior (a hole) or exterior before performing the triangulation.
00672     // Afterwards we can compare the normals of the created triangles with the z-direction of our local coordinate system.
00673     // If the scalar product is positive it was a hole, otherwise not.
00674     cTria.SetPolygon(polygon);
00675 
00676     std::vector<Base::Vector3f> surf_pts = cTria.GetPolygon();
00677     if (pP2FStructure && level > 0) {
00678         std::set<unsigned long> index = pP2FStructure->NeighbourPoints(boundary, level);
00679         for (std::set<unsigned long>::iterator it = index.begin(); it != index.end(); ++it) {
00680             Base::Vector3f pt(_rclMesh._aclPointArray[*it]);
00681             surf_pts.push_back(pt);
00682         }
00683     }
00684 
00685     if (cTria.TriangulatePolygon()) {
00686         // if we have enough points then we fit a surface through the points and project
00687         // the added points onto this surface
00688         cTria.ProjectOntoSurface(surf_pts);
00689         // get the facets and add the additional points to the array
00690         rFaces.insert(rFaces.end(), cTria.GetFacets().begin(), cTria.GetFacets().end());
00691         std::vector<Base::Vector3f> newVertices = cTria.AddedPoints();
00692         for (std::vector<Base::Vector3f>::iterator pt = newVertices.begin(); pt != newVertices.end(); ++pt) {
00693             rPoints.push_back((*pt));
00694         }
00695 
00696         // Unfortunately, the CDT algorithm does not care about the orientation of the polygon so we cannot rely on the normal
00697         // criterion to decide whether it's a hole or not.
00698         //
00699         std::vector<MeshFacet> faces = cTria.GetFacets();
00700         
00701         // Special case handling for a hole with three edges: the resulting facet might be coincident with the 
00702         // reference facet
00703         if (faces.size()==1){
00704             MeshFacet first;
00705             first._aulPoints[0] = boundary[faces.front()._aulPoints[0]];
00706             first._aulPoints[1] = boundary[faces.front()._aulPoints[1]];
00707             first._aulPoints[2] = boundary[faces.front()._aulPoints[2]];
00708             if ( first.IsEqual(rFace) ) {
00709                 rFaces.clear();
00710                 rPoints.clear();
00711                 cTria.Discard();
00712                 return false;
00713             }
00714         }
00715 
00716         // Get the new neighbour to our reference facet
00717         MeshFacet facet;
00718         unsigned short ref_side = rFace.Side(refPoint0, refPoint1);
00719         unsigned short tri_side = USHRT_MAX;
00720         if (ref_side < USHRT_MAX) {
00721             for (std::vector<MeshFacet>::iterator it = faces.begin(); it != faces.end(); ++it) {
00722                 tri_side = it->Side(0, 1);
00723                 if (tri_side < USHRT_MAX) {
00724                     facet = *it;
00725                     break;
00726                 }
00727             }
00728         }
00729 
00730         // in case the reference facet has not an open edge print a log message 
00731         if (ref_side == USHRT_MAX || tri_side == USHRT_MAX) {
00732             Base::Console().Log("MeshAlgorithm::FillupHole: Expected open edge for facet <%d, %d, %d>\n", 
00733                 rFace._aulPoints[0], rFace._aulPoints[1], rFace._aulPoints[2]);
00734             rFaces.clear();
00735             rPoints.clear();
00736             cTria.Discard();
00737             return false;
00738         }
00739 
00740         MeshGeomFacet triangle;
00741         triangle._aclPoints[0] = rPoints[facet._aulPoints[0]];
00742         triangle._aclPoints[1] = rPoints[facet._aulPoints[1]];
00743         triangle._aclPoints[2] = rPoints[facet._aulPoints[2]];
00744 
00745         // Now we have two adjacent triangles which we check for overlaps.
00746         // Therefore we build a separation plane that must separate the two diametrically opposed points.
00747         Base::Vector3f planeNormal = rTriangle.GetNormal() % (rTriangle._aclPoints[(ref_side+1)%3]-rTriangle._aclPoints[ref_side]);
00748         Base::Vector3f planeBase = rTriangle._aclPoints[ref_side];
00749         Base::Vector3f ref_point = rTriangle._aclPoints[(ref_side+2)%3];
00750         Base::Vector3f tri_point = triangle._aclPoints[(tri_side+2)%3];
00751 
00752         float ref_dist = (ref_point - planeBase) * planeNormal;
00753         float tri_dist = (tri_point - planeBase) * planeNormal;
00754         if (ref_dist * tri_dist > 0.0f) {
00755             rFaces.clear();
00756             rPoints.clear();
00757             cTria.Discard();
00758             return false;
00759         }
00760 
00761         // we know to have filled a polygon, now check for the orientation
00762         if ( triangle.GetNormal() * rTriangle.GetNormal() <= 0.0f ) {
00763             for (MeshFacetArray::_TIterator it = rFaces.begin(); it != rFaces.end(); ++it)
00764                 it->FlipNormal();
00765         }
00766 
00767         return true;
00768     }
00769 
00770     return false;
00771 }
00772 
00773 void MeshAlgorithm::SetFacetsProperty(const std::vector<unsigned long> &raulInds, const std::vector<unsigned long> &raulProps) const
00774 {
00775     if (raulInds.size() != raulProps.size()) return;
00776 
00777     std::vector<unsigned long>::const_iterator iP = raulProps.begin();
00778     for (std::vector<unsigned long>::const_iterator i = raulInds.begin(); i != raulInds.end(); i++, iP++)
00779         _rclMesh._aclFacetArray[*i].SetProperty(*iP);
00780 }
00781 
00782 void MeshAlgorithm::SetFacetsFlag (const std::vector<unsigned long> &raulInds, MeshFacet::TFlagType tF) const
00783 {
00784     for (std::vector<unsigned long>::const_iterator i = raulInds.begin(); i != raulInds.end(); i++)
00785         _rclMesh._aclFacetArray[*i].SetFlag(tF);
00786 }
00787 
00788 void MeshAlgorithm::SetPointsFlag (const std::vector<unsigned long> &raulInds, MeshPoint::TFlagType tF) const 
00789 {
00790     for (std::vector<unsigned long>::const_iterator i = raulInds.begin(); i != raulInds.end(); i++)
00791         _rclMesh._aclPointArray[*i].SetFlag(tF);
00792 }
00793 
00794 void MeshAlgorithm::GetFacetsFlag (std::vector<unsigned long> &raulInds, MeshFacet::TFlagType tF) const
00795 {
00796     raulInds.reserve(raulInds.size() + CountFacetFlag(tF));
00797     MeshFacetArray::_TConstIterator beg = _rclMesh._aclFacetArray.begin();
00798     MeshFacetArray::_TConstIterator end = _rclMesh._aclFacetArray.end();
00799     for (MeshFacetArray::_TConstIterator it = beg; it != end; ++it) {
00800         if (it->IsFlag(tF))
00801             raulInds.push_back(it-beg);
00802     }
00803 }
00804 
00805 void MeshAlgorithm::GetPointsFlag (std::vector<unsigned long> &raulInds, MeshPoint::TFlagType tF) const
00806 {
00807     raulInds.reserve(raulInds.size() + CountPointFlag(tF));
00808     MeshPointArray::_TConstIterator beg = _rclMesh._aclPointArray.begin();
00809     MeshPointArray::_TConstIterator end = _rclMesh._aclPointArray.end();
00810     for (MeshPointArray::_TConstIterator it = beg; it != end; ++it) {
00811         if (it->IsFlag(tF))
00812             raulInds.push_back(it-beg);
00813     }
00814 }
00815 
00816 void MeshAlgorithm::ResetFacetsFlag (const std::vector<unsigned long> &raulInds, MeshFacet::TFlagType tF) const
00817 {
00818     for (std::vector<unsigned long>::const_iterator i = raulInds.begin(); i != raulInds.end(); i++)
00819         _rclMesh._aclFacetArray[*i].ResetFlag(tF);
00820 }
00821 
00822 void MeshAlgorithm::ResetPointsFlag (const std::vector<unsigned long> &raulInds, MeshPoint::TFlagType tF) const
00823 {
00824     for (std::vector<unsigned long>::const_iterator i = raulInds.begin(); i != raulInds.end(); i++)
00825         _rclMesh._aclPointArray[*i].ResetFlag(tF);
00826 }
00827 
00828 void MeshAlgorithm::SetFacetFlag (MeshFacet::TFlagType tF) const
00829 {
00830     _rclMesh._aclFacetArray.SetFlag(tF);
00831 }
00832 
00833 void MeshAlgorithm::SetPointFlag (MeshPoint::TFlagType tF) const
00834 {
00835     _rclMesh._aclPointArray.SetFlag(tF);
00836 }
00837 
00838 void MeshAlgorithm::ResetFacetFlag (MeshFacet::TFlagType tF) const
00839 {
00840     _rclMesh._aclFacetArray.ResetFlag(tF);
00841 }
00842 
00843 void MeshAlgorithm::ResetPointFlag (MeshPoint::TFlagType tF) const
00844 {
00845     _rclMesh._aclPointArray.ResetFlag(tF);
00846 }
00847 
00848 unsigned long MeshAlgorithm::CountFacetFlag (MeshFacet::TFlagType tF) const
00849 {
00850     return std::count_if(_rclMesh._aclFacetArray.begin(), _rclMesh._aclFacetArray.end(),
00851                     std::bind2nd(MeshIsFlag<MeshFacet>(), tF));
00852 }
00853 
00854 unsigned long MeshAlgorithm::CountPointFlag (MeshPoint::TFlagType tF) const
00855 {
00856     return std::count_if(_rclMesh._aclPointArray.begin(), _rclMesh._aclPointArray.end(),
00857                     std::bind2nd(MeshIsFlag<MeshPoint>(), tF));
00858 }
00859 
00860 void MeshAlgorithm::GetFacetsFromToolMesh( const MeshKernel& rToolMesh, const Base::Vector3f& rcDir, std::vector<unsigned long> &raclCutted ) const
00861 {
00862     MeshFacetIterator cFIt(_rclMesh);
00863     MeshFacetIterator cTIt(rToolMesh);
00864 
00865     BoundBox3f cBB = rToolMesh.GetBoundBox();
00866 
00867     Base::SequencerLauncher seq("Check facets...", _rclMesh.CountFacets());
00868 
00869     // check all facets
00870     Base::Vector3f tmp;
00871     for (cFIt.Init(); cFIt.More(); cFIt.Next()) {
00872         // check each point of each facet
00873         for (int i=0; i<3; i++) {
00874             // at least the point must be inside the bounding box of the tool mesh
00875             if (cBB.IsInBox( cFIt->_aclPoints[i])) {
00876                 // should not cause performance problems since the tool mesh is usually rather lightweight
00877                 int ct=0;
00878                 for (cTIt.Init(); cTIt.More(); cTIt.Next()) {
00879                     if (cTIt->IsPointOfFace( cFIt->_aclPoints[i], MeshPoint::epsilon())) {
00880                         ct=1;
00881                         break; // the point lies on the tool mesh
00882                     }
00883                     else if (cTIt->Foraminate( cFIt->_aclPoints[i], rcDir, tmp)) {
00884                         // check if the intersection point lies in direction rcDir of the considered point
00885                         if ((tmp - cFIt->_aclPoints[i]) * rcDir > 0)
00886                             ct++;
00887                     }
00888                 }
00889 
00890                 // odd number => point is inside the tool mesh
00891                 if (ct % 2 == 1) {
00892                     raclCutted.push_back( cFIt.Position() );
00893                     break;
00894                 }
00895             }
00896         }
00897 
00898         seq.next();
00899     }
00900 }
00901 
00902 void MeshAlgorithm::GetFacetsFromToolMesh(const MeshKernel& rToolMesh, const Base::Vector3f& rcDir,
00903                                           const MeshFacetGrid& rGrid, std::vector<unsigned long> &raclCutted) const
00904 {
00905     // iterator over grid structure
00906     MeshGridIterator clGridIter(rGrid);
00907     BoundBox3f cBB = rToolMesh.GetBoundBox();
00908     Base::Vector3f tmp;
00909 
00910     MeshFacetIterator cFIt(_rclMesh);
00911     MeshFacetIterator cTIt(rToolMesh);
00912     MeshAlgorithm cToolAlg(rToolMesh);
00913 
00914     // To speed up the algorithm we use the grid built up from the associated mesh. For each grid
00915     // element we check whether it lies completely inside or outside the toolmesh or even intersect
00916     // with the toolmesh. So we can reduce the number of facets with further tests dramatically.
00917     // If the grid box is outside the toolmesh all the facets inside can be skipped. If the grid
00918     // box is inside the toolmesh all facets are stored with no further tests because they must
00919     // also lie inside the toolmesh. Finally, if the grid box intersect with the toolmesh we must
00920     // also check for each whether it intersect we the toolmesh as well.
00921     std::vector<unsigned long> aulInds;
00922     for (clGridIter.Init(); clGridIter.More(); clGridIter.Next()) {
00923         int ret = cToolAlg.Surround(clGridIter.GetBoundBox(), rcDir);
00924 
00925         // the box is completely inside the toolmesh
00926         if (ret == 1) {
00927             // these facets can be removed without more checks
00928             clGridIter.GetElements(raclCutted);
00929         } 
00930         // the box intersect with toolmesh
00931         else if (ret == 0) {
00932             // these facets must be tested for intersectons with the toolmesh
00933             clGridIter.GetElements(aulInds);
00934         }
00935         // the box is outside the toolmesh but this could still mean that the triangles
00936         // inside the grid intersect with the toolmesh
00937         else if (ret == -1) {
00938             // these facets must be tested for intersectons with the toolmesh
00939             clGridIter.GetElements(aulInds);
00940         }
00941     }
00942 
00943     // remove duplicates
00944     std::sort(aulInds.begin(), aulInds.end());
00945     aulInds.erase(std::unique(aulInds.begin(), aulInds.end()), aulInds.end());
00946     std::sort(raclCutted.begin(), raclCutted.end());
00947     raclCutted.erase(std::unique(raclCutted.begin(), raclCutted.end()), raclCutted.end());
00948 
00949     Base::SequencerLauncher seq("Check facets...", aulInds.size());
00950 
00951     // check all facets
00952     for (std::vector<unsigned long>::iterator it = aulInds.begin(); it != aulInds.end(); ++it) {
00953         cFIt.Set(*it);
00954 
00955         // check each point of each facet
00956         for (int i=0; i<3; i++) {
00957             // at least the point must be inside the bounding box of the tool mesh
00958             if (cBB.IsInBox(cFIt->_aclPoints[i])) {
00959                 // should not cause performance problems since the tool mesh is usually rather lightweight
00960                 int ct=0;
00961                 for (cTIt.Init(); cTIt.More(); cTIt.Next()) {
00962                     if (cTIt->IsPointOfFace(cFIt->_aclPoints[i], MeshPoint::epsilon())) {
00963                         ct=1;
00964                         break; // the point lies on the tool mesh
00965                     }
00966                     else if (cTIt->Foraminate(cFIt->_aclPoints[i], rcDir, tmp)) {
00967                         // check if the intersection point lies in direction rcDir of the considered point
00968                         if ((tmp - cFIt->_aclPoints[i]) * rcDir > 0)
00969                             ct++;
00970                     }
00971                 }
00972 
00973                 // odd number => point is inside the tool mesh
00974                 if (ct % 2 == 1) {
00975                     raclCutted.push_back(cFIt.Position());
00976                     break;
00977                 }
00978             }
00979         }
00980 
00981         seq.next();
00982     }
00983 
00984     // remove duplicates
00985     std::sort(raclCutted.begin(), raclCutted.end());
00986     raclCutted.erase(std::unique(raclCutted.begin(), raclCutted.end()), raclCutted.end());
00987 }
00988 
00989 int MeshAlgorithm::Surround(const Base::BoundBox3f& rBox, const Base::Vector3f& rcDir)
00990 {
00991     Base::Vector3f pt1, pt2, tmp;
00992     const BoundBox3f& cBB = _rclMesh.GetBoundBox();
00993 
00994     // at least both boxes intersect
00995     if (cBB && rBox) {
00996         // check for intersections with the actual mesh
00997         Base::Vector3f cCorner[8] = {
00998         Base::Vector3f(rBox.MinX,rBox.MinY,rBox.MinZ), Base::Vector3f(rBox.MaxX,rBox.MinY,rBox.MinZ), 
00999         Base::Vector3f(rBox.MaxX,rBox.MaxY,rBox.MinZ), Base::Vector3f(rBox.MinX,rBox.MaxY,rBox.MinZ),
01000         Base::Vector3f(rBox.MinX,rBox.MinY,rBox.MaxZ), Base::Vector3f(rBox.MaxX,rBox.MinY,rBox.MaxZ), 
01001         Base::Vector3f(rBox.MaxX,rBox.MaxY,rBox.MaxZ), Base::Vector3f(rBox.MinX,rBox.MaxY,rBox.MaxZ)};
01002 
01003         MeshFacetIterator cTIt(_rclMesh);
01004 
01005         // triangulation of the box
01006         int triangles[36] = {
01007             0,1,2,0,2,3,
01008             0,1,5,0,5,4,
01009             0,4,7,0,7,3,
01010             6,7,4,6,4,5,
01011             6,2,3,6,3,7,
01012             6,1,2,6,5,1
01013         };
01014 
01015         std::vector<MeshGeomFacet> cFacet(12);
01016         int id=0;
01017         for (int ii=0; ii<12; ii++) {
01018             cFacet[ii]._aclPoints[0]=cCorner[triangles[id++]]; 
01019             cFacet[ii]._aclPoints[1]=cCorner[triangles[id++]]; 
01020             cFacet[ii]._aclPoints[2]=cCorner[triangles[id++]];
01021         }
01022 
01023         // check for intersections of the box with the mesh
01024         for (std::vector<MeshGeomFacet>::iterator it = cFacet.begin(); it != cFacet.end(); ++it) {
01025             for (cTIt.Init(); cTIt.More(); cTIt.Next()) {
01026                 int ret = cTIt->IntersectWithFacet(*it,pt1, pt2);
01027 
01028                 // the box intersects the mesh?
01029                 if (ret != 0)
01030                     return 0; // => no more investigations required
01031             }
01032         }
01033 
01034         // Now we know that the box doesn't intersect with the mesh. This means that either the box
01035         // is completely inside or outside the mesh. To check this we test one point of the box
01036         // whether it is inside or outside.
01037         int ct=0;
01038         for (cTIt.Init(); cTIt.More(); cTIt.Next()) {
01039             if (cTIt->IsPointOfFace(cCorner[0], MeshPoint::epsilon())) {
01040                 ct=1;
01041                 break; // the point lies on the tool mesh
01042             }
01043             else if (cTIt->Foraminate(cCorner[0], rcDir, tmp)) {
01044                 // check if the intersection point lies in direction rcDir of the considered point
01045                 if ((tmp - cCorner[0]) * rcDir > 0)
01046                     ct++;
01047             }
01048         }
01049 
01050         // odd number => point (i.e. the box) is inside the mesh, even number => point is outside the mesh
01051         return (ct % 2 == 1) ? 1 : -1;
01052     }
01053 
01054     // no intersection the box is outside the mesh
01055     return -1;
01056 }
01057 
01058 void MeshAlgorithm::CheckFacets(const MeshFacetGrid& rclGrid, const Base::ViewProjMethod* pclProj, const Base::Polygon2D& rclPoly,
01059                                 bool bInner, std::vector<unsigned long> &raulFacets) const
01060 {
01061     std::vector<unsigned long>::iterator it;
01062     MeshFacetIterator clIter(_rclMesh, 0);
01063     Base::Vector3f clPt2d;
01064     Base::Vector3f clGravityOfFacet;
01065     bool bNoPointInside;
01066 
01067     // Falls true, verwende Grid auf Mesh, um Suche zu beschleunigen
01068     if (bInner)
01069     {
01070         BoundBox3f clBBox3d;
01071         BoundBox2D clViewBBox, clPolyBBox;
01072         std::vector<unsigned long> aulAllElements;
01073 
01074         //B-Box des Polygons
01075         clPolyBBox = rclPoly.CalcBoundBox();
01076         // Iterator fuer die zu durchsuchenden B-Boxen des Grids
01077         MeshGridIterator clGridIter(rclGrid);
01078         // alle B-Boxen durchlaufen und die Facets speichern
01079         for (clGridIter.Init(); clGridIter.More(); clGridIter.Next())
01080         {
01081             clBBox3d = clGridIter.GetBoundBox();
01082             clViewBBox = clBBox3d.ProjectBox(pclProj);
01083             if (clViewBBox || clPolyBBox)
01084             {
01085                 // alle Elemente in AllElements sammeln
01086                 clGridIter.GetElements(aulAllElements);
01087             }
01088         }
01089         // doppelte Elemente wieder entfernen
01090         std::sort(aulAllElements.begin(), aulAllElements.end());
01091         aulAllElements.erase(std::unique(aulAllElements.begin(), aulAllElements.end()), aulAllElements.end());
01092 
01093         Base::SequencerLauncher seq( "Check facets", aulAllElements.size() );
01094 
01095         for (it = aulAllElements.begin(); it != aulAllElements.end(); it++)
01096         {
01097             bNoPointInside = true;
01098             clGravityOfFacet.Set(0.0f, 0.0f, 0.0f);
01099             MeshGeomFacet rclFacet = _rclMesh.GetFacet(*it);
01100             for (int j=0; j<3; j++)
01101             {
01102                 clPt2d = pclProj->operator()(rclFacet._aclPoints[j]);
01103                 clGravityOfFacet += clPt2d;
01104                 if (rclPoly.Contains(Base::Vector2D(clPt2d.x, clPt2d.y)) == bInner)
01105                 {
01106                     raulFacets.push_back(*it);
01107                     bNoPointInside = false;
01108                     break;
01109                 }
01110             }
01111 
01112             // if no facet point is inside the polygon then check also the gravity
01113             if (bNoPointInside == true)
01114             {
01115               clGravityOfFacet *= 1.0f/3.0f;
01116 
01117               if (rclPoly.Contains(Base::Vector2D(clGravityOfFacet.x, clGravityOfFacet.y)) == bInner)
01118                  raulFacets.push_back(*it);
01119             }
01120 
01121             seq.next();
01122         }
01123     }
01124     // Dreiecke ausserhalb schneiden, dann alles durchsuchen
01125     else
01126     {
01127       Base::SequencerLauncher seq("Check facets", _rclMesh.CountFacets());
01128       for (clIter.Init(); clIter.More(); clIter.Next())
01129       {
01130           for (int j=0; j<3; j++)
01131           {
01132               clPt2d = pclProj->operator()(clIter->_aclPoints[j]);
01133               if (rclPoly.Contains(Base::Vector2D(clPt2d.x, clPt2d.y)) == bInner)
01134               {
01135                   raulFacets.push_back(clIter.Position());
01136                   break;
01137               }
01138           }
01139           seq.next();
01140       }
01141     }
01142 }
01143 
01144 void MeshAlgorithm::CheckFacets(const Base::ViewProjMethod* pclProj, const Base::Polygon2D& rclPoly,
01145                                 bool bInner, std::vector<unsigned long> &raulFacets) const
01146 {
01147     const MeshPointArray& p = _rclMesh.GetPoints();
01148     const MeshFacetArray& f = _rclMesh.GetFacets();
01149     Base::SequencerLauncher seq("Check facets", f.size());
01150     Base::Vector3f pt2d;
01151     unsigned long index=0;
01152     for (MeshFacetArray::_TConstIterator it = f.begin(); it != f.end(); ++it,++index) {
01153         for (int i = 0; i < 3; i++) {
01154             pt2d = (*pclProj)(p[it->_aulPoints[i]]);
01155             if (rclPoly.Contains(Base::Vector2D(pt2d.x, pt2d.y)) == bInner) {
01156                 raulFacets.push_back(index);
01157                 break;
01158             }
01159         }
01160         seq.next();
01161     }
01162 }
01163 
01164 float MeshAlgorithm::Surface (void) const
01165 {
01166   float              fTotal = 0.0f;
01167   MeshFacetIterator clFIter(_rclMesh);
01168 
01169   for (clFIter.Init(); clFIter.More(); clFIter.Next())
01170     fTotal +=  clFIter->Area();
01171   
01172   return fTotal;
01173 }
01174 
01175 void MeshAlgorithm::SubSampleByDist (float fDist, std::vector<Base::Vector3f> &rclPoints) const
01176 {
01177     rclPoints.clear();
01178     MeshFacetIterator clFIter(_rclMesh);
01179     for (clFIter.Init(); clFIter.More(); clFIter.Next()) {
01180         size_t k = rclPoints.size();
01181         clFIter->SubSample(fDist, rclPoints);
01182         if (rclPoints.size() == k)
01183             rclPoints.push_back(clFIter->GetGravityPoint()); // min. add middle point
01184     }
01185 }
01186 
01187 void MeshAlgorithm::SubSampleAllPoints (std::vector<Base::Vector3f> &rclPoints) const
01188 {
01189     rclPoints.clear();
01190 
01191     // Add all Points
01192     //
01193     MeshPointIterator clPIter(_rclMesh);
01194     for (clPIter.Init(); clPIter.More(); clPIter.Next()) {
01195         rclPoints.push_back(*clPIter);
01196     }
01197 }
01198 
01199 void MeshAlgorithm::SubSampleByCount (unsigned long ulCtPoints, std::vector<Base::Vector3f> &rclPoints) const
01200 {
01201     float fDist = float(sqrt(Surface() / float(ulCtPoints)));
01202     SubSampleByDist(fDist, rclPoints);
01203 }
01204 
01205 void MeshAlgorithm::SearchFacetsFromPolyline (const std::vector<Base::Vector3f> &rclPolyline, float fRadius,
01206                                               const MeshFacetGrid& rclGrid, std::vector<unsigned long> &rclResultFacetsIndices) const
01207 {
01208   rclResultFacetsIndices.clear();
01209   if ( rclPolyline.size() < 3 )
01210     return; // no polygon defined
01211 
01212   std::set<unsigned long>  aclFacets;
01213   for (std::vector<Base::Vector3f>::const_iterator pV = rclPolyline.begin(); pV < (rclPolyline.end() - 1); pV++)
01214   {
01215     const Base::Vector3f &rclP0 = *pV, &rclP1 = *(pV + 1);
01216    
01217     // BB eines Polyline-Segments
01218     BoundBox3f clSegmBB(rclP0.x, rclP0.y, rclP0.z, rclP0.x, rclP0.y, rclP0.z);
01219     clSegmBB &= rclP1;
01220     clSegmBB.Enlarge(fRadius);  // BB um Suchradius vergroessern
01221 
01222     std::vector<unsigned long> aclBBFacets;
01223     unsigned long k = rclGrid.Inside(clSegmBB, aclBBFacets, false);
01224     for (unsigned long i = 0; i < k; i++)
01225     {
01226       if (_rclMesh.GetFacet(aclBBFacets[i]).DistanceToLineSegment(rclP0, rclP1) < fRadius)
01227         aclFacets.insert(aclBBFacets[i]);
01228     }
01229   }
01230 
01231   rclResultFacetsIndices.insert(rclResultFacetsIndices.begin(), aclFacets.begin(), aclFacets.end());
01232 }
01233 
01234 void MeshAlgorithm::CutBorderFacets (std::vector<unsigned long> &raclFacetIndices, unsigned short usLevel) const
01235 {
01236   std::vector<unsigned long> aclToDelete;
01237 
01238   CheckBorderFacets(raclFacetIndices, aclToDelete, usLevel);
01239 
01240   // alle gefunden "Rand"-Facetsindizes" aus dem Array loeschen
01241   std::vector<unsigned long>  aclResult;
01242   std::set<unsigned long>     aclTmp(aclToDelete.begin(), aclToDelete.end());
01243 
01244   for (std::vector<unsigned long>::iterator pI = raclFacetIndices.begin(); pI != raclFacetIndices.end(); pI++)
01245   {
01246     if (aclTmp.find(*pI) == aclTmp.end())
01247       aclResult.push_back(*pI);
01248   }
01249 
01250   raclFacetIndices = aclResult;
01251 }
01252 
01253 unsigned long MeshAlgorithm::CountBorderEdges() const
01254 {
01255     unsigned long cnt=0;
01256     const MeshFacetArray &rclFAry = _rclMesh._aclFacetArray;
01257     MeshFacetArray::_TConstIterator end = rclFAry.end();
01258     for (MeshFacetArray::_TConstIterator it = rclFAry.begin(); it != end; ++it) {
01259         for (int i=0; i<3; i++) {
01260             if (it->_aulNeighbours[i] == ULONG_MAX)
01261                 cnt++;
01262         }
01263     }
01264 
01265     return cnt;
01266 }
01267 
01268 void MeshAlgorithm::CheckBorderFacets (const std::vector<unsigned long> &raclFacetIndices, std::vector<unsigned long> &raclResultIndices, unsigned short usLevel) const
01269 {
01270   ResetFacetFlag(MeshFacet::TMP0);
01271   SetFacetsFlag(raclFacetIndices, MeshFacet::TMP0);
01272 
01273   const MeshFacetArray &rclFAry = _rclMesh._aclFacetArray;
01274 
01275   for (unsigned short usL = 0; usL < usLevel; usL++)
01276   {
01277     for (std::vector<unsigned long>::const_iterator pF = raclFacetIndices.begin(); pF != raclFacetIndices.end(); pF++)
01278     {
01279       for (int i = 0; i < 3; i++)
01280       {
01281         unsigned long ulNB = rclFAry[*pF]._aulNeighbours[i];
01282         if (ulNB == ULONG_MAX)
01283         {
01284           raclResultIndices.push_back(*pF);
01285           rclFAry[*pF].ResetFlag(MeshFacet::TMP0);
01286           continue;
01287         }
01288         if (rclFAry[ulNB].IsFlag(MeshFacet::TMP0) == false)
01289         {
01290           raclResultIndices.push_back(*pF);
01291           rclFAry[*pF].ResetFlag(MeshFacet::TMP0);
01292           continue;
01293         }
01294       }
01295     }
01296   }
01297 }
01298 
01299 void MeshAlgorithm::GetBorderPoints (const std::vector<unsigned long> &raclFacetIndices, std::set<unsigned long> &raclResultPointsIndices) const
01300 {
01301   ResetFacetFlag(MeshFacet::TMP0);
01302   SetFacetsFlag(raclFacetIndices, MeshFacet::TMP0);
01303 
01304   const MeshFacetArray &rclFAry = _rclMesh._aclFacetArray;
01305 
01306   for (std::vector<unsigned long>::const_iterator pF = raclFacetIndices.begin(); pF != raclFacetIndices.end(); pF++)
01307   {
01308     for (int i = 0; i < 3; i++)
01309     {
01310       const MeshFacet &rclFacet = rclFAry[*pF];
01311       unsigned long      ulNB     = rclFacet._aulNeighbours[i];
01312       if (ulNB == ULONG_MAX)
01313       {
01314         raclResultPointsIndices.insert(rclFacet._aulPoints[i]);
01315         raclResultPointsIndices.insert(rclFacet._aulPoints[(i+1)%3]);
01316         continue;
01317       }
01318       if (rclFAry[ulNB].IsFlag(MeshFacet::TMP0) == false)
01319       {
01320         raclResultPointsIndices.insert(rclFacet._aulPoints[i]);
01321         raclResultPointsIndices.insert(rclFacet._aulPoints[(i+1)%3]);
01322         continue;
01323       }
01324     }
01325   }
01326 }
01327 
01328 bool MeshAlgorithm::NearestPointFromPoint (const Base::Vector3f &rclPt, unsigned long &rclResFacetIndex, Base::Vector3f &rclResPoint) const
01329 {
01330   if (_rclMesh.CountFacets() == 0)
01331     return false;
01332 
01333   // calc each facet
01334   float fMinDist = FLOAT_MAX;
01335   unsigned long ulInd   = ULONG_MAX;
01336   MeshFacetIterator pF(_rclMesh);
01337   for (pF.Init(); pF.More(); pF.Next())
01338   {
01339     float fDist = pF->DistanceToPoint(rclPt);
01340     if (fDist < fMinDist)
01341     {
01342       fMinDist = fDist;
01343       ulInd    = pF.Position();
01344     }
01345   }
01346 
01347   MeshGeomFacet rclSFacet = _rclMesh.GetFacet(ulInd);
01348   rclSFacet.DistanceToPoint(rclPt, rclResPoint);
01349   rclResFacetIndex = ulInd;
01350 
01351   return true;
01352 }
01353 
01354 bool MeshAlgorithm::NearestPointFromPoint (const Base::Vector3f &rclPt, const MeshFacetGrid& rclGrid, unsigned long &rclResFacetIndex, Base::Vector3f &rclResPoint) const
01355 {
01356   unsigned long ulInd = rclGrid.SearchNearestFromPoint(rclPt);
01357 
01358   if (ulInd == ULONG_MAX)
01359   {
01360     return false;
01361   }
01362 
01363   MeshGeomFacet rclSFacet = _rclMesh.GetFacet(ulInd);
01364   rclSFacet.DistanceToPoint(rclPt, rclResPoint);
01365   rclResFacetIndex = ulInd;
01366 
01367   return true;
01368 }
01369 
01370 bool MeshAlgorithm::NearestPointFromPoint (const Base::Vector3f &rclPt, const MeshFacetGrid& rclGrid, float fMaxSearchArea,
01371                                            unsigned long &rclResFacetIndex, Base::Vector3f &rclResPoint) const
01372 {
01373   unsigned long ulInd = rclGrid.SearchNearestFromPoint(rclPt, fMaxSearchArea);
01374 
01375   if (ulInd == ULONG_MAX)
01376     return false;  // no facets inside BoundingBox
01377 
01378   MeshGeomFacet rclSFacet = _rclMesh.GetFacet(ulInd);
01379   rclSFacet.DistanceToPoint(rclPt, rclResPoint);
01380   rclResFacetIndex = ulInd;
01381 
01382   return true;
01383 }
01384 
01385 bool MeshAlgorithm::CutWithPlane (const Base::Vector3f &clBase, const Base::Vector3f &clNormal, const MeshFacetGrid &rclGrid,
01386                                   std::list<std::vector<Base::Vector3f> > &rclResult, float fMinEps, bool bConnectPolygons) const
01387 {
01388   std::vector<unsigned long>  aulFacets;  
01389 
01390   // Grid durschsuchen
01391   MeshGridIterator clGridIter(rclGrid);
01392   for (clGridIter.Init(); clGridIter.More(); clGridIter.Next())
01393   {  
01394     // wenn Gridvoxel Ebene schneidet: alle Facets aufnehmen zum Schneiden
01395     if (clGridIter.GetBoundBox().IsCutPlane(clBase, clNormal) == true)
01396       clGridIter.GetElements(aulFacets);
01397   }
01398 
01399   // mehrfach vorhanden Dreiecke entfernen
01400   std::sort(aulFacets.begin(), aulFacets.end());
01401   aulFacets.erase(std::unique(aulFacets.begin(), aulFacets.end()), aulFacets.end());  
01402 
01403   // alle Facets mit Ebene schneiden
01404   std::list<std::pair<Base::Vector3f, Base::Vector3f> > clTempPoly;  // Feld mit Schnittlinien (unsortiert, nicht verkettet)
01405 
01406   for (std::vector<unsigned long>::iterator pF = aulFacets.begin(); pF != aulFacets.end(); pF++)
01407   {
01408     Base::Vector3f  clE1, clE2;
01409     const MeshGeomFacet clF(_rclMesh.GetFacet(*pF));
01410 
01411     // Facet schneiden und Schnittstrecke ablegen
01412     if (clF.IntersectWithPlane(clBase, clNormal, clE1, clE2) == true)
01413       clTempPoly.push_back(std::pair<Base::Vector3f, Base::Vector3f>(clE1, clE2));
01414   }
01415 
01416   if(bConnectPolygons)
01417   {
01418       //std::list<std::pair<Base::Vector3f, Base::Vector3f> > rclTempLines;
01419       std::list<std::pair<Base::Vector3f, Base::Vector3f> > rclResultLines(clTempPoly.begin(),clTempPoly.end());
01420       std::list<std::vector<Base::Vector3f> > tempList;
01421       ConnectLines(clTempPoly, tempList, fMinEps);
01422       ConnectPolygons(tempList, clTempPoly);
01423 
01424       for(std::list<std::pair<Base::Vector3f, Base::Vector3f> >::iterator iter = clTempPoly.begin(); iter != clTempPoly.end(); iter++)
01425       {
01426         rclResultLines.push_front(*iter);
01427       }
01428 
01429       return ConnectLines(rclResultLines, rclResult, fMinEps);
01430   }
01431 
01432   return ConnectLines(clTempPoly, rclResult, fMinEps);
01433 }
01434 
01435 bool MeshAlgorithm::ConnectLines (std::list<std::pair<Base::Vector3f, Base::Vector3f> > &rclLines,
01436                                   std::list<std::vector<Base::Vector3f> > &rclPolylines, float fMinEps) const
01437 {
01438   typedef std::list<std::pair<Base::Vector3f, Base::Vector3f> >::iterator  TCIter;
01439 
01440   // square search radius
01441   // const float fMinEps = 1.0e-2f; // := 10 mirometer distance
01442   fMinEps = fMinEps * fMinEps;  
01443 
01444   // remove all lines which distance is smaller than epsilon
01445   std::list<TCIter> _clToDelete;
01446   float fToDelDist = fMinEps / 10.0f;
01447   for (TCIter pF = rclLines.begin(); pF != rclLines.end(); pF++)
01448   {
01449     if (Base::DistanceP2(pF->first, pF->second) < fToDelDist)
01450       _clToDelete.push_back(pF);
01451   }
01452   for (std::list<TCIter>::iterator pI = _clToDelete.begin(); pI != _clToDelete.end(); pI++)
01453     rclLines.erase(*pI);
01454 
01455 
01456   while (rclLines.size() > 0)
01457   {
01458     TCIter pF;
01459     
01460     // new polyline
01461     std::list<Base::Vector3f> clPoly;
01462 
01463     // add first line and delete from the list
01464     Base::Vector3f clFront = rclLines.begin()->first;  // current start point of the polyline
01465     Base::Vector3f clEnd   = rclLines.begin()->second; // current end point of the polyline
01466     clPoly.push_back(clFront);
01467     clPoly.push_back(clEnd);
01468     rclLines.erase(rclLines.begin());
01469 
01470     // search for the next line on the begin/end of the polyline and add it
01471     TCIter pFront, pEnd;
01472     do
01473     {
01474       float  fFrontMin = fMinEps, fEndMin = fMinEps;
01475       bool   bFrontFirst=false, bEndFirst=false;
01476 
01477       pFront = rclLines.end();
01478       pEnd   = rclLines.end();
01479       for (pF = rclLines.begin(); pF != rclLines.end(); pF++)
01480       {
01481         if (Base::DistanceP2(clFront, pF->first) < fFrontMin) 
01482         {
01483           fFrontMin   = Base::DistanceP2(clFront, pF->first);
01484           pFront      = pF;
01485           bFrontFirst = true;
01486         }
01487         else if (Base::DistanceP2(clEnd, pF->first) < fEndMin)
01488         {
01489           fEndMin     = Base::DistanceP2(clEnd, pF->first);
01490           pEnd        = pF;
01491           bEndFirst   = true;
01492         }
01493         else if (Base::DistanceP2(clFront, pF->second) < fFrontMin)
01494         {
01495           fFrontMin   = Base::DistanceP2(clFront, pF->second);
01496           pFront      = pF;
01497           bFrontFirst = false;
01498         }
01499         else if (Base::DistanceP2(clEnd, pF->second) < fEndMin)
01500         {
01501           fEndMin     = Base::DistanceP2(clEnd, pF->second);
01502           pEnd        = pF;
01503           bEndFirst   = false;
01504         }        
01505       }
01506 
01507       if (pFront != rclLines.end())
01508       {
01509         if (bFrontFirst == true)
01510         {
01511           clPoly.push_front(pFront->second);
01512           clFront = pFront->second;
01513         }
01514         else
01515         {
01516           clPoly.push_front(pFront->first);
01517           clFront = pFront->first;
01518         }
01519         rclLines.erase(pFront);
01520       }
01521 
01522       if (pEnd != rclLines.end())
01523       {
01524         if (bEndFirst == true)
01525         {
01526           clPoly.push_back(pEnd->second);
01527           clEnd = pEnd->second;
01528         }
01529         else
01530         {
01531           clPoly.push_back(pEnd->first);
01532           clEnd = pEnd->first;
01533         }
01534         rclLines.erase(pEnd);
01535       }
01536     }
01537     while ((pFront != rclLines.end()) || (pEnd != rclLines.end()));
01538 
01539     rclPolylines.push_back(std::vector<Base::Vector3f>(clPoly.begin(), clPoly.end()));
01540   }  
01541   
01542   // remove all polylines with too few length
01543   typedef std::list<std::vector<Base::Vector3f> >::iterator TPIter;
01544   std::list<TPIter> _clPolyToDelete;
01545   for (TPIter pJ = rclPolylines.begin(); pJ != rclPolylines.end(); pJ++)
01546   {
01547     if (pJ->size() == 2) // only one line segment
01548     {
01549       if (Base::DistanceP2(*pJ->begin(), *(pJ->begin() + 1)) <= fMinEps)
01550         _clPolyToDelete.push_back(pJ);
01551     }
01552   }
01553   for (std::list<TPIter>::iterator pK = _clPolyToDelete.begin(); pK != _clPolyToDelete.end(); pK++)
01554     rclPolylines.erase(*pK);
01555 
01556   return true;
01557 }
01558 
01559 bool MeshAlgorithm::ConnectPolygons(std::list<std::vector<Base::Vector3f> > &clPolyList,
01560                                     std::list<std::pair<Base::Vector3f, Base::Vector3f> > &rclLines) const
01561 {
01562 
01563   for(std::list< std::vector<Base::Vector3f> >::iterator OutIter = clPolyList.begin(); OutIter != clPolyList.end() ; OutIter++)
01564   {
01565     std::pair<Base::Vector3f,Base::Vector3f> currentSort;
01566     float fDist = Base::Distance(OutIter->front(),OutIter->back());
01567     currentSort.first = OutIter->front();
01568     currentSort.second = OutIter->back();
01569 
01570     for(std::list< std::vector<Base::Vector3f> >::iterator InnerIter = clPolyList.begin(); InnerIter != clPolyList.end(); InnerIter++)
01571     {
01572       if(OutIter == InnerIter) continue;
01573 
01574       if(Base::Distance(OutIter->front(),InnerIter->front()) < fDist)
01575       {
01576         currentSort.second = InnerIter->front();
01577         fDist = Base::Distance(OutIter->front(),InnerIter->front());
01578       }
01579       if(Base::Distance(OutIter->front(),InnerIter->back()) < fDist)
01580       {
01581         currentSort.second = InnerIter->back();
01582         fDist = Base::Distance(OutIter->front(),InnerIter->back());
01583       }
01584     }
01585 
01586     rclLines.push_front(currentSort);
01587 
01588   }
01589 
01590   return true;
01591 }
01592 
01593 void MeshAlgorithm::GetFacetsFromPlane (const MeshFacetGrid &rclGrid, const Base::Vector3f& clNormal, float d, const Base::Vector3f &rclLeft,
01594                                         const Base::Vector3f &rclRight, std::vector<unsigned long> &rclRes) const
01595 {
01596   std::vector<unsigned long> aulFacets;
01597 
01598   Base::Vector3f clBase = d * clNormal;
01599 
01600   Base::Vector3f clPtNormal(rclLeft - rclRight);
01601   clPtNormal.Normalize();
01602 
01603   // search grid 
01604   MeshGridIterator clGridIter(rclGrid);
01605   for (clGridIter.Init(); clGridIter.More(); clGridIter.Next())
01606   {  
01607     // add facets from grid if the plane if cut the grid-voxel
01608     if (clGridIter.GetBoundBox().IsCutPlane(clBase, clNormal) == true)
01609       clGridIter.GetElements(aulFacets);
01610   }
01611 
01612   // testing facet against planes
01613   for (std::vector<unsigned long>::iterator pI = aulFacets.begin(); pI != aulFacets.end(); pI++)
01614   {
01615     MeshGeomFacet clSFacet = _rclMesh.GetFacet(*pI);
01616     if (clSFacet.IntersectWithPlane(clBase, clNormal) == true)
01617     {
01618       bool bInner = false;
01619       for (int i = 0; (i < 3) && (bInner == false); i++)
01620       {
01621         Base::Vector3f clPt = clSFacet._aclPoints[i];
01622         if ((clPt.DistanceToPlane(rclLeft, clPtNormal) <= 0.0f) && (clPt.DistanceToPlane(rclRight, clPtNormal) >= 0.0f))
01623           bInner = true;
01624       }
01625 
01626       if (bInner == true)
01627         rclRes.push_back(*pI);
01628     }
01629   }
01630 }
01631 
01632 void MeshAlgorithm::PointsFromFacetsIndices (const std::vector<unsigned long> &rvecIndices, std::vector<Base::Vector3f> &rvecPoints) const
01633 {
01634   const MeshFacetArray &rclFAry = _rclMesh._aclFacetArray;
01635   const MeshPointArray &rclPAry = _rclMesh._aclPointArray;
01636 
01637   std::set<unsigned long> setPoints;
01638 
01639   for (std::vector<unsigned long>::const_iterator itI = rvecIndices.begin(); itI != rvecIndices.end(); itI++)
01640   {
01641     for (int i = 0; i < 3; i++)
01642       setPoints.insert(rclFAry[*itI]._aulPoints[i]);
01643   }
01644 
01645   rvecPoints.clear();
01646   for (std::set<unsigned long>::iterator itP = setPoints.begin(); itP != setPoints.end(); itP++)
01647     rvecPoints.push_back(rclPAry[*itP]);
01648 }
01649 
01650 bool MeshAlgorithm::Distance (const Base::Vector3f &rclPt, unsigned long ulFacetIdx, float fMaxDistance, float &rfDistance) const
01651 {
01652   const MeshFacetArray &rclFAry = _rclMesh._aclFacetArray;
01653   const MeshPointArray &rclPAry = _rclMesh._aclPointArray;
01654   const unsigned long *pulIdx = rclFAry[ulFacetIdx]._aulPoints;
01655 
01656   BoundBox3f clBB;
01657   clBB &= rclPAry[*(pulIdx++)];
01658   clBB &= rclPAry[*(pulIdx++)];
01659   clBB &= rclPAry[*pulIdx];
01660   clBB.Enlarge(fMaxDistance);
01661 
01662   if (clBB.IsInBox(rclPt) == false)
01663     return false;
01664 
01665   rfDistance = _rclMesh.GetFacet(ulFacetIdx).DistanceToPoint(rclPt);
01666 
01667   return rfDistance < fMaxDistance;
01668 }
01669 
01670 float MeshAlgorithm::CalculateMinimumGridLength(float fLength, const Base::BoundBox3f& rBBox, unsigned long maxElements) const
01671 {
01672   // Max. limit of grid elements
01673   float fMaxGridElements=(float)maxElements;
01674 
01675   // estimate the minimum allowed grid length
01676   float fMinGridLen = (float)pow((rBBox.LengthX()*rBBox.LengthY()*rBBox.LengthZ()/fMaxGridElements), 0.3333f);
01677   return std::max<float>(fMinGridLen, fLength);
01678 }
01679 
01680 // ----------------------------------------------------
01681 
01682 void MeshRefPointToFacets::Rebuild (void)
01683 {
01684     _map.clear();
01685 
01686     const MeshPointArray& rPoints = _rclMesh.GetPoints();
01687     const MeshFacetArray& rFacets = _rclMesh.GetFacets();
01688     _map.resize(rPoints.size());
01689 
01690     MeshFacetArray::_TConstIterator pFBegin = rFacets.begin();
01691     for (MeshFacetArray::_TConstIterator pFIter = rFacets.begin(); pFIter != rFacets.end(); ++pFIter) {
01692         _map[pFIter->_aulPoints[0]].insert(pFIter - pFBegin);
01693         _map[pFIter->_aulPoints[1]].insert(pFIter - pFBegin);
01694         _map[pFIter->_aulPoints[2]].insert(pFIter - pFBegin);
01695     }
01696 }
01697 
01698 Base::Vector3f MeshRefPointToFacets::GetNormal(unsigned long pos) const
01699 {
01700     const std::set<unsigned long>& n = _map[pos];
01701     Base::Vector3f normal;
01702     MeshGeomFacet f;
01703     for (std::set<unsigned long>::const_iterator it = n.begin(); it != n.end(); ++it) {
01704         f = _rclMesh.GetFacet(*it);
01705         normal += f.Area() * f.GetNormal();
01706     }
01707 
01708     normal.Normalize();
01709     return normal;
01710 }
01711 
01712 std::set<unsigned long> MeshRefPointToFacets::NeighbourPoints(const std::vector<unsigned long>& pt, int level) const
01713 {
01714     std::set<unsigned long> cp,nb,lp;
01715     cp.insert(pt.begin(), pt.end());
01716     lp.insert(pt.begin(), pt.end());
01717     MeshFacetArray::_TConstIterator f_it = _rclMesh.GetFacets().begin();
01718     for (int i=0; i < level; i++) {
01719         std::set<unsigned long> cur;
01720         for (std::set<unsigned long>::iterator it = lp.begin(); it != lp.end(); ++it) {
01721             const std::set<unsigned long>& ft = (*this)[*it];
01722             for (std::set<unsigned long>::const_iterator jt = ft.begin(); jt != ft.end(); ++jt) {
01723                 for (int j = 0; j < 3; j++) {
01724                     unsigned long index = f_it[*jt]._aulPoints[j];
01725                     if (cp.find(index) == cp.end() && nb.find(index) == nb.end()) {
01726                         nb.insert(index);
01727                         cur.insert(index);
01728                     }
01729                 }
01730             }
01731         }
01732 
01733         lp = cur;
01734         if (lp.empty())
01735             break;
01736     }
01737     return nb;
01738 }
01739 
01740 // ermittelt alle Nachbarn zum Facet deren Schwerpunkt unterhalb der mpx. Distanz befindet. 
01741 // Facet deren VISIT-Flag gesetzt ist werden nicht beruecksichtig. 
01743 void MeshRefPointToFacets::Neighbours (unsigned long ulFacetInd, float fMpxDist, std::vector<MeshFacetArray::_TConstIterator> &rclNb)
01744 {
01745     rclNb.clear();
01746     Base::Vector3f  clCenter = _rclMesh.GetFacet(ulFacetInd).GetGravityPoint();
01747 
01748     const MeshFacetArray& rFacets = _rclMesh.GetFacets();
01749     SearchNeighbours(rFacets.begin() + ulFacetInd, clCenter, fMpxDist * fMpxDist, rclNb);
01750 
01751     for (std::vector<MeshFacetArray::_TConstIterator>::iterator i = rclNb.begin(); i != rclNb.end(); i++)
01752         (*i)->ResetFlag(MeshFacet::VISIT);  
01753 }
01754 
01756 void MeshRefPointToFacets::SearchNeighbours(MeshFacetArray::_TConstIterator f_it, const Base::Vector3f &rclCenter, float fMpxDist, std::vector<MeshFacetArray::_TConstIterator> &rclNb)
01757 {
01758     if (f_it->IsFlag(MeshFacet::VISIT) == true)
01759         return;
01760 
01761     if (Base::DistanceP2(rclCenter, _rclMesh.GetFacet(*f_it).GetGravityPoint()) > fMpxDist)
01762         return;
01763 
01764     rclNb.push_back(f_it);
01765     f_it->SetFlag(MeshFacet::VISIT);
01766 
01767     MeshFacetArray::_TConstIterator f_beg = _rclMesh.GetFacets().begin();
01768 
01769     for (int i = 0; i < 3; i++) {
01770         const std::set<unsigned long> &f = (*this)[f_it->_aulPoints[i]];
01771 
01772         for (std::set<unsigned long>::const_iterator j = f.begin(); j != f.end(); ++j) {
01773             SearchNeighbours(f_beg+*j, rclCenter, fMpxDist, rclNb);
01774         }
01775     }
01776 }
01777 
01778 MeshFacetArray::_TConstIterator
01779 MeshRefPointToFacets::getFacet (unsigned long index) const
01780 {
01781     return _rclMesh.GetFacets().begin() + index;
01782 }
01783 
01784 const std::set<unsigned long>&
01785 MeshRefPointToFacets::operator[] (unsigned long pos) const
01786 {
01787     return _map[pos];
01788 }
01789 
01790 //----------------------------------------------------------------------------
01791 
01792 void MeshRefFacetToFacets::Rebuild (void)
01793 {
01794     _map.clear();
01795 
01796     const MeshFacetArray& rFacets = _rclMesh.GetFacets();
01797     _map.resize(rFacets.size());
01798 
01799     MeshRefPointToFacets  vertexFace(_rclMesh);
01800     MeshFacetArray::_TConstIterator pFBegin = rFacets.begin();
01801     for (MeshFacetArray::_TConstIterator pFIter = pFBegin; pFIter != rFacets.end(); ++pFIter) {
01802         for (int i = 0; i < 3; i++) {
01803             const std::set<unsigned long>& faces = vertexFace[pFIter->_aulPoints[i]];
01804             for (std::set<unsigned long>::const_iterator it = faces.begin(); it != faces.end(); ++it)
01805                 _map[pFIter - pFBegin].insert(*it);
01806         }
01807     }
01808 }
01809 
01810 const std::set<unsigned long>&
01811 MeshRefFacetToFacets::operator[] (unsigned long pos) const
01812 {
01813     return _map[pos];
01814 }
01815 
01816 //----------------------------------------------------------------------------
01817 
01818 void MeshRefPointToPoints::Rebuild (void)
01819 {
01820     _map.clear();
01821 
01822     const MeshPointArray& rPoints = _rclMesh.GetPoints();
01823     _map.resize(rPoints.size());
01824 
01825     const MeshFacetArray& rFacets = _rclMesh.GetFacets();
01826     for (MeshFacetArray::_TConstIterator pFIter = rFacets.begin(); pFIter != rFacets.end(); ++pFIter) {
01827         unsigned long ulP0 = pFIter->_aulPoints[0];
01828         unsigned long ulP1 = pFIter->_aulPoints[1];
01829         unsigned long ulP2 = pFIter->_aulPoints[2];
01830 
01831         _map[ulP0].insert(ulP1);
01832         _map[ulP0].insert(ulP2);
01833         _map[ulP1].insert(ulP0);
01834         _map[ulP1].insert(ulP2);
01835         _map[ulP2].insert(ulP0);
01836         _map[ulP2].insert(ulP1);
01837     }
01838 }
01839 
01840 Base::Vector3f MeshRefPointToPoints::GetNormal(unsigned long pos) const
01841 {
01842     const MeshPointArray& rPoints = _rclMesh.GetPoints();
01843     MeshCore::PlaneFit pf;
01844     pf.AddPoint(rPoints[pos]);
01845     MeshCore::MeshPoint center = rPoints[pos];
01846     const std::set<unsigned long>& cv = _map[pos];
01847     for (std::set<unsigned long>::const_iterator cv_it = cv.begin(); cv_it !=cv.end(); ++cv_it) {
01848         pf.AddPoint(rPoints[*cv_it]);
01849         center += rPoints[*cv_it];
01850     }
01851 
01852     pf.Fit();
01853 
01854     Base::Vector3f normal = pf.GetNormal();
01855     normal.Normalize();
01856     return normal;
01857 }
01858 
01859 float MeshRefPointToPoints::GetAverageEdgeLength(unsigned long index) const
01860 {
01861     const MeshPointArray& rPoints = _rclMesh.GetPoints();
01862     float len=0.0f;
01863     const std::set<unsigned long>& n = (*this)[index];
01864     const Base::Vector3f& p = rPoints[index];
01865     for (std::set<unsigned long>::const_iterator it = n.begin(); it != n.end(); ++it) {
01866         len += Base::Distance(p, rPoints[*it]);
01867     }
01868     return (len/n.size());
01869 }
01870 
01871 const std::set<unsigned long>&
01872 MeshRefPointToPoints::operator[] (unsigned long pos) const
01873 {
01874     return _map[pos];
01875 }
01876 
01877 //----------------------------------------------------------------------------
01878 
01879 void MeshRefEdgeToFacets::Rebuild (void)
01880 {
01881     _map.clear();
01882 
01883     const MeshFacetArray& rFacets = _rclMesh.GetFacets();
01884     unsigned long index = 0;
01885     for (MeshFacetArray::_TConstIterator it = rFacets.begin(); it != rFacets.end(); ++it, ++index) {
01886         for (int i=0; i<3; i++) {
01887             MeshEdge e;
01888             e.first = it->_aulPoints[i];
01889             e.second = it->_aulPoints[(i+1)%3];
01890             std::map<MeshEdge, MeshFacetPair, EdgeOrder>::iterator jt =
01891                 _map.find(e);
01892             if (jt == _map.end()) {
01893                 _map[e].first = index;
01894                 _map[e].second = ULONG_MAX;
01895             }
01896             else {
01897                 _map[e].second = index;
01898             }
01899         }
01900     }
01901 }
01902 
01903 const std::pair<unsigned long, unsigned long>&
01904 MeshRefEdgeToFacets::operator[] (const MeshEdge& edge) const
01905 {
01906     return _map.find(edge)->second;
01907 }
01908 
01909 //----------------------------------------------------------------------------
01910 
01911 void MeshRefNormalToPoints::Rebuild (void)
01912 {
01913     _norm.clear();
01914 
01915     const MeshPointArray& rPoints = _rclMesh.GetPoints();
01916     _norm.resize(rPoints.size());
01917 
01918     const MeshFacetArray& rFacets = _rclMesh.GetFacets();
01919     for (MeshFacetArray::_TConstIterator pF = rFacets.begin(); pF != rFacets.end(); ++pF) {
01920         const MeshPoint &p0 = rPoints[pF->_aulPoints[0]];
01921         const MeshPoint &p1 = rPoints[pF->_aulPoints[1]];
01922         const MeshPoint &p2 = rPoints[pF->_aulPoints[2]];
01923         float l2p01 = Base::DistanceP2(p0,p1);
01924         float l2p12 = Base::DistanceP2(p1,p2);
01925         float l2p20 = Base::DistanceP2(p2,p0);
01926 
01927         Base::Vector3f facenormal = _rclMesh.GetFacet(*pF).GetNormal();
01928         _norm[pF->_aulPoints[0]] += facenormal * (1.0f / (l2p01 * l2p20));
01929         _norm[pF->_aulPoints[1]] += facenormal * (1.0f / (l2p12 * l2p01));
01930         _norm[pF->_aulPoints[2]] += facenormal * (1.0f / (l2p20 * l2p12));
01931     }
01932     for (std::vector<Base::Vector3f>::iterator it = _norm.begin(); it != _norm.end(); ++it)
01933         it->Normalize();
01934 }
01935 
01936 const Base::Vector3f&
01937 MeshRefNormalToPoints::operator[] (unsigned long pos) const
01938 {
01939     return _norm[pos];
01940 }

Generated on Wed Nov 23 18:59:54 2011 for FreeCAD by  doxygen 1.6.1