src/VTKViewer/VTKViewer_ConvexTool.cxx

   1 //  This library is distributed in the hope that it will be useful,
   2 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
   3 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
   4 //  Lesser General Public License for more details.
   5 //
   6 //  You should have received a copy of the GNU Lesser General Public
   7 //  License along with this library; if not, write to the Free Software
   8 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
   9 //
  10 //  See http://www.opencascade.org/SALOME/ or email : webmaster.salome@opencascade.org
  11
  12 #include "VTKViewer_ConvexTool.h"
  13
  14 #include <vtkUnstructuredGrid.h>
  15 #include <vtkTriangle.h>
  16 #include <vtkConvexPointSet.h>
  17 #include <vtkMath.h>
  18
  19 #include <set>
  20 #include <iterator>
  21 #include <algorithm>
  22 #include <math.h>
  23
  24 static float FACE_ANGLE_TOLERANCE=1.5;
  25
  26 typedef std::set<vtkIdType> TUIDS; // unique ids
  27 typedef std::map<vtkIdType,TUIDS> TPTOIDS; // id points -> unique ids
  28
  29 namespace CONVEX_TOOL
  30 {
  31 #ifdef _DEBUG_
  32   static int MYDEBUG = 0;
  33 #else
  34   static int MYDEBUG = 0;
  35 #endif
  36
  37 static void GetCenter(vtkUnstructuredGrid* theGrid,TCell theptIds,float *center)
  38 {
  39   float *p;
  40   center[0] = center[1] = center[2] = 0.0;
  41
  42   int numIds=theptIds.size();
  43
  44   // get the center of the cell
  45   for (int i=0; i < numIds; i++)
  46     {
  47       p = theGrid->GetPoint(theptIds[i]);
  48       for (int j=0; j < 3; j++)
  49         {
  50           center[j] += p[j];
  51         }
  52     }
  53   for (int j=0; j<3; j++)
  54     {
  55       center[j] /= numIds;
  56     }
  57 }
  58
  59 static void ReverseIds(TCell &theIds)
  60 {
  61   int i;
  62   vtkIdType tmp;
  63   vtkIdType npts=theIds.size();
  64
  65   for(i=0;i<(npts/2);i++){
  66     tmp = theIds[i];
  67     theIds[i] = theIds[npts-i-1];
  68     theIds[npts-i-1] = tmp;
  69   }
  70 }
  71
  72 // caclulation of connected faces (faceId -> (faceId1,faceId2, ...))
  73 void GetFriends(const TPTOIDS p2faces,const TCellArray f2points,TPTOIDS& face2face_output)
  74 {
  75   TCellArray::const_iterator f2pIter = f2points.begin();
  76
  77   for( ; f2pIter!=f2points.end() ; f2pIter++ ){
  78     vtkIdType faceId = f2pIter->first;
  79     TCell face_points = f2pIter->second;
  80     int nb_face_points = face_points.size();
  81
  82     vtkIdType id1;
  83     vtkIdType id2;
  84     TPTOIDS::const_iterator faces1;
  85     TPTOIDS::const_iterator faces2;
  86
  87     id1 = face_points[0];
  88     faces1 = p2faces.find(id1);
  89
  90     TUIDS output_faces;
  91
  92     for(int i=1 ; i<nb_face_points ; i++ ){
  93
  94       id2 = face_points[i];
  95
  96       faces2 = p2faces.find(id2);
  97
  98       std::set_intersection(faces1->second.begin(), faces1->second.end(), faces2->second.begin(), faces2->second.end(),
  99         std::inserter(output_faces,output_faces.begin()));
 100
 101       id1 = id2;
 102       faces1 = faces2;
 103     }
 104     id1 = face_points[0];
 105     faces1 = p2faces.find(id1);
 106     std::set_intersection(faces1->second.begin(), faces1->second.end(), faces2->second.begin(), faces2->second.end(),
 107       std::inserter(output_faces,output_faces.begin()));
 108
 109     output_faces.erase(faceId); // erase the face id for which we found friends
 110
 111     if(MYDEBUG){
 112       cout << "fId[" << faceId <<"]: ";
 113       std::copy(output_faces.begin(), output_faces.end(), std::ostream_iterator<vtkIdType>(cout, " "));
 114       cout << endl;
 115     }
 116
 117     face2face_output[faceId] = output_faces;
 118   }
 119 }
 120
 121 bool IsConnectedFacesOnOnePlane( vtkUnstructuredGrid* theGrid,
 122                                  vtkIdType theFId1, vtkIdType theFId2,
 123                                                          TUIDS FpIds1, TUIDS FpIds2 )
 124 {
 125   bool status = false;
 126   TUIDS common_ids;
 127   std::set_intersection(FpIds1.begin(), FpIds1.end(), FpIds2.begin(), FpIds2.end(),
 128                         std::inserter(common_ids,common_ids.begin()));
 129
 130   /*           Number of common ids = 2 (A1,A2)
 131
 132
 133                 _ _ _ _ _      _ _ _ _    vectors:
 134                |         \   /         |   v1 {A2,A1}
 135                           \ /              v2 {A1,B1}
 136                |           | A2        |   v3 {A1,C1}
 137                            |
 138                |           |           |
 139                            |
 140                |           | A1        |
 141                           / \
 142                |_ _ _ _ _/   \_ _ _ _ _|
 143                B2        B1   C1        C2
 144
 145    */
 146   TUIDS::iterator common_iter = common_ids.begin();
 147   if(common_ids.size() == 2){
 148     TUIDS::iterator loc_id1_0 = FpIds1.find(*(common_iter));
 149     common_iter++;
 150     TUIDS::iterator loc_id1_1 = FpIds1.find(*(common_iter));
 151
 152     TUIDS::iterator loc_id2_0 = FpIds1.begin();
 153     TUIDS::iterator loc_id2_1 = FpIds2.begin();
 154
 155     vtkIdType A1 = *loc_id1_0;
 156     vtkIdType A2 = *loc_id1_1;
 157     vtkIdType B1;
 158     vtkIdType C1;
 159
 160     for(;loc_id2_0!=FpIds1.end();loc_id2_0++)
 161       if(*loc_id2_0 != A1 && *loc_id2_0!= A2){
 162         B1 = *loc_id2_0;
 163         break;
 164       }
 165     for(;loc_id2_1!=FpIds2.end();loc_id2_1++)
 166       if(*loc_id2_1 != A1 && *loc_id2_1!= A2){
 167         C1 = *loc_id2_1;
 168         break;
 169       }
 170     if(MYDEBUG) cout <<endl;
 171     if(MYDEBUG) cout << "FId_1="<<theFId1<<" FId_2="<<theFId2<<endl;
 172     if(MYDEBUG) cout << "   A1="<<A1<<" A2="<<A2<<" B1="<<B1<<" C1="<<C1<<" ->";
 173     float *p[4];
 174     float v1[3],v2[3],v3[3];
 175     p[0] = theGrid->GetPoint(A1);
 176     p[1] = theGrid->GetPoint(A2);
 177     p[2] = theGrid->GetPoint(B1);
 178     p[3] = theGrid->GetPoint(C1);
 179
 180     for(int i=0;i<3;i++){
 181       v1[i] = p[1][i] - p[0][i];
 182       v2[i] = p[2][i] - p[0][i];
 183       v3[i] = p[3][i] - p[0][i];
 184     }
 185
 186
 187     float vec_b1[3];
 188     vtkMath::Cross(v2,v1,vec_b1);
 189     float vec_b2[3];
 190     vtkMath::Cross(v1,v3,vec_b2);
 191
 192     float b1 = vtkMath::Norm(vec_b1);
 193
 194     float b2 = vtkMath::Norm(vec_b2);
 195
 196     float angle=180*acosf(vtkMath::Dot(vec_b1,vec_b2)/(b1*b2))/vtkMath::Pi();
 197
 198     if( angle <= FACE_ANGLE_TOLERANCE )
 199       status = true;
 200     if (MYDEBUG){
 201       for(int k=0;k<4;k++){
 202         cout << " (";
 203         for(int j=0;j<2;j++){
 204           cout << p[k][j] << ",";
 205         }
 206         cout << p[k][2] << ")   ";
 207       }
 208       cout << "angle="<<angle<<endl;
 209     }
 210
 211   } else if (common_ids.size() >2){
 212     // not implemented yet
 213     if(MYDEBUG) cout << "Warning! VTKViewer_ConvexTool::IsConnectedFacesOnOnePlane()";
 214   } else {
 215     // one or no connection ... continue
 216   }
 217
 218   return status;
 219 }
 220
 221 void GetAllFacesOnOnePlane( TPTOIDS theFaces, vtkIdType faceId,
 222                             TUIDS &new_faces, TCell &new_faces_v2 )
 223 {
 224   if (new_faces.find(faceId) != new_faces.end()) return;
 225
 226   new_faces.insert(new_faces.begin(),faceId);
 227   new_faces_v2.push_back(faceId);
 228
 229   TPTOIDS::const_iterator aIter1 = theFaces.find(faceId);
 230   if(aIter1!=theFaces.end()){
 231     TUIDS::const_iterator aIter2 = (aIter1->second).begin();
 232     for(;aIter2!=(aIter1->second).end();aIter2++){
 233       if (new_faces.find(*aIter2) != new_faces.end()) continue;
 234       GetAllFacesOnOnePlane(theFaces,*aIter2,
 235                             new_faces,new_faces_v2); // recurvise
 236     }
 237   }
 238   return;
 239 }
 240
 241 void GetSumm(TCell v1,TCell v2,TCell &output)
 242 {
 243   output.clear();
 244
 245   if(MYDEBUG) cout << "========================================="<<endl;
 246   if(MYDEBUG) cout << "v1:";
 247   if(MYDEBUG) std::copy(v1.begin(), v1.end(), std::ostream_iterator<vtkIdType>(cout, " "));
 248   if(MYDEBUG) cout << "\tv2:";
 249   if(MYDEBUG) std::copy(v2.begin(), v2.end(), std::ostream_iterator<vtkIdType>(cout, " "));
 250   if(MYDEBUG) cout << endl;
 251
 252   TUIDS v1_set;
 253   std::copy(v1.begin(), v1.end(), std::inserter(v1_set,v1_set.begin()));
 254   TUIDS v2_set;
 255   std::copy(v2.begin(), v2.end(), std::inserter(v2_set,v2_set.begin()));
 256   TUIDS tmpIntersection;
 257   std::set_intersection(v1_set.begin(),v1_set.end(),v2_set.begin(),v2_set.end(), std::inserter(tmpIntersection,tmpIntersection.begin()));
 258   if(MYDEBUG) std::copy(tmpIntersection.begin(),tmpIntersection.end(), std::ostream_iterator<vtkIdType>(cout, " "));
 259   if(MYDEBUG) cout << endl;
 260
 261   if(tmpIntersection.size() < 2)
 262     if(MYDEBUG) cout << __FILE__ << "[" << __LINE__ << "]: Warning ! Wrong ids" << endl;
 263
 264   TCell::iterator v1_iter = v1.begin();
 265
 266   for(;v1_iter!=v1.end();v1_iter++){
 267
 268     vtkIdType curr_id = *v1_iter;
 269
 270     output.push_back(curr_id);
 271
 272     if(tmpIntersection.find(curr_id) != tmpIntersection.end()){
 273       TCell::iterator v1_iter_tmp;
 274       v1_iter_tmp = v1_iter;
 275       v1_iter++;
 276
 277       if(v1_iter==v1.end()) v1_iter=v1.begin();
 278
 279       curr_id = *v1_iter;
 280
 281       if(tmpIntersection.find(curr_id) != tmpIntersection.end()){
 282         TCell::iterator v2_iter = v2.begin();
 283         for(;v2_iter!=v2.end();v2_iter++){
 284           vtkIdType v2_id = *v2_iter;
 285           if(tmpIntersection.find(v2_id) == tmpIntersection.end())
 286             output.push_back(v2_id);
 287         }
 288       }
 289
 290       v1_iter = v1_iter_tmp;
 291       curr_id = *v1_iter;
 292
 293     }
 294   }
 295
 296   if(MYDEBUG) cout << "Result: " ;
 297   if(MYDEBUG) std::copy(output.begin(),output.end(),std::ostream_iterator<vtkIdType>(cout, " "));
 298   if(MYDEBUG) cout << endl;
 299 }
 300
 301 void GetPolygonalFaces(vtkUnstructuredGrid* theGrid,int cellId,TCellArray &outputCellArray)
 302 {
 303   if (theGrid->GetCellType(cellId) == VTK_CONVEX_POINT_SET){
 304     // get vtkCell type = VTK_CONVEX_POINT_SET
 305     if(vtkConvexPointSet* convex = static_cast<vtkConvexPointSet*>(theGrid->GetCell(cellId))){
 306       TCellArray f2points;
 307       float convex_center[3]; // convex center point coorinat
 308       int aNbFaces = convex->GetNumberOfFaces();
 309       int numPts = convex->GetNumberOfPoints();
 310       TCell convex_ids;
 311       TPTOIDS p2faces; // key=pointId , value facesIds set
 312
 313       for(int i=0;i<numPts;i++)
 314           convex_ids.push_back(convex->GetPointId(i));
 315
 316       GetCenter(theGrid,convex_ids,convex_center);
 317
 318       for (vtkIdType faceId=0; faceId < aNbFaces; faceId++){
 319         vtkCell *aFace = convex->GetFace(faceId);
 320         int numFacePts = aFace->GetNumberOfPoints();
 321         TCell aIds;
 322
 323   int i = 0;
 324         for(i=0;i<numFacePts;i++)
 325           aIds.push_back(aFace->GetPointId(i));
 326
 327         float v_a[3],v_b[3],v_convex2face[3]; // vectors
 328         float *id_0,*id_1,*id_n;
 329         /*=============================================
 330                         ,+- - - -  _
 331                    _   / id_n   |  v_b {id_0,id_n}
 332                   v_b /            _
 333                      /          |  v_a {id_0,id_1}
 334                     /
 335                    /            |
 336                   + id_0
 337                    \
 338                   _ \           |
 339                  v_a \
 340                       \ id_1    |
 341                        "+- - - -
 342
 343          ============================================*/
 344         id_0 = theGrid->GetPoint(aIds[0]);
 345         id_1 = theGrid->GetPoint(aIds[1]);
 346         id_n = theGrid->GetPoint(aIds[numFacePts-1]);
 347
 348         for(i=0;i<3;i++){
 349           v_a[i] = id_1[i] - id_0[i];
 350           v_b[i] = id_n[i] - id_0[i];
 351           v_convex2face[i] = id_0[i] - convex_center[i];
 352         }
 353
 354         if (vtkMath::Determinant3x3(v_a,v_b,v_convex2face) < 0){
 355           ReverseIds(aIds);
 356         }
 357
 358         for(i=0;i<(int)aIds.size();i++){
 359           TUIDS &acell = p2faces[aIds[i]];
 360           acell.insert(faceId);
 361         }
 362
 363         f2points[faceId] = aIds;
 364
 365       }
 366
 367       TPTOIDS face2face;
 368       GetFriends(p2faces,f2points,face2face);
 369
 370       TPTOIDS face2points;
 371
 372       // copy TCellArray::f2points to TPTOIDS::face2points
 373       for(TCellArray::iterator f2points_iter=f2points.begin();
 374           f2points_iter!=f2points.end();
 375           f2points_iter++){
 376
 377         TUIDS tmp;
 378         for(TCell::iterator points_iter=(f2points_iter->second).begin();
 379             points_iter!=(f2points_iter->second).end();
 380             points_iter++)
 381           tmp.insert(*points_iter);
 382
 383         face2points[f2points_iter->first] = tmp;
 384       } // end copy
 385
 386
 387       TPTOIDS new_face2faces; // which connected and in one plane
 388
 389       TPTOIDS::const_iterator aF2FIter = face2face.begin();
 390       for(;aF2FIter!=face2face.end();aF2FIter++){
 391         vtkIdType f_key = aF2FIter->first;
 392         TUIDS &faces = new_face2faces[f_key];
 393         //faces.insert(f_key);
 394         TUIDS f_friends = aF2FIter->second;
 395         TUIDS::const_iterator a_friends_iter = f_friends.begin();
 396         for(;a_friends_iter!=f_friends.end();a_friends_iter++){
 397           vtkIdType friend_id = *a_friends_iter;
 398           if( IsConnectedFacesOnOnePlane(theGrid,f_key,friend_id,
 399                                         (face2points.find(f_key))->second,
 400                                         (face2points.find(friend_id))->second)){
 401             faces.insert(friend_id);
 402           } // end if
 403
 404         } // end a_friends_iter
 405       } // end aF2FIter
 406
 407       if(MYDEBUG)
 408       {
 409         TPTOIDS::const_iterator new_face2face_iter = new_face2faces.begin();
 410         cout << "Connected faces and on plane:" << endl;
 411         for(;new_face2face_iter!=new_face2faces.end();new_face2face_iter++){
 412           cout << "Group ["<<new_face2face_iter->first<<"] :";
 413           TUIDS::const_iterator new_faces_iter = (new_face2face_iter->second).begin();
 414           for(;new_faces_iter!=(new_face2face_iter->second).end();new_faces_iter++)
 415             cout << " " << *new_faces_iter ;
 416           cout << endl;
 417         }
 418       }
 419
 420       TPTOIDS output_newid2face;
 421       TCellArray output_newid2face_v2;
 422       {
 423         TUIDS already_in;
 424         TUIDS already_in_tmp;
 425         int k=0;
 426         TPTOIDS::const_iterator new_face2face_iter = new_face2faces.begin();
 427         for(;new_face2face_iter!=new_face2faces.end();new_face2face_iter++){
 428           if(already_in.find(new_face2face_iter->first) != already_in.end())
 429             continue;
 430           if(new_face2face_iter->second.size() > 1)
 431             continue;
 432
 433           TCell &tmp_v2 = output_newid2face_v2[k];
 434           tmp_v2.push_back(new_face2face_iter->first);
 435           already_in.insert(new_face2face_iter->first);
 436
 437           TUIDS::const_iterator new_faces_iter = (new_face2face_iter->second).begin();
 438           for(;new_faces_iter!=(new_face2face_iter->second).end();new_faces_iter++){
 439             if(already_in.find(*new_faces_iter) != already_in.end()) continue;
 440             already_in.insert(*new_faces_iter);
 441
 442             already_in_tmp.clear();
 443             already_in_tmp.insert(new_face2face_iter->first);
 444
 445             TUIDS &tmp = output_newid2face[k];
 446             GetAllFacesOnOnePlane(new_face2faces,*new_faces_iter,
 447                                   already_in_tmp,tmp_v2);
 448
 449             for(TUIDS::const_iterator aIter=already_in_tmp.begin();
 450                 aIter!=already_in_tmp.end();
 451                 aIter++)
 452               {
 453                 already_in.insert(*aIter);
 454                 tmp.insert(*aIter);
 455               }
 456           }
 457           k++;
 458         }
 459       }
 460
 461       if(MYDEBUG) {
 462         cout << "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
 463         {
 464           TPTOIDS::const_iterator new_face2face_iter = output_newid2face.begin();
 465           for(;new_face2face_iter!=output_newid2face.end();new_face2face_iter++){
 466             cout << "Group ["<<new_face2face_iter->first<<"] :";
 467             TUIDS::const_iterator new_faces_iter = (new_face2face_iter->second).begin();
 468             for(;new_faces_iter!=(new_face2face_iter->second).end();new_faces_iter++)
 469               cout << " " << *new_faces_iter ;
 470             cout << endl;
 471           }
 472         }
 473         cout << "++++++++++++++++++++++++++ +++++++++++++++++++++++++++++++"<<endl;
 474         cout << "++++++++++++++++++++++++++ +++++++++++++++++++++++++++++++"<<endl;
 475         cout << "+++++++++++++++++++++++ ++ ++ ++++++++++++++++++++++++++++"<<endl;
 476         cout << "+++++++++++++++++++++++++   ++++++++++++++++++++++++++++++"<<endl;
 477         cout << "++++++++++++++++++++++++++ +++++++++++++++++++++++++++++++"<<endl;
 478         {
 479           TCellArray::const_iterator new_face2face_iter = output_newid2face_v2.begin();
 480           for(;new_face2face_iter!=output_newid2face_v2.end();new_face2face_iter++){
 481             cout << "Group ["<<new_face2face_iter->first<<"] :";
 482             TCell::const_iterator new_faces_iter = (new_face2face_iter->second).begin();
 483             for(;new_faces_iter!=(new_face2face_iter->second).end();new_faces_iter++)
 484               cout << " " << *new_faces_iter ;
 485             cout << endl;
 486           }
 487         }
 488       }
 489       TCellArray output_new_face2ids;
 490 //       {
 491 //      TPTOIDS::const_iterator new_face2face_iter = output_newid2face.begin();
 492 //      for(;new_face2face_iter!=output_newid2face.end();new_face2face_iter++){
 493
 494 //        vtkIdType new_faceId = new_face2face_iter->first;
 495 //        TUIDS::const_iterator new_faces_iter = (new_face2face_iter->second).begin();
 496 //        vtkIdType fId0 = *new_faces_iter;
 497 //        TCellArray::const_iterator pIds0_iter = f2points.find(fId0);
 498 //        TCell pIds0 = pIds0_iter->second;
 499 //        TCell &output = output_new_face2ids[new_faceId];
 500 //        new_faces_iter++;
 501 //        if(new_face2face_iter->second.size() > 2 ){}
 502 //        for(;new_faces_iter!=(new_face2face_iter->second).end();new_faces_iter++){
 503
 504 //          vtkIdType faceId = *new_faces_iter;
 505 //          // find how much nodes in face (f2points)
 506 //          TCellArray::const_iterator pIds_iter = f2points.find(faceId);
 507 //          TCell pIds = pIds_iter->second;
 508
 509 //          GetSumm(pIds0,pIds,output);
 510 //          pIds0 = output;
 511
 512 //        } // end new_faces_iter
 513
 514 //      } // new_face2face_iter
 515 //       }
 516
 517       {
 518         TCellArray::const_iterator new_face2face_iter = output_newid2face_v2.begin();
 519         for(;new_face2face_iter!=output_newid2face_v2.end();new_face2face_iter++){
 520
 521           vtkIdType new_faceId = new_face2face_iter->first;
 522           TCell::const_iterator new_faces_iter = (new_face2face_iter->second).begin();
 523           vtkIdType fId0 = *new_faces_iter;
 524           TCellArray::const_iterator pIds0_iter = f2points.find(fId0);
 525           TCell pIds0 = pIds0_iter->second;
 526           TCell &output = output_new_face2ids[new_faceId];
 527           new_faces_iter++;
 528           if(new_face2face_iter->second.size() == 1 ){
 529             TCellArray::const_iterator pIds_iter = f2points.find(fId0);
 530             TCell pIds = pIds_iter->second;
 531             output = pIds;
 532             continue;
 533           }
 534           for(;new_faces_iter!=(new_face2face_iter->second).end();new_faces_iter++){
 535
 536             vtkIdType faceId = *new_faces_iter;
 537             // find how much nodes in face (f2points)
 538             TCellArray::const_iterator pIds_iter = f2points.find(faceId);
 539             TCell pIds = pIds_iter->second;
 540
 541             GetSumm(pIds0,pIds,output);
 542             pIds0 = output;
 543
 544           } // end new_faces_iter
 545
 546         } // new_face2face_iter
 547       }
 548
 549       if(MYDEBUG) {
 550         cout << "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
 551         cout << "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
 552         cout << "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
 553       }
 554       outputCellArray = output_new_face2ids;//f2points;
 555     }
 556   } else {
 557     // not implemented
 558   }
 559 }
 560 }