Fix MA construction

src/SMESHUtils/SMESH_MAT2d.cxx

@@ -47,6 +47,7 @@
 #include <TopoDS_Wire.hxx>
 
 #ifdef _DEBUG_
+//#define _MYDEBUG_
 #include "SMESH_File.hxx"
 #include "SMESH_Comment.hxx"
 #endif
@@ -66,8 +67,8 @@ namespace
 
   struct InPoint
   {
-    int _a, _b;
+    int _a, _b;    // coordinates
-    double _param;
+    double _param; // param on EDGE
     InPoint(int x, int y, double param) : _a(x), _b(y), _param(param) {}
     InPoint() : _a(0), _b(0), _param(0) {}
 
@@ -76,6 +77,8 @@ namespace
 
     size_t index( const vector< InPoint >& inPoints ) const { return this - &inPoints[0]; }
     bool operator==( const InPoint& other ) const { return _a == other._a && _b == other._b; }
+    bool operator==( const TVDVertex* v ) const { return ( Abs( _a - v->x() ) < 1. &&
+                                                           Abs( _b - v->y() ) < 1. ); }
   };
   // -------------------------------------------------------------------------------------
 
@@ -105,10 +108,12 @@ namespace
   // check  if a TVDEdge begins at my end or ends at my start
   inline bool InSegment::isConnected( const TVDEdge* edge )
   {
-    return ((Abs( edge->vertex0()->x() - _p1->_a ) < 1.&&
+    return (( edge->vertex0() && edge->vertex1() )
-             Abs( edge->vertex0()->y() - _p1->_b ) < 1.  ) ||
+            &&
-            (Abs( edge->vertex1()->x() - _p0->_a ) < 1.&&
+            ((Abs( edge->vertex0()->x() - _p1->_a ) < 1.&&
-             Abs( edge->vertex1()->y() - _p0->_b ) < 1.  ));
+              Abs( edge->vertex0()->y() - _p1->_b ) < 1.  ) ||
+             (Abs( edge->vertex1()->x() - _p0->_a ) < 1.&&
+              Abs( edge->vertex1()->y() - _p0->_b ) < 1.  )));
   }
 
   // check if a MA TVDEdge is outside of a domain
@@ -147,7 +152,7 @@ namespace
   // }
 
   // -------------------------------------------------------------------------------------
-#ifdef _DEBUG_
+#ifdef _MYDEBUG_
   // writes segments into a txt file readable by voronoi_visualizer
   void inSegmentsToFile( vector< InSegment>& inSegments)
   {
@@ -155,6 +160,7 @@ namespace
       return;
     const char* fileName = "/misc/dn25/salome/eap/salome/misc/Code/C++/MAdebug.txt";
     SMESH_File file(fileName, false );
+    file.remove();
     file.openForWriting();
     SMESH_Comment text;
     text << "0\n"; // nb points
@@ -180,7 +186,7 @@ namespace
     if ( !edge->vertex1() )
       cout << ") -> ( INF, INF";
     else
-      cout << ") -> (" << edge->vertex1()->x() << ", " << edge->vertex1()->y();
+      cout << ") -> ( " << edge->vertex1()->x() << ", " << edge->vertex1()->y();
     cout << ")\t cell=" << edge->cell()
          << " iBnd=" << edge->color()
          << " twin=" << edge->twin()
@@ -253,6 +259,7 @@ namespace boost {
 namespace
 {
   const int theNoBrachID = 0; // std::numeric_limits<int>::max();
+  double theScale[2]; // scale used in bndSegsToMesh()
 
   // -------------------------------------------------------------------------------------
   /*!
@@ -298,7 +305,7 @@ namespace
       if ( _edge ) // pass branch to an opposite BndSeg
       {
         size_t oppSegIndex = SMESH_MAT2d::Branch::getBndSegment( _edge->twin() );
-        if ( oppSegIndex < bndSegs.size() /*&& bndSegs[ oppSegIndex ]._branchID == theNoBrachID*/ )
+        if ( oppSegIndex < bndSegs.size() && bndSegs[ oppSegIndex ]._branchID == theNoBrachID )
           bndSegs[ oppSegIndex ]._branchID = -branchID;
       }
     }
@@ -349,7 +356,64 @@ namespace
 
   //================================================================================
   /*!
-   * \brief Computes length of of TVDEdge
+   * \brief debug: to visually check found MA edges
+   */
+  //================================================================================
+
+  void bndSegsToMesh( const vector< BndSeg >& bndSegs )
+  {
+#ifdef _MYDEBUG_
+    if ( !getenv("bndSegsToMesh")) return;
+    map< const TVDVertex *, int > v2Node;
+    map< const TVDVertex *, int >::iterator v2n;
+    set< const TVDEdge* > addedEdges;
+
+    const char* fileName = "/misc/dn25/salome/eap/salome/misc/Code/C++/MAedges.py";
+    SMESH_File file(fileName, false );
+    file.remove();
+    file.openForWriting();
+    SMESH_Comment text;
+    text << "import salome, SMESH\n";
+    text << "salome.salome_init()\n";
+    text << "from salome.smesh import smeshBuilder\n";
+    text << "smesh = smeshBuilder.New(salome.myStudy)\n";
+    text << "m=smesh.Mesh()\n";
+    for ( size_t i = 0; i < bndSegs.size(); ++i )
+    {
+      if ( !bndSegs[i]._edge )
+        text << "# " << i << " NULL edge";
+      else if ( !bndSegs[i]._edge->vertex0() ||
+                !bndSegs[i]._edge->vertex1() )
+        text << "# " << i << " INFINITE edge";
+      else if ( addedEdges.insert( bndSegs[i]._edge ).second &&
+                addedEdges.insert( bndSegs[i]._edge->twin() ).second )
+      {
+        v2n = v2Node.insert( make_pair( bndSegs[i]._edge->vertex0(), v2Node.size() + 1 )).first;
+        int n0 = v2n->second;
+        if ( n0 == v2Node.size() )
+          text << "n" << n0 << " = m.AddNode( "
+               << bndSegs[i]._edge->vertex0()->x() / theScale[0] << ", "
+               << bndSegs[i]._edge->vertex0()->y() / theScale[1] << ", 0 )\n";
+
+        v2n = v2Node.insert( make_pair( bndSegs[i]._edge->vertex1(), v2Node.size() + 1 )).first;
+        int n1 = v2n->second;
+        if ( n1 == v2Node.size() )
+          text << "n" << n1 << " = m.AddNode( "
+               << bndSegs[i]._edge->vertex1()->x() / theScale[0] << ", "
+               << bndSegs[i]._edge->vertex1()->y() / theScale[1] << ", 0 )\n";
+
+        text << "e" << i << " = m.AddEdge([ n" << n0 << ", n" << n1 << " ])\n";
+      }
+    }
+    text << "\n";
+    file.write( text.c_str(), text.size() );
+    cout << "Write " << fileName << endl;
+#endif
+  }
+
+  //================================================================================
+  /*!
+   * \brief Computes length of a TVDEdge
    */
   //================================================================================
 
@@ -459,24 +523,25 @@ namespace
       gp_Pnt p, pPrev;
       if ( !c3d.IsNull() )
         pPrev = c3d->Value( f );
-      for ( int i = 2; i <= discret.NbPoints(); i++ ) // skip the 1st point
+      if ( discret.NbPoints() > 2 )
-      {
+        for ( int i = 2; i <= discret.NbPoints(); i++ ) // skip the 1st point
-        double u = discret.Parameter(i);
-        if ( !c3d.IsNull() )
         {
-          p = c3d->Value( u );
+          double u = discret.Parameter(i);
-          int nbDiv = int( p.Distance( pPrev ) / minSegLen / theDiscrCoef );
+          if ( !c3d.IsNull() )
-          double dU = ( u - points.back()._u ) / nbDiv;
-          for ( int iD = 1; iD < nbDiv; ++iD )
           {
-            double uD = points.back()._u + dU;
+            p = c3d->Value( u );
-            points.push_back( UVU( c2d->Value( uD ), uD ));
+            int nbDiv = int( p.Distance( pPrev ) / minSegLen / theDiscrCoef );
+            double dU = ( u - points.back()._u ) / nbDiv;
+            for ( int iD = 1; iD < nbDiv; ++iD )
+            {
+              double uD = points.back()._u + dU;
+              points.push_back( UVU( c2d->Value( uD ), uD ));
+            }
+            pPrev = p;
           }
-          pPrev = p;
+          points.push_back( UVU( c2d->Value( u ), u ));
+          uvBox.Add( points.back()._uv );
         }
-        points.push_back( UVU( c2d->Value( u ), u ));
-        uvBox.Add( points.back()._uv );
-      }
       // if ( !c3d.IsNull() )
       // {
       //   vector<double> params;
@@ -592,9 +657,52 @@ namespace
         }
       }
     }
+    // debug
+    theScale[0] = scale[0];
+    theScale[1] = scale[1];
+
     return true;
   }
 
+  //================================================================================
+  /*!
+   * \brief Update a branch joined to another one
+   */
+  //================================================================================
+
+  void updateJoinedBranch( vector< const TVDEdge* > & branchEdges,
+                           const size_t               newID,
+                           vector< BndSeg > &         bndSegs,
+                           const bool                 reverse)
+  {
+    if ( reverse )
+    {
+      for ( size_t i = 0; i < branchEdges.size(); ++i )
+      {
+        size_t seg1 = SMESH_MAT2d::Branch::getBndSegment( branchEdges[i] );
+        size_t seg2 = SMESH_MAT2d::Branch::getBndSegment( branchEdges[i]->twin() );
+        bndSegs[ seg1 ]._branchID /= bndSegs[ seg1 ].branchID();
+        bndSegs[ seg2 ]._branchID /= bndSegs[ seg2 ].branchID();
+        bndSegs[ seg1 ]._branchID *= -newID;
+        bndSegs[ seg2 ]._branchID *= -newID;
+        branchEdges[i] = branchEdges[i]->twin();
+      }
+      std::reverse( branchEdges.begin(), branchEdges.end() );
+    }
+    else
+    {
+      for ( size_t i = 0; i < branchEdges.size(); ++i )
+      {
+        size_t seg1 = SMESH_MAT2d::Branch::getBndSegment( branchEdges[i] );
+        size_t seg2 = SMESH_MAT2d::Branch::getBndSegment( branchEdges[i]->twin() );
+        bndSegs[ seg1 ]._branchID /= bndSegs[ seg1 ].branchID();
+        bndSegs[ seg2 ]._branchID /= bndSegs[ seg2 ].branchID();
+        bndSegs[ seg1 ]._branchID *= newID;
+        bndSegs[ seg2 ]._branchID *= newID;
+      }
+    }
+  }
+
   //================================================================================
   /*!
    * \brief Create MA branches and FACE boundary data
@@ -608,6 +716,7 @@ namespace
   //================================================================================
 
   void makeMA( const TVD&                               vd,
+               const bool                               ignoreCorners,
                vector< InPoint >&                       inPoints,
                vector< InSegment > &                    inSegments,
                vector< SMESH_MAT2d::Branch >&           branch,
@@ -688,9 +797,12 @@ namespace
           continue;
         inPntChecked[ pInd ] = true;
 
-        const TVDEdge* edge =  // a TVDEdge passing through an end of inSeg
+        const TVDEdge* maE = is2nd ? maEdges.front() : maEdges.back();
-          is2nd ? maEdges.front()->prev() : maEdges.back()->next();
+        if ( inPnt == ( is2nd ? maE->vertex0() : maE->vertex1() ))
-        while ( true )
+          continue;
+        const TVDEdge* edge =  // a secondary TVDEdge connecting inPnt and maE
+          is2nd ? maE->prev() : maE->next();
+        while ( inSeg.isConnected( edge ))
         {
           if ( edge->is_primary() ) break; // this should not happen
           const TVDEdge* edge2 = edge->twin(); // we are in a neighbor cell, add MA edges to inPnt
@@ -732,7 +844,7 @@ namespace
       inPoints[0]._edges.clear();
     }
 
-    // Divide InSegment's into BndSeg's
+    // Divide InSegment's into BndSeg's (so that each BndSeg corresponds to one MA edge)
 
     vector< BndSeg > bndSegs;
     bndSegs.reserve( inSegments.size() * 3 );
@@ -750,25 +862,26 @@ namespace
       inPntChecked[ ip0 ] = false;
 
       // segments of InSegment's
-      size_t nbMaEdges = inSeg._edges.size();
+      const size_t nbMaEdges = inSeg._edges.size();
       switch ( nbMaEdges ) {
       case 0: // "around" circle center
         bndSegs.push_back( BndSeg( &inSeg, 0, inSeg._p1->_param )); break;
       case 1:
         bndSegs.push_back( BndSeg( &inSeg, inSeg._edges.back(), inSeg._p1->_param )); break;
       default:
-        vector< double > len;
+        gp_XY inSegDir( inSeg._p1->_a - inSeg._p0->_a,
-        len.push_back(0);
+                        inSeg._p1->_b - inSeg._p0->_b );
-        for ( e = inSeg._edges.rbegin(); e != inSeg._edges.rend(); ++e )
+        const double inSegLen2 = inSegDir.SquareModulus();
-          len.push_back( len.back() + length( *e ));
-
         e = inSeg._edges.rbegin();
-        for ( size_t l = 1; l < len.size(); ++e, ++l )
+        for ( size_t iE = 1; iE < nbMaEdges; ++e, ++iE )
         {
-          double dl = len[l] / len.back();
+          gp_XY toMA( (*e)->vertex0()->x() - inSeg._p0->_a,
-          double u  = dl * inSeg._p1->_param + ( 1. - dl ) * inSeg._p0->_param;
+                      (*e)->vertex0()->y() - inSeg._p0->_b );
+          double r = toMA * inSegDir / inSegLen2;
+          double u = r * inSeg._p1->_param + ( 1. - r ) * inSeg._p0->_param;
           bndSegs.push_back( BndSeg( &inSeg, *e, u ));
         }
+        bndSegs.push_back( BndSeg( &inSeg, *e, inSeg._p1->_param ));
       }
       // segments around 2nd concave point
       size_t ip1 = inSeg._p1->index( inPoints );
@@ -783,6 +896,8 @@ namespace
     for ( size_t i = 0; i < bndSegs.size(); ++i )
       bndSegs[i].setIndexToEdge( i );
 
+    bndSegsToMesh( bndSegs ); // debug: visually check found MA edges
+
 
     // Find TVDEdge's of Branches and associate them with bndSegs
 
@@ -797,7 +912,7 @@ namespace
     size_t i1st = 0;
     while ( i1st < bndSegs.size() && !bndSegs[i1st].hasOppositeEdge( noEdgeID ))
       ++i1st;
-    bndSegs[i1st].setBranch( branchID, bndSegs ); // set to i-th and the opposite bndSeg
+    bndSegs[i1st].setBranch( branchID, bndSegs ); // set to the i-th and to the opposite bndSeg
     branchEdges[ branchID ].push_back( bndSegs[i1st]._edge );
 
     for ( size_t i = i1st+1; i < bndSegs.size(); ++i )
@@ -824,7 +939,7 @@ namespace
           endType.insert( make_pair( bndSegs[i]._edge->vertex1(),
                                      SMESH_MAT2d::BE_BRANCH_POINT ));
       }
-      bndSegs[i].setBranch( branchID, bndSegs ); // set to i-th and the opposite bndSeg
+      bndSegs[i].setBranch( branchID, bndSegs ); // set to i-th and to the opposite bndSeg
       if ( bndSegs[i].hasOppositeEdge( noEdgeID ))
         branchEdges[ bndSegs[i].branchID() ].push_back( bndSegs[i]._edge );
     }
@@ -858,17 +973,92 @@ namespace
       br2.clear();
     }
 
+    // remove branches ending at BE_ON_VERTEX
+
+    vector<bool> isBranchRemoved( branchEdges.size(), false );
+
+    if ( ignoreCorners && branchEdges.size() > 2 && !branchEdges[2].empty() )
+    {
+      // find branches to remove
+      map< const TVDVertex*, SMESH_MAT2d::BranchEndType >::iterator v2et;
+      for ( size_t iB = 1; iB < branchEdges.size(); ++iB )
+      {
+        if ( branchEdges[iB].empty() )
+          continue;
+        const TVDVertex* v0 = branchEdges[iB][0]->vertex1();
+        const TVDVertex* v1 = branchEdges[iB].back()->vertex0();
+        v2et = endType.find( v0 );
+        if ( v2et != endType.end() && v2et->second == SMESH_MAT2d::BE_ON_VERTEX )
+          isBranchRemoved[ iB ] = true;
+        v2et = endType.find( v1 );
+        if ( v2et != endType.end() && v2et->second == SMESH_MAT2d::BE_ON_VERTEX )
+          isBranchRemoved[ iB ] = true;
+      }
+      // try to join not removed branches into one
+      for ( size_t iB = 1; iB < branchEdges.size(); ++iB )
+      {
+        if ( branchEdges[iB].empty() || isBranchRemoved[iB] )
+          continue;
+        const TVDVertex* v0 = branchEdges[iB][0]->vertex1();
+        const TVDVertex* v1 = branchEdges[iB].back()->vertex0();
+        v2et = endType.find( v0 );
+        if ( v2et == endType.end() || v2et->second != SMESH_MAT2d::BE_BRANCH_POINT )
+          v0 = 0;
+        v2et = endType.find( v1 );
+        if ( v2et == endType.end() || v2et->second != SMESH_MAT2d::BE_BRANCH_POINT )
+          v1 = 0;
+        if ( !v0 && !v1 )
+          continue;
+
+        size_t iBrToJoin = 0;
+        for ( size_t iB2 = 1; iB2 < branchEdges.size(); ++iB2 )
+        {
+          if ( branchEdges[iB2].empty() || isBranchRemoved[iB2] || iB == iB2 )
+            continue;
+          const TVDVertex* v02 = branchEdges[iB2][0]->vertex1();
+          const TVDVertex* v12 = branchEdges[iB2].back()->vertex0();
+          if ( v0 == v02 || v0 == v12 || v1 == v02 || v1 == v12 )
+          {
+            if ( iBrToJoin > 0 )
+            {
+              iBrToJoin = 0;
+              break; // more than 2 not removed branches meat at a TVDVertex
+            }
+            iBrToJoin = iB2;
+          }
+        }
+        if ( iBrToJoin > 0 )
+        {
+          vector<const TVDEdge*>& branch = branchEdges[ iBrToJoin ];
+          const TVDVertex* v02 = branch[0]->vertex1();
+          const TVDVertex* v12 = branch.back()->vertex0();
+          updateJoinedBranch( branch, iB, bndSegs, /*reverse=*/(v0 == v02 || v1 == v12 ));
+          if ( v0 == v02 || v0 == v12 )
+            branchEdges[iB].insert( branchEdges[iB].begin(), branch.begin(), branch.end() );
+          else
+            branchEdges[iB].insert( branchEdges[iB].end(),   branch.begin(), branch.end() );
+          branch.clear();
+        }
+      } // loop on branchEdges
+    } // if ( ignoreCorners )
+
     // associate branchIDs and the input branch vector (arg)
-    vector< const SMESH_MAT2d::Branch* > branchByID( branchEdges.size() );
+    vector< SMESH_MAT2d::Branch* > branchByID( branchEdges.size(), 0 );
     int nbBranches = 0;
     for ( size_t i = 0; i < branchEdges.size(); ++i )
     {
       nbBranches += ( !branchEdges[i].empty() );
     }
     branch.resize( nbBranches );
-    for ( size_t iBr = 0, brID = 0; brID < branchEdges.size(); ++brID )
+    size_t iBr = 0;
+    for ( size_t brID = 1; brID < branchEdges.size(); ++brID ) // 1st - not removed
     {
-      if ( !branchEdges[ brID ].empty() )
+      if ( !branchEdges[ brID ].empty() && !isBranchRemoved[ brID ])
+        branchByID[ brID ] = & branch[ iBr++ ];
+    }
+    for ( size_t brID = 1; brID < branchEdges.size(); ++brID ) // then - removed
+    {
+      if ( !branchEdges[ brID ].empty() && isBranchRemoved[ brID ])
         branchByID[ brID ] = & branch[ iBr++ ];
     }
 
@@ -912,28 +1102,24 @@ namespace
         {
           edgeInd += dInd;
 
-          if ( edgeInd < 0 ||
+          if (( edgeInd < 0 ||
-               edgeInd >= (int) branchEdges[ brID ].size() ||
+                edgeInd >= (int) branchEdges[ brID ].size() ) ||
-               branchEdges[ brID ][ edgeInd ] != bndSegs[ i ]._edge )
+              ( branchEdges[ brID ][ edgeInd ]         != bndSegs[ i ]._edge &&
+                branchEdges[ brID ][ edgeInd ]->twin() != bndSegs[ i ]._edge ))
           {
-            if ( bndSegs[ i ]._branchID < 0 &&
+            if ( bndSegs[ i ]._branchID < 0 )
-                 branchEdges[ brID ].back() == bndSegs[ i ]._edge )
             {
-              edgeInd = branchEdges[ brID ].size() - 1;
               dInd = -1;
+              for ( edgeInd = branchEdges[ brID ].size() - 1; edgeInd > 0; --edgeInd )
+                if ( branchEdges[ brID ][ edgeInd ]->twin() == bndSegs[ i ]._edge )
+                  break;
             }
-            else if ( bndSegs[ i ]._branchID > 0 &&
+            else // bndSegs[ i ]._branchID > 0
-                      branchEdges[ brID ].front() == bndSegs[ i ]._edge )
             {
-              edgeInd = 0;
               dInd = +1;
-            }
-            else
-            {
               for ( edgeInd = 0; edgeInd < branchEdges[ brID ].size(); ++edgeInd )
                 if ( branchEdges[ brID ][ edgeInd ] == bndSegs[ i ]._edge )
                   break;
-              dInd = bndSegs[ i ]._branchID > 0 ? +1 : -1;
             }
           }
         }
@@ -941,15 +1127,10 @@ namespace
         {
           // no MA edge, bndSeg corresponds to an end point of a branch
           if ( bndPoints._maEdges.empty() )
-          {
-            // should not get here according to algo design
             edgeInd = 0;
-          }
           else
-          {
             edgeInd = branchEdges[ brID ].size();
-            dInd = bndSegs[ i ]._branchID > 0 ? +1 : -1;
+          dInd = bndSegs[ i ]._branchID > 0 ? +1 : -1;
-          }
         }
 
         bndPoints._maEdges.push_back( make_pair( br, ( 1 + edgeInd ) * dInd ));
@@ -958,19 +1139,40 @@ namespace
 
       iSeg = iSegEnd;
 
-    } // loop on all bndSegs
+    } // loop on all bndSegs to construct Boundary
 
+    // Initialize branches
 
+    // find a not removed branch
+    size_t iBrNorRemoved = 0;
+    for ( size_t brID = 1; brID < branchEdges.size(); ++brID )
+      if ( !branchEdges[brID].empty() && !isBranchRemoved[brID] )
+      {
+        iBrNorRemoved = brID;
+        break;
+      }
     // fill the branches with MA edges
-    for ( size_t iBr = 0, brID = 0; brID < branchEdges.size(); ++brID )
+    for ( size_t brID = 1; brID < branchEdges.size(); ++brID )
       if ( !branchEdges[brID].empty() )
       {
-        branch[ iBr ].init( branchEdges[brID], & boundary, endType );
+        branchByID[ brID ]->init( branchEdges[brID], & boundary, endType );
-        iBr++;
+      }
+    // mark removed branches
+    for ( size_t brID = 1; brID < branchEdges.size(); ++brID )
+      if ( isBranchRemoved[brID] && iBrNorRemoved > 0 )
+      {
+        SMESH_MAT2d::Branch* branch     = branchByID[ brID ];
+        SMESH_MAT2d::Branch* mainBranch = branchByID[ iBrNorRemoved ];
+        bool is1stBrPnt = ( branch->getEnd(0)->_type == SMESH_MAT2d::BE_BRANCH_POINT );
+        const TVDVertex* branchVextex = 
+          is1stBrPnt ? branch->getEnd(0)->_vertex : branch->getEnd(1)->_vertex;
+        SMESH_MAT2d::BranchPoint bp = mainBranch->getPoint( branchVextex );
+        branch->setRemoved( bp );
       }
     // set branches to branch ends
     for ( size_t i = 0; i < branch.size(); ++i )
-      branch[i].setBranchesToEnds( branch );
+      if ( !branch[i].isRemoved() )
+        branch[i].setBranchesToEnds( branch );
 
     // fill branchPnt arg
     map< const TVDVertex*, const SMESH_MAT2d::BranchEnd* > v2end;
@@ -1013,13 +1215,30 @@ SMESH_MAT2d::MedialAxis::MedialAxis(const TopoDS_Face&                face,
   if ( !makeInputData( face, edges, minSegLen, inPoints, inSegments, _scale ))
     return;
 
-  //inSegmentsToFile( inSegments );
+  inSegmentsToFile( inSegments );
 
   // build voronoi diagram
   construct_voronoi( inSegments.begin(), inSegments.end(), &_vd );
 
   // make MA data
-  makeMA( _vd, inPoints, inSegments, _branch, _branchPnt, _boundary );
+  makeMA( _vd, ignoreCorners, inPoints, inSegments, _branch, _branchPnt, _boundary );
+
+  // count valid branches
+  _nbBranches = _branch.size();
+  for ( size_t i = 0; i < _branch.size(); ++i )
+    if ( _branch[i].isRemoved() )
+      --_nbBranches;
+}
+
+//================================================================================
+/*!
+ * \brief Returns the i-th branch
+ */
+//================================================================================
+
+const SMESH_MAT2d::Branch* SMESH_MAT2d::MedialAxis::getBranch(size_t i) const
+{
+  return i < _nbBranches ? &_branch[i] : 0;
 }
 
 //================================================================================
@@ -1028,16 +1247,16 @@ SMESH_MAT2d::MedialAxis::MedialAxis(const TopoDS_Face&                face,
  */
 //================================================================================
 
-void SMESH_MAT2d::MedialAxis::getPoints( const Branch&         branch,
+void SMESH_MAT2d::MedialAxis::getPoints( const Branch*         branch,
                                          std::vector< gp_XY >& points) const
 {
-  branch.getPoints( points, _scale );
+  branch->getPoints( points, _scale );
 }
 
 //================================================================================
 /*!
  * \brief Returns a BranchPoint corresponding to a given point on a geom EDGE
- *  \param [in] iGeomEdge - index of geom EDGE within a vector passed at construction
+ *  \param [in] iEdge - index of geom EDGE within a vector passed at MA construction
  *  \param [in] u - parameter of the point on EDGE curve
  *  \param [out] p - the found BranchPoint
  *  \return bool - is OK
@@ -1052,7 +1271,7 @@ bool SMESH_MAT2d::Boundary::getBranchPoint( const std::size_t iEdge,
     return false;
 
   const BndPoints& points = _pointsPerEdge[ iEdge ];
-  const bool edgeReverse = ( points._params[0] > points._params.back() );
+  const bool  edgeReverse = ( points._params[0] > points._params.back() );
 
   if ( u < ( edgeReverse ? points._params.back() : points._params[0] ))
     u = edgeReverse ? points._params.back() : points._params[0];
@@ -1071,14 +1290,30 @@ bool SMESH_MAT2d::Boundary::getBranchPoint( const std::size_t iEdge,
     while ( points._params[i  ] > u ) --i;
     while ( points._params[i+1] < u ) ++i;
   }
+
   double edgeParam = ( u - points._params[i] ) / ( points._params[i+1] - points._params[i] );
 
+  if ( !points._maEdges[ i ].second ) // no branch at the EDGE end, look for a closest branch
+  {
+    if ( i < points._maEdges.size() / 2 ) // near 1st point
+    {
+      while ( i < points._maEdges.size()-1 && !points._maEdges[ i ].second )
+        ++i;
+      edgeParam = edgeReverse;
+    }
+    else // near last point
+    {
+      while ( i > 0 && !points._maEdges[ i ].second )
+        --i;
+      edgeParam = !edgeReverse;
+    }
+  }
   const std::pair< const Branch*, int >& maE = points._maEdges[ i ];
   bool maReverse = ( maE.second < 0 );
 
-  p._branch = maE.first;
+  p._branch    = maE.first;
-  p._iEdge  = maE.second - 1; // countered from 1 to store sign
+  p._iEdge     = ( maReverse ? -maE.second : maE.second ) - 1; // countered from 1 to store sign
-  p._edgeParam = maReverse ? ( 1. - edgeParam ) : edgeParam;
+  p._edgeParam = ( maE.first && maReverse ) ? ( 1. - edgeParam ) : edgeParam;
 
   return true;
 }
@@ -1093,16 +1328,36 @@ bool SMESH_MAT2d::Boundary::getBranchPoint( const std::size_t iEdge,
  */
 //================================================================================
 
-bool SMESH_MAT2d::Boundary::IsConcaveSegment( std::size_t iEdge, std::size_t iSeg ) const
+bool SMESH_MAT2d::Boundary::isConcaveSegment( std::size_t iEdge, std::size_t iSeg ) const
 {
   if ( iEdge >= _pointsPerEdge.size() || _pointsPerEdge[iEdge]._params.empty() )
     return false;
 
   const BndPoints& points = _pointsPerEdge[ iEdge ];
-  if ( points._params.size() >= iSeg+1 )
+  if ( points._params.size() <= iSeg+1 )
     return false;
 
-  return Abs( points._params[ iEdge ] - points._params[ iEdge+1 ]) < 1e-20;
+  return Abs( points._params[ iSeg ] - points._params[ iSeg+1 ]) < 1e-20;
+}
+
+//================================================================================
+/*!
+ * \brief Moves (changes _param) a given BoundaryPoint to a closest EDGE end
+ */
+//================================================================================
+
+bool SMESH_MAT2d::Boundary::moveToClosestEdgeEnd( BoundaryPoint& bp ) const
+{
+  if ( bp._edgeIndex >= _pointsPerEdge.size() )
+    return false;
+
+  const BndPoints& points = _pointsPerEdge[ bp._edgeIndex ];
+  if ( bp._param - points._params[0] < points._params.back() - bp._param )
+    bp._param = points._params[0];
+  else 
+    bp._param = points._params.back();
+
+  return true;
 }
 
 //================================================================================
@@ -1179,7 +1434,7 @@ void SMESH_MAT2d::Branch::init( vector<const TVDEdge*>&                maEdges,
 
 //================================================================================
 /*!
- * \brief fill BranchEnd::_branches of its ends
+ * \brief fills BranchEnd::_branches of its ends
  */
 //================================================================================
 
@@ -1197,6 +1452,47 @@ void SMESH_MAT2d::Branch::setBranchesToEnds( const vector< Branch >& branches )
   }
 }
 
+//================================================================================
+/*!
+ * \brief returns a BranchPoint corresponding to a TVDVertex
+ */
+//================================================================================
+
+SMESH_MAT2d::BranchPoint SMESH_MAT2d::Branch::getPoint( const TVDVertex* vertex ) const
+{
+  BranchPoint p;
+  p._branch = this;
+  p._iEdge  = 0;
+
+  if ( vertex == _maEdges[0]->vertex1() )
+  {
+    p._edgeParam = 0;
+  }
+  else
+  {
+    for ( ; p._iEdge < _maEdges.size(); ++p._iEdge )
+      if ( vertex == _maEdges[ p._iEdge ]->vertex0() )
+      {
+        p._edgeParam = _params[ p._iEdge ];
+        break;
+      }
+  }
+  return p;
+}
+
+//================================================================================
+/*!
+ * \brief Sets a proxy point for a removed branch
+ *  \param [in] proxyPoint - a point of another branch to which all points of this
+ *         branch are mapped
+ */
+//================================================================================
+
+void SMESH_MAT2d::Branch::setRemoved( const BranchPoint& proxyPoint )
+{
+  _proxyPoint = proxyPoint;
+}
+
 //================================================================================
 /*!
  * \brief Returns points on two EDGEs, equidistant from a given point of this Branch
@@ -1213,7 +1509,7 @@ bool SMESH_MAT2d::Branch::getBoundaryPoints(double         param,
 {
   if ( param < _params[0] || param > _params.back() )
     return false;
-  
+
   // look for an index of a MA edge by param
   double ip = param * _params.size();
   size_t  i = size_t( Min( int( _maEdges.size()-1), int( ip )));
@@ -1242,8 +1538,13 @@ bool SMESH_MAT2d::Branch::getBoundaryPoints(std::size_t    iMAEdge,
                                             BoundaryPoint& bp1,
                                             BoundaryPoint& bp2 ) const
 {
+  if ( isRemoved() )
+    return _proxyPoint._branch->getBoundaryPoints( _proxyPoint, bp1, bp2 );
+
   if ( iMAEdge > _maEdges.size() )
     return false;
+  if ( iMAEdge == _maEdges.size() )
+    iMAEdge = _maEdges.size() - 1;
 
   size_t iGeom1 = getGeomEdge( _maEdges[ iMAEdge ] );
   size_t iGeom2 = getGeomEdge( _maEdges[ iMAEdge ]->twin() );
@@ -1275,8 +1576,16 @@ bool SMESH_MAT2d::Branch::getBoundaryPoints(const BranchPoint& p,
 
 bool SMESH_MAT2d::Branch::getParameter(const BranchPoint & p, double & u ) const
 {
+  if ( this != p._branch && p._branch )
+    return p._branch->getParameter( p, u );
+
+  if ( isRemoved() )
+    return _proxyPoint._branch->getParameter( _proxyPoint, u );
+
   if ( p._iEdge > _params.size()-1 )
     return false;
+  if ( p._iEdge == _params.size()-1 )
+    return u = 1.;
 
   u = ( _params[ p._iEdge   ] * ( 1 - p._edgeParam ) +
         _params[ p._iEdge+1 ] * p._edgeParam );
@@ -1366,8 +1675,8 @@ bool SMESH_MAT2d::Branch::addDivPntForConcaVertex( std::vector< std::size_t >&
 
   size_t iSeg1  = getBndSegment( maEdges[ i-1 ] );
   size_t iSeg2  = getBndSegment( maEdges[ i ] );
-  bool isConcaPrev = _boundary->IsConcaveSegment( edgeIDs1.back(), iSeg1 );
+  bool isConcaPrev = _boundary->isConcaveSegment( edgeIDs1.back(), iSeg1 );
-  bool isConcaNext = _boundary->IsConcaveSegment( ie1,             iSeg2 );
+  bool isConcaNext = _boundary->isConcaveSegment( ie1,             iSeg2 );
   if ( !isConcaNext && !isConcaPrev )
     return false;
 
@@ -1381,7 +1690,7 @@ bool SMESH_MAT2d::Branch::addDivPntForConcaVertex( std::vector< std::size_t >&
     while ( iNext < maEdges.size() )
     {
       iSeg2 = getBndSegment( maEdges[ iNext ] );
-      if ( _boundary->IsConcaveSegment( ie1, iSeg2 ))
+      if ( _boundary->isConcaveSegment( ie1, iSeg2 ))
         ++iNext;
       else
         break;
@@ -1403,7 +1712,7 @@ bool SMESH_MAT2d::Branch::addDivPntForConcaVertex( std::vector< std::size_t >&
     }
     if ( vertexFound )
     {
-      i = --iNext;
+      iPrev = i = --iNext; // not to add a BP in the moddle
       isConcaveV = true;
     }
   }
@@ -1413,7 +1722,7 @@ bool SMESH_MAT2d::Branch::addDivPntForConcaVertex( std::vector< std::size_t >&
     while ( iPrev-1 >= 0 )
     {
       iSeg1 = getBndSegment( maEdges[ iPrev-1 ] );
-      if ( _boundary->IsConcaveSegment( edgeIDs1.back(), iSeg1 ))
+      if ( _boundary->isConcaveSegment( edgeIDs1.back(), iSeg1 ))
         --iPrev;
       else
         break;
@@ -1423,7 +1732,7 @@ bool SMESH_MAT2d::Branch::addDivPntForConcaVertex( std::vector< std::size_t >&
   if ( iPrev < i-1 || iNext > i )
   {
     // no VERTEX on the opposite EDGE, put the Branch Point in the middle
-    double par1 = _params[ iPrev ], par2 = _params[ iNext ];
+    double par1 = _params[ iPrev+1 ], par2 = _params[ iNext ];
     double midPar = 0.5 * ( par1 + par2 );
     divisionPnt._iEdge = iPrev;
     while ( _params[ divisionPnt._iEdge + 1 ] < midPar )

src/SMESHUtils/SMESH_MAT2d.hxx

@@ -88,6 +88,7 @@ namespace SMESH_MAT2d
   /*!
    * \brief Branch is a set of MA edges enclosed between branch points and/or MA ends.
    *        It's main feature is to return two BoundaryPoint's per a point on it.
+   *        Points on a Branch are defined by [0,1] parameter 
    */
   class SMESHUtils_EXPORT Branch
   {
@@ -114,15 +115,20 @@ namespace SMESH_MAT2d
                                std::vector< std::size_t >& edgeIDs2,
                                std::vector< BranchPoint >& divPoints) const;
 
-    // construction
+    bool isRemoved() const { return _proxyPoint._branch; }
+
+  public: // internal: construction
+
     void init( std::vector<const TVDEdge*>&                maEdges,
                const Boundary*                             boundary,
                std::map< const TVDVertex*, BranchEndType > endType);
     void setBranchesToEnds( const std::vector< Branch >&   branches);
+    BranchPoint getPoint( const TVDVertex* vertex ) const;
+    void setRemoved( const BranchPoint& proxyPoint );
 
-    static void setGeomEdge( std::size_t geomIndex, const TVDEdge* maEdge );
+    static void        setGeomEdge  ( std::size_t geomIndex, const TVDEdge* maEdge );
-    static std::size_t getGeomEdge( const TVDEdge* maEdge );
+    static std::size_t getGeomEdge  ( const TVDEdge* maEdge );
-    static void setBndSegment( std::size_t segIndex, const TVDEdge* maEdge );
+    static void        setBndSegment( std::size_t segIndex, const TVDEdge* maEdge );
     static std::size_t getBndSegment( const TVDEdge* maEdge );
 
   private:
@@ -142,6 +148,7 @@ namespace SMESH_MAT2d
     const Boundary*      _boundary; // face boundary
     BranchEnd            _endPoint1;
     BranchEnd            _endPoint2;
+    BranchPoint          _proxyPoint;
   };
 
   //-------------------------------------------------------------------------------------
@@ -173,7 +180,9 @@ namespace SMESH_MAT2d
 
     bool getBranchPoint( const std::size_t iEdge, double u, BranchPoint& p ) const;
 
-    bool IsConcaveSegment( std::size_t iEdge, std::size_t iSeg ) const;
+    bool isConcaveSegment( std::size_t iEdge, std::size_t iSeg ) const;
+
+    bool moveToClosestEdgeEnd( BoundaryPoint& bp ) const;
 
   private:
     std::vector< BndPoints > _pointsPerEdge;
@@ -201,11 +210,12 @@ namespace SMESH_MAT2d
                const std::vector< TopoDS_Edge >& edges,
                const double                      minSegLen,
                const bool                        ignoreCorners = false );
-    const Boundary&                        getBoundary() const { return _boundary; }
+    std::size_t                            nbBranches() const { return _nbBranches; }
-    const std::vector< Branch >&           getBranches() const { return _branch; }
+    const Branch*                          getBranch(size_t i) const;
     const std::vector< const BranchEnd* >& getBranchPoints() const { return _branchPnt; }
+    const Boundary&                        getBoundary() const { return _boundary; }
 
-    void getPoints( const Branch& branch, std::vector< gp_XY >& points) const;
+    void getPoints( const Branch* branch, std::vector< gp_XY >& points) const;
     Adaptor3d_Curve* make3DCurve(const Branch& branch) const;
 
   private:
@@ -214,6 +224,7 @@ namespace SMESH_MAT2d
     TopoDS_Face                     _face;
     TVD                             _vd;
     std::vector< Branch >           _branch;
+    std::size_t                     _nbBranches; // removed branches ignored
     std::vector< const BranchEnd* > _branchPnt;
     Boundary                        _boundary;
     double                          _scale[2];

src/StdMeshers/StdMeshers_QuadFromMedialAxis_1D2D.cxx

@@ -542,11 +542,11 @@ namespace
   TopoDS_Edge makeEdgeFromMA( SMESH_MesherHelper&            theHelper,
                               const SMESH_MAT2d::MedialAxis& theMA )
   {
-    if ( theMA.getBranches().size() != 1 )
+    if ( theMA.nbBranches() != 1 )
       return TopoDS_Edge();
 
     vector< gp_XY > uv;
-    theMA.getPoints( theMA.getBranches()[0], uv );
+    theMA.getPoints( theMA.getBranch(0), uv );
     if ( uv.size() < 2 )
       return TopoDS_Edge();
 
@@ -775,7 +775,7 @@ namespace
 
     double uMA;
     SMESH_MAT2d::BoundaryPoint bp[2];
-    const SMESH_MAT2d::Branch& branch = theMA.getBranches()[0];
+    const SMESH_MAT2d::Branch& branch = *theMA.getBranch(0);
 
     for ( size_t i = 0; i < theDivPoints.size(); ++i )
     {
@@ -830,7 +830,7 @@ namespace
                          const vector<TopoDS_Edge>& theSinuEdges,
                          const size_t               theSinuSide0Size)
   {
-    if ( theMA.getBranches().size() != 1 )
+    if ( theMA.nbBranches() != 1 )
       return false;
 
     // normalize theMAParams
@@ -858,7 +858,7 @@ namespace
       //   hasComputed = true;
     }
 
-    const SMESH_MAT2d::Branch& branch = theMA.getBranches()[0];
+    const SMESH_MAT2d::Branch& branch = *theMA.getBranch(0);
     SMESH_MAT2d::BoundaryPoint bp[2];
 
     vector< std::size_t > edgeIDs1, edgeIDs2;
@@ -1248,4 +1248,3 @@ bool StdMeshers_QuadFromMedialAxis_1D2D::Evaluate(SMESH_Mesh &         theMesh,
 {
   return StdMeshers_Quadrangle_2D::Evaluate(theMesh,theShape,theResMap);
 }
-