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