src/VTKViewer/VTKViewer_ConvexTool.cxx

   1 // Copyright (C) 2005  OPEN CASCADE, CEA/DEN, EDF R&D, PRINCIPIA R&D
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License.
   7 //
   8 // This library is distributed in the hope that it will be useful
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/
  18 //
  19 //  This library is distributed in the hope that it will be useful,
  20 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
  21 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  22 //  Lesser General Public License for more details.
  23 //
  24 //  You should have received a copy of the GNU Lesser General Public
  25 //  License along with this library; if not, write to the Free Software
  26 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  27 //
  28 //  See http://www.opencascade.org/SALOME/ or email : webmaster.salome@opencascade.org
  29
  30 #include "VTKViewer_ConvexTool.h"
  31
  32 #include <vtkUnstructuredGrid.h>
  33 #include <vtkTriangle.h>
  34 #include <vtkConvexPointSet.h>
  35 #include <vtkUnstructuredGridWriter.h>
  36 #include <vtkMath.h>
  37 #include <vtkSmartPointer.h>
  38
  39 #include <set>
  40 #include <iterator>
  41 #include <algorithm>
  42 #include <math.h>
  43
  44 typedef vtkUnstructuredGrid TInput;
  45
  46 typedef std::set<vtkIdType> TUIDS; // unique ids
  47 typedef std::map<vtkIdType,TUIDS> TPTOIDS; // id points -> unique ids
  48
  49 namespace CONVEX_TOOL
  50 {
  51   // all pairs
  52   typedef std::pair<vtkIdType,vtkIdType> TPair;
  53   typedef std::set<TPair> TSet;
  54
  55   void
  56   WriteToFile(vtkUnstructuredGrid* theDataSet, const std::string& theFileName)
  57   {
  58     vtkUnstructuredGridWriter* aWriter = vtkUnstructuredGridWriter::New();
  59     aWriter->SetFileName(theFileName.c_str());
  60     aWriter->SetInput(theDataSet);
  61     aWriter->Write();
  62     aWriter->Delete();
  63   }
  64
  65 static float FACE_ANGLE_TOLERANCE=1.5;
  66 #define EPS 1.0e-38
  67 #define EPS_T 1.0e-3
  68
  69 #ifdef _DEBUG_
  70   static int MYDEBUG = 0;
  71   static int MYDEBUG_REMOVE = 0;
  72 #else
  73   static int MYDEBUG = 0;
  74   static int MYDEBUG_REMOVE = 0;
  75 #endif
  76
  77 /*! \fn static void GetCenter(TInput* theGrid,TCell theptIds,float *center)
  78  * \brief Calculation of geometry center.
  79  * \param theGrid - TInput cell.
  80  * \param theptIds - point ids.
  81  * \retval center - output array[3] with coordinates of geometry center.
  82  */
  83 static void GetCenter(vtkPoints* thePoints,float center[3])
  84 {
  85   float p[3];
  86   center[0] = center[1] = center[2] = 0.0;
  87
  88   int numPts = thePoints->GetNumberOfPoints();
  89   if (numPts == 0) return;
  90
  91   // get the center of the cell
  92   for (int i = 0; i < numPts; i++)
  93   {
  94     thePoints->GetPoint(i, p);
  95     for (int j = 0; j < 3; j++)
  96     {
  97       center[j] += p[j];
  98     }
  99   }
 100   for (int j = 0; j < 3; j++)
 101   {
 102     center[j] /= numPts;
 103   }
 104 }
 105
 106 /*! \fn static void ReverseIds(TCell &theIds)
 107  * \brief Reverse ids.
 108  * \param theIds - points ids.
 109  * \retval theIds - example input:(1,2,3,4) -> output:(4,3,2,1)
 110  */
 111 static void ReverseIds(TCell &theIds)
 112 {
 113   int i;
 114   vtkIdType tmp;
 115   vtkIdType npts=theIds.size();
 116
 117   for(i=0;i<(npts/2);i++){
 118     tmp = theIds[i];
 119     theIds[i] = theIds[npts-i-1];
 120     theIds[npts-i-1] = tmp;
 121   }
 122 }
 123
 124 /*! \fn void GetFriends(const TPTOIDS p2faces,const TCellArray f2points,TPTOIDS& face2face_output)
 125  * \brief Caclulation of connected faces (faceId -> (faceId1,faceId2, ...))
 126  * \param p2faces - point to faces ids map.
 127  * \param f2points - faces to points ids map.
 128  * \retval face2face_output - faces to faces ids map.
 129  */
 130 void GetFriends(const TPTOIDS p2faces,const TCellArray f2points,TPTOIDS& face2face_output)
 131 {
 132   TCellArray::const_iterator f2pIter = f2points.begin();
 133
 134   for( ; f2pIter!=f2points.end() ; f2pIter++ ){
 135     vtkIdType faceId = f2pIter->first;
 136     TCell face_points = f2pIter->second;
 137     int nb_face_points = face_points.size();
 138
 139     vtkIdType id1;
 140     vtkIdType id2;
 141     TPTOIDS::const_iterator faces1;
 142     TPTOIDS::const_iterator faces2;
 143
 144     id1 = face_points[0];
 145     faces1 = p2faces.find(id1);
 146
 147     TUIDS output_faces;
 148
 149     for(int i=1 ; i<nb_face_points ; i++ ){
 150
 151       id2 = face_points[i];
 152
 153       faces2 = p2faces.find(id2);
 154
 155       std::set_intersection(faces1->second.begin(), faces1->second.end(), faces2->second.begin(), faces2->second.end(),
 156         std::inserter(output_faces,output_faces.begin()));
 157
 158       id1 = id2;
 159       faces1 = faces2;
 160     }
 161     id1 = face_points[0];
 162     faces1 = p2faces.find(id1);
 163     std::set_intersection(faces1->second.begin(), faces1->second.end(), faces2->second.begin(), faces2->second.end(),
 164       std::inserter(output_faces,output_faces.begin()));
 165
 166     output_faces.erase(faceId); // erase the face id for which we found friends
 167
 168     if(MYDEBUG){
 169       cout << "fId[" << faceId <<"]: ";
 170       std::copy(output_faces.begin(), output_faces.end(), std::ostream_iterator<vtkIdType>(cout, " "));
 171       cout << endl;
 172     }
 173
 174     face2face_output[faceId] = output_faces;
 175   }
 176 }
 177
 178 /*! \fn bool IsConnectedFacesOnOnePlane( TInput* theGrid,vtkIdType theFId1, vtkIdType theFId2,TUIDS FpIds1, TUIDS FpIds2 )
 179  * \brief Check is connected faces on one plane.
 180  * \param theGrid - TInput
 181  * \param theFId1 - id of first face
 182  * \param theFId2 - id of second face
 183  * \param FpIds1  - first face points ids.
 184  * \param FpIds2  - second face points ids.
 185  * \return TRUE if two faces on one plane, else FALSE.
 186  */
 187 bool IsConnectedFacesOnOnePlane( TInput* theGrid,
 188                                  vtkIdType theFId1, vtkIdType theFId2,
 189                                  TUIDS FpIds1, TUIDS FpIds2 )
 190 {
 191   bool status = false;
 192   TUIDS common_ids;
 193   std::set_intersection(FpIds1.begin(), FpIds1.end(), FpIds2.begin(), FpIds2.end(),
 194                         std::inserter(common_ids,common_ids.begin()));
 195
 196   /*           Number of common ids = 2 (A1,A2)
 197
 198
 199                 _ _ _ _ _      _ _ _ _    vectors:
 200                |         \   /         |   v1 {A2,A1}
 201                           \ /              v2 {A1,B1}
 202                |           | A2        |   v3 {A1,C1}
 203                            |
 204                |           |           |
 205                            |
 206                |           | A1        |
 207                           / \
 208                |_ _ _ _ _/   \_ _ _ _ _|
 209                B2        B1   C1        C2
 210
 211    */
 212   TUIDS::iterator common_iter = common_ids.begin();
 213   if(common_ids.size() == 2){
 214     TUIDS::iterator loc_id1_0 = FpIds1.find(*(common_iter));
 215     common_iter++;
 216     TUIDS::iterator loc_id1_1 = FpIds1.find(*(common_iter));
 217
 218     TUIDS::iterator loc_id2_0 = FpIds1.begin();
 219     TUIDS::iterator loc_id2_1 = FpIds2.begin();
 220
 221     vtkIdType A1 = *loc_id1_0;
 222     vtkIdType A2 = *loc_id1_1;
 223     vtkIdType B1;
 224     vtkIdType C1;
 225
 226     for(;loc_id2_0!=FpIds1.end();loc_id2_0++)
 227       if(*loc_id2_0 != A1 && *loc_id2_0!= A2){
 228         B1 = *loc_id2_0;
 229         break;
 230       }
 231     for(;loc_id2_1!=FpIds2.end();loc_id2_1++)
 232       if(*loc_id2_1 != A1 && *loc_id2_1!= A2){
 233         C1 = *loc_id2_1;
 234         break;
 235       }
 236     if(MYDEBUG) cout <<endl;
 237     if(MYDEBUG) cout << "FId_1="<<theFId1<<" FId_2="<<theFId2<<endl;
 238     if(MYDEBUG) cout << "   A1="<<A1<<" A2="<<A2<<" B1="<<B1<<" C1="<<C1<<" ->";
 239     float *p[4];
 240     float v1[3],v2[3],v3[3];
 241     p[0] = theGrid->GetPoint(A1);
 242     p[1] = theGrid->GetPoint(A2);
 243     p[2] = theGrid->GetPoint(B1);
 244     p[3] = theGrid->GetPoint(C1);
 245
 246     for(int i=0;i<3;i++){
 247       v1[i] = p[1][i] - p[0][i];
 248       v2[i] = p[2][i] - p[0][i];
 249       v3[i] = p[3][i] - p[0][i];
 250     }
 251
 252
 253     float vec_b1[3];
 254     vtkMath::Cross(v2,v1,vec_b1);
 255     float vec_b2[3];
 256     vtkMath::Cross(v1,v3,vec_b2);
 257
 258     float b1 = vtkMath::Norm(vec_b1);
 259
 260     float b2 = vtkMath::Norm(vec_b2);
 261     float aCos = vtkMath::Dot(vec_b1,vec_b2)/(b1*b2);
 262
 263     float angle=90.0;
 264     angle = aCos>=1.0 ? 0.0 : 180*acosf(aCos)/vtkMath::Pi();
 265
 266     if( angle <= FACE_ANGLE_TOLERANCE)
 267       status = true;
 268     if (MYDEBUG){
 269       for(int k=0;k<4;k++){
 270         cout << " (";
 271         for(int j=0;j<2;j++){
 272           cout << p[k][j] << ",";
 273         }
 274         cout << p[k][2] << ")   ";
 275       }
 276       if(status) cout << "angle="<<angle<<" status="<<status<<endl;
 277     }
 278
 279   } else if (common_ids.size() >2){
 280     // not implemented yet
 281     if(MYDEBUG) cout << "Warning! VTKViewer_ConvexTool::IsConnectedFacesOnOnePlane()";
 282   } else {
 283     // one or no connection ... continue
 284   }
 285
 286   return status;
 287 }
 288
 289 /*! \fn void GetAllFacesOnOnePlane( TPTOIDS theFaces, vtkIdType faceId,TUIDS &new_faces, TCell &new_faces_v2 )
 290  * \brief Calculate faces which on one plane.
 291  * \param theFaces -
 292  * \param faceId -
 293  * \param new_faces -
 294  * \param new_faces_v2 -
 295  */
 296 void GetAllFacesOnOnePlane( TPTOIDS theFaces, vtkIdType faceId,
 297                             TUIDS &new_faces, TCell &new_faces_v2 )
 298 {
 299   if (new_faces.find(faceId) != new_faces.end()) return;
 300
 301   new_faces.insert(new_faces.begin(),faceId);
 302   new_faces_v2.push_back(faceId);
 303
 304   TPTOIDS::const_iterator aIter1 = theFaces.find(faceId);
 305   if(aIter1!=theFaces.end()){
 306     TUIDS::const_iterator aIter2 = (aIter1->second).begin();
 307     for(;aIter2!=(aIter1->second).end();aIter2++){
 308       if (new_faces.find(*aIter2) != new_faces.end()) continue;
 309       GetAllFacesOnOnePlane(theFaces,*aIter2,
 310                             new_faces,new_faces_v2); // recurvise
 311     }
 312   }
 313   return;
 314 }
 315
 316 /*! \fn void GetSumm(TCell v1,TCell v2,TCell &output)
 317  * \brief Gluing two faces (gluing points ids)
 318  * \param v1 - first face
 319  * \param v2 - second face
 320  * \param output - output face.
 321  */
 322 void GetSumm(TCell v1,TCell v2,TCell &output)
 323 {
 324   output.clear();
 325
 326   if(MYDEBUG) cout << "========================================="<<endl;
 327   if(MYDEBUG) cout << "v1:";
 328   if(MYDEBUG) std::copy(v1.begin(), v1.end(), std::ostream_iterator<vtkIdType>(cout, " "));
 329   if(MYDEBUG) cout << "\tv2:";
 330   if(MYDEBUG) std::copy(v2.begin(), v2.end(), std::ostream_iterator<vtkIdType>(cout, " "));
 331   if(MYDEBUG) cout << endl;
 332
 333   TUIDS v1_set;
 334   std::copy(v1.begin(), v1.end(), std::inserter(v1_set,v1_set.begin()));
 335   TUIDS v2_set;
 336   std::copy(v2.begin(), v2.end(), std::inserter(v2_set,v2_set.begin()));
 337   TUIDS tmpIntersection;
 338   std::set_intersection(v1_set.begin(),v1_set.end(),v2_set.begin(),v2_set.end(), std::inserter(tmpIntersection,tmpIntersection.begin()));
 339   if(MYDEBUG) std::copy(tmpIntersection.begin(),tmpIntersection.end(), std::ostream_iterator<vtkIdType>(cout, " "));
 340   if(MYDEBUG) cout << endl;
 341
 342   if(tmpIntersection.size() < 2)
 343     if(MYDEBUG) cout << __FILE__ << "[" << __LINE__ << "]: Warning ! Wrong ids" << endl;
 344
 345   TCell::iterator v1_iter = v1.begin();
 346
 347   for(;v1_iter!=v1.end();v1_iter++){
 348
 349     vtkIdType curr_id = *v1_iter;
 350
 351     output.push_back(curr_id);
 352
 353     if(tmpIntersection.find(curr_id) != tmpIntersection.end()){
 354       TCell::iterator v1_iter_tmp;
 355       v1_iter_tmp = v1_iter;
 356       v1_iter++;
 357
 358       if(v1_iter==v1.end()) v1_iter=v1.begin();
 359
 360       curr_id = *v1_iter;
 361
 362       if(tmpIntersection.find(curr_id) != tmpIntersection.end()){
 363         TCell::iterator v2_iter = v2.begin();
 364         for(;v2_iter!=v2.end();v2_iter++){
 365           vtkIdType v2_id = *v2_iter;
 366           if(tmpIntersection.find(v2_id) == tmpIntersection.end())
 367             output.push_back(v2_id);
 368         }
 369       }
 370
 371       v1_iter = v1_iter_tmp;
 372       curr_id = *v1_iter;
 373
 374     }
 375   }
 376
 377   if(MYDEBUG) cout << "Result: " ;
 378   if(MYDEBUG) std::copy(output.begin(),output.end(),std::ostream_iterator<vtkIdType>(cout, " "));
 379   if(MYDEBUG) cout << endl;
 380 }
 381
 382 static void GetAndRemoveIdsOnOneLine(vtkPoints* points,
 383                                      TUIDS input_points_ids,
 384                                      TUIDS input_two_points_ids,
 385                                      TUIDS& out_two_points_ids,
 386                                      TUIDS& removed_points_ids){
 387   if (MYDEBUG_REMOVE) cout << EPS <<endl;
 388   float P[3][3];
 389   vtkIdType current_points_ids[2];
 390   if(MYDEBUG_REMOVE)
 391     if(input_two_points_ids.size()!=2) cout << "Error. Must be two ids in variable input_two_points_ids="<<input_two_points_ids.size()<<endl;
 392   TUIDS::const_iterator aInPointsIter = input_two_points_ids.begin();
 393   for(int i=0;i<2 && aInPointsIter!=input_two_points_ids.end();aInPointsIter++,i++){
 394     current_points_ids[i] = *aInPointsIter;
 395     if (MYDEBUG_REMOVE) cout << " " << *aInPointsIter;
 396   }
 397   if (MYDEBUG_REMOVE) cout << endl;
 398   bool iscurrent_points_changed = false;
 399   points->GetPoint(current_points_ids[0],P[0]);
 400   points->GetPoint(current_points_ids[1],P[1]);
 401   TUIDS::iterator aPointsIter = input_points_ids.begin();
 402   for(;aPointsIter!=input_points_ids.end();aPointsIter++){
 403     if(iscurrent_points_changed){
 404       points->GetPoint(current_points_ids[0],P[0]);
 405       points->GetPoint(current_points_ids[1],P[1]);
 406       iscurrent_points_changed = false;
 407       if (MYDEBUG_REMOVE)
 408         cout << " " << current_points_ids[0] << " " << current_points_ids[1] << endl;
 409     }
 410     // check: is point on line input_two_points_ids
 411     points->GetPoint(*aPointsIter,P[2]);
 412     if (MYDEBUG_REMOVE) {
 413       cout << "\t" << current_points_ids[0] << ":"<<P[0][0]<<","<<P[0][1]<<","<<P[0][2]<<endl;
 414       cout << "\t" << current_points_ids[1] << ":"<<P[1][0]<<","<<P[1][1]<<","<<P[1][2]<<endl;
 415       cout << "\t" << *aPointsIter << ":"<<P[2][0]<<","<<P[2][1]<<","<<P[2][2]<<endl;
 416     }
 417
 418     // x-x1=(x2-x1)*t -> coeff[0][0] = (x-x1), coeff[0][1] = x2-x1
 419     // y-y1=(y2-y1)*t -> coeff[1][0] = (y-y1), coeff[1][1] = y2-y1
 420     // z-z1=(z2-z1)*t -> coeff[2][0] = (z-z1), coeff[2][1] = z2-z1
 421     float coeff[3][2];
 422     for(int i=0;i<3;i++){
 423       coeff[i][0] = P[2][i]-P[0][i];
 424       coeff[i][1] = P[1][i]-P[0][i];
 425     }
 426     bool isok_coord[3];
 427     bool isok = true;
 428     float t[3];
 429     for(int i=0;i<3;i++){
 430       isok_coord[i] = false;
 431       if( fabs(coeff[i][0]) <= EPS && fabs(coeff[i][1]) <= EPS) {
 432         isok_coord[i] = true;
 433         continue;
 434       }
 435       if( fabs(coeff[i][1]) <= EPS && fabs(coeff[i][0]) > EPS) {isok = false;t[i]=1.0/EPS;break;}
 436       t[i] = (coeff[i][0])/(coeff[i][1]);
 437     }
 438     for(int i=0;i<3;i++)
 439       if (MYDEBUG_REMOVE)
 440         cout << __LINE__ << "          "
 441              << coeff[i][0] << ","<<coeff[i][1]
 442              <<"   t="<<t[i]<<" isok_coord="<<isok_coord[i]<<endl;
 443     if(!isok) continue;
 444
 445     if (!isok_coord[0] && !isok_coord[1]){
 446       if (fabs(t[0]-t[1]) <= EPS_T) isok = true;
 447       else isok = false;
 448     }
 449     if (MYDEBUG_REMOVE) cout << __LINE__ << "          1000 " << isok << endl;
 450     if(!isok) continue;
 451     if (!isok_coord[1] && !isok_coord[2]){
 452       if (fabs(t[1] - t[2]) <= EPS_T) isok = true;
 453       else isok = false;
 454     }
 455     if (MYDEBUG_REMOVE) cout << __LINE__ << "          2000 " << isok << endl;
 456     if(!isok) continue;
 457     if (!isok_coord[0] && !isok_coord[2]){
 458       if (fabs(t[0] - t[2]) <= EPS_T) isok = true;
 459       else isok = false;
 460     }
 461     if (MYDEBUG_REMOVE) cout << __LINE__ << "          3000 " << isok<<endl;
 462     if(!isok) continue;
 463
 464     float param[3]; // parametric coord for P[0],P[1],P[2] <--->t[0],t[1],t[2]
 465     // anilize bounds of line
 466     for(int i=0;i<3;i++){
 467       for(int j=0;j<3;j++)
 468         if(!isok_coord[j]) param[i] = (P[i][j]-P[0][j])/(P[1][j]-P[0][j]);
 469     }
 470     if (MYDEBUG_REMOVE) cout << "Params: " << param[0] << "," << param[1] << "," << param[2] << endl;
 471     vtkIdType imax,imin;
 472     float min,max;
 473     for(vtkIdType i=0;i<3;i++)
 474       if(!isok_coord[i]){
 475         min = param[0];imin=0;
 476         max = param[0];imax=0;
 477         break;
 478       }
 479     for(vtkIdType i=0;i<3;i++){
 480         if(min > param[i]) {min = param[i]; imin=i;}
 481         if(max < param[i]) {max = param[i]; imax=i;}
 482     }
 483     if (MYDEBUG_REMOVE)
 484       cout << "\t min="<<min<<"  max="<<max<<" - "<<"imin="<<imin<<"  imax="<<imax<<endl;
 485     // imin - index of left point
 486     // imax - index of right point
 487
 488     // add id to removed point
 489     vtkIdType rem_id,id1,id2;
 490     for(vtkIdType i=0;i<3;i++)
 491       if (i!=imin && i!=imax) {rem_id = i; break;}
 492
 493     if(rem_id == 0) {
 494       rem_id = current_points_ids[0];
 495       id1=current_points_ids[1];
 496       id2=*aPointsIter;
 497     }
 498     else if (rem_id == 1) {
 499       rem_id = current_points_ids[1];
 500       id1=current_points_ids[0];
 501       id2=*aPointsIter;
 502     }
 503     else if (rem_id == 2) {
 504       rem_id = *aPointsIter;
 505       id1=current_points_ids[0];
 506       id2=current_points_ids[1];
 507     }
 508     if (MYDEBUG_REMOVE)
 509       cout << " " << current_points_ids[0] <<","<<current_points_ids[1]<<"---->"<<id1<<","<<id2<<endl;
 510     if((current_points_ids[0] == id1 && current_points_ids[1] == id2) ||
 511        (current_points_ids[0] == id2 && current_points_ids[1] == id1))
 512       {}
 513     else {
 514       iscurrent_points_changed = true;
 515       current_points_ids[0] = id1;
 516       current_points_ids[1] = id2;
 517     }
 518
 519     removed_points_ids.insert(rem_id);
 520   }
 521   out_two_points_ids.insert(current_points_ids[0]);
 522   out_two_points_ids.insert(current_points_ids[1]);
 523 }
 524
 525 static vtkSmartPointer<vtkConvexPointSet> RemoveAllUnneededPoints(vtkConvexPointSet* convex){
 526   vtkSmartPointer<vtkConvexPointSet> out = vtkConvexPointSet::New();
 527
 528   TUIDS two_points,input_points,out_two_points_ids,removed_points_ids,loc_removed_points_ids;
 529   vtkIdList* aPointIds = convex->GetPointIds();
 530   int numIds = aPointIds->GetNumberOfIds();
 531   if (numIds<2) return out;
 532   TSet good_point_ids;
 533   TSet aLists[numIds-2];
 534   for(int i=0;i<numIds-2;i++){
 535     for(int j=i+1;j<numIds-1;j++){
 536       TPair aPair(i,j);
 537       aLists[i].insert(aPair);
 538     }
 539   }
 540   for(vtkIdType i=0;i<numIds-2;i++){
 541     TUIDS::iterator aRemIter=removed_points_ids.find(i);
 542     if(aRemIter!=removed_points_ids.end()) continue;
 543     TSet::iterator aPairIter=aLists[i].begin();
 544     loc_removed_points_ids.clear();
 545     out_two_points_ids.clear();
 546     input_points.clear();
 547     two_points.clear();
 548     for(;aPairIter!=aLists[i].end();aPairIter++){
 549       vtkIdType aFirId = (*aPairIter).first;
 550       vtkIdType aSecId = (*aPairIter).second;
 551       aRemIter=removed_points_ids.find(aSecId);
 552       if(aRemIter!=removed_points_ids.end()) continue;
 553       TPair aPair1(aFirId,aSecId);
 554       TPair aPair2(aSecId,aFirId);
 555       TSet::iterator aGoodIter=good_point_ids.find(aPair1);
 556       if(aGoodIter!=good_point_ids.end()) continue;
 557       aGoodIter=good_point_ids.find(aPair2);
 558       if(aGoodIter!=good_point_ids.end()) continue;
 559       two_points.insert(aFirId);
 560       two_points.insert(aSecId);
 561       if (MYDEBUG_REMOVE)
 562         cout << "\nInput: " << aFirId<<":"<<aPointIds->GetId(aFirId) << "," << aSecId <<":"<<aPointIds->GetId(aSecId)<< "  --- ";
 563       for(vtkIdType k=aSecId+1;k<numIds;k++) {
 564         input_points.insert(k);
 565         if (MYDEBUG_REMOVE) cout << k<<":"<<aPointIds->GetId(k) << ",";
 566       }
 567       if (MYDEBUG_REMOVE) {
 568         cout << endl;
 569         cout << "\t";
 570         for(TUIDS::iterator aDelIter = loc_removed_points_ids.begin();aDelIter!=loc_removed_points_ids.end();aDelIter++)
 571           cout << *aDelIter<<",";
 572         cout << endl;
 573       }
 574       GetAndRemoveIdsOnOneLine(convex->Points,
 575                                input_points,
 576                                two_points,
 577                                out_two_points_ids,
 578                                loc_removed_points_ids);
 579       TUIDS::iterator aOutIter = out_two_points_ids.begin();
 580       vtkIdType aFirst=*aOutIter;aOutIter++;vtkIdType aSecond=*aOutIter;
 581       TPair aPair(aFirst,aSecond);
 582       good_point_ids.insert(aPair);
 583       if (MYDEBUG_REMOVE){
 584         cout << "Output: ";
 585         TUIDS::iterator aIter = out_two_points_ids.begin();
 586         for(;aIter!=out_two_points_ids.end();aIter++)
 587           cout << *aIter << ",";
 588         cout << " --- ";
 589       }
 590       TUIDS::iterator aDelIter = loc_removed_points_ids.begin();
 591       for(;aDelIter!=loc_removed_points_ids.end();aDelIter++){
 592         removed_points_ids.insert(*aDelIter);
 593         if (MYDEBUG_REMOVE) cout << *aDelIter << ",";
 594       }
 595       if (MYDEBUG_REMOVE) cout << endl;
 596     }
 597   }
 598   if (MYDEBUG_REMOVE) {
 599     cout << "============ Resultat ================" <<endl;
 600     cout << "Removed:";
 601     for(TUIDS::iterator aIter=removed_points_ids.begin();aIter!=removed_points_ids.end();aIter++)
 602       cout << *aIter << ",";
 603     cout << endl;
 604   }
 605
 606   TUIDS result_ids,all_ids;
 607   for(vtkIdType i=0;i<numIds;i++) all_ids.insert(i);
 608   std::set_difference(all_ids.begin(),
 609                       all_ids.end(),
 610                       removed_points_ids.begin(),
 611                       removed_points_ids.end(),
 612                       std::inserter(result_ids,result_ids.begin()));
 613
 614   out->Points->SetNumberOfPoints(result_ids.size());
 615   out->PointIds->SetNumberOfIds(result_ids.size());
 616   int aId=0;
 617   if(MYDEBUG_REMOVE) cout << "Result out:";
 618   for(TUIDS::iterator aIter=result_ids.begin();aIter!=result_ids.end();aIter++,aId++){
 619     float P[3];
 620     convex->Points->GetPoint(*aIter,P);
 621     out->Points->SetPoint(aId,P);
 622     out->PointIds->SetId(aId,aPointIds->GetId(*aIter));
 623     if (MYDEBUG_REMOVE) cout << *aIter << ":" << aPointIds->GetId(*aIter) << " , ";
 624   }
 625   if(MYDEBUG_REMOVE) cout << endl;
 626   out->Modified();
 627   out->Initialize();
 628
 629   return out;
 630 }
 631
 632 void GetPolygonalFaces(vtkUnstructuredGrid* theGrid,int cellId,TCellArray &outputCellArray)
 633 {
 634   if (theGrid->GetCellType(cellId) == VTK_CONVEX_POINT_SET){
 635     // get vtkCell type = VTK_CONVEX_POINT_SET
 636     if(vtkConvexPointSet* convex_in = static_cast<vtkConvexPointSet*>(theGrid->GetCell(cellId))){
 637       vtkSmartPointer<vtkConvexPointSet> convex = RemoveAllUnneededPoints(convex_in);
 638       TCellArray f2points;
 639       float convex_center[3]; // convex center point coorinat
 640       int aNbFaces = convex->GetNumberOfFaces();
 641       int numPts = convex->GetNumberOfPoints();
 642       if(MYDEBUG_REMOVE) cout << "aNbFaces="<<aNbFaces<<endl;
 643       if(MYDEBUG_REMOVE) cout << "numPts="<<numPts<<endl;
 644       TPTOIDS p2faces; // key=pointId , value facesIds set
 645
 646       GetCenter(convex->Points,convex_center);
 647
 648       for (vtkIdType faceId=0; faceId < aNbFaces; faceId++){
 649         vtkCell *aFace = convex->GetFace(faceId);
 650         int numFacePts = aFace->GetNumberOfPoints();
 651         TCell aIds;
 652
 653         int i = 0;
 654         for(i=0;i<numFacePts;i++)
 655           aIds.push_back(aFace->GetPointId(i));
 656
 657         float v_a[3],v_b[3],v_convex2face[3]; // vectors
 658         float *id_0,*id_1,*id_n;
 659         /*=============================================
 660                         ,+- - - -  _
 661                    _   / id_n   |  v_b {id_0,id_n}
 662                   v_b /            _
 663                      /          |  v_a {id_0,id_1}
 664                     /
 665                    /            |
 666                   + id_0
 667                    \
 668                   _ \           |
 669                  v_a \
 670                       \ id_1    |
 671                        "+- - - -
 672
 673          ============================================*/
 674         id_0 = theGrid->GetPoint(aIds[0]);
 675         id_1 = theGrid->GetPoint(aIds[1]);
 676         id_n = theGrid->GetPoint(aIds[numFacePts-1]);
 677
 678         for(i=0;i<3;i++){
 679           v_a[i] = id_1[i] - id_0[i];
 680           v_b[i] = id_n[i] - id_0[i];
 681           v_convex2face[i] = id_0[i] - convex_center[i];
 682         }
 683
 684         if (vtkMath::Determinant3x3(v_a,v_b,v_convex2face) < 0){
 685           ReverseIds(aIds);
 686         }
 687
 688         for(i=0;i<(int)aIds.size();i++){
 689           TUIDS &acell = p2faces[aIds[i]];
 690           acell.insert(faceId);
 691         }
 692
 693         f2points[faceId] = aIds;
 694
 695       }
 696
 697       TPTOIDS face2face;
 698       GetFriends(p2faces,f2points,face2face);
 699
 700       TPTOIDS face2points;
 701
 702       // copy TCellArray::f2points to TPTOIDS::face2points
 703       for(TCellArray::iterator f2points_iter=f2points.begin();
 704           f2points_iter!=f2points.end();
 705           f2points_iter++){
 706
 707         TUIDS tmp;
 708         for(TCell::iterator points_iter=(f2points_iter->second).begin();
 709             points_iter!=(f2points_iter->second).end();
 710             points_iter++)
 711           tmp.insert(*points_iter);
 712
 713         face2points[f2points_iter->first] = tmp;
 714       } // end copy
 715
 716
 717       TPTOIDS new_face2faces; // which connected and in one plane
 718
 719       TPTOIDS::const_iterator aF2FIter = face2face.begin();
 720       for(;aF2FIter!=face2face.end();aF2FIter++){
 721         vtkIdType f_key = aF2FIter->first;
 722         TUIDS &faces = new_face2faces[f_key];
 723         //faces.insert(f_key);
 724         TUIDS f_friends = aF2FIter->second;
 725         TUIDS::const_iterator a_friends_iter = f_friends.begin();
 726         for(;a_friends_iter!=f_friends.end();a_friends_iter++){
 727           vtkIdType friend_id = *a_friends_iter;
 728           if( IsConnectedFacesOnOnePlane(theGrid,f_key,friend_id,
 729                                         (face2points.find(f_key))->second,
 730                                         (face2points.find(friend_id))->second)){
 731             faces.insert(friend_id);
 732           } // end if
 733
 734         } // end a_friends_iter
 735       } // end aF2FIter
 736
 737       if(MYDEBUG)
 738       {
 739         TPTOIDS::const_iterator new_face2face_iter = new_face2faces.begin();
 740         cout << "Connected faces and on plane:" << endl;
 741         for(;new_face2face_iter!=new_face2faces.end();new_face2face_iter++){
 742           cout << "Group ["<<new_face2face_iter->first<<"] :";
 743           TUIDS::const_iterator new_faces_iter = (new_face2face_iter->second).begin();
 744           for(;new_faces_iter!=(new_face2face_iter->second).end();new_faces_iter++)
 745             cout << " " << *new_faces_iter ;
 746           cout << endl;
 747         }
 748       }
 749
 750       TPTOIDS output_newid2face;
 751       TCellArray output_newid2face_v2;
 752       {
 753         TUIDS already_in;
 754         TUIDS already_in_tmp;
 755         int k=0;
 756         TPTOIDS::const_iterator new_face2face_iter = new_face2faces.begin();
 757         for(;new_face2face_iter!=new_face2faces.end();new_face2face_iter++){
 758           if(already_in.find(new_face2face_iter->first) != already_in.end())
 759             continue;
 760           if(new_face2face_iter->second.size() > 1)
 761             continue;
 762
 763           TCell &tmp_v2 = output_newid2face_v2[k];
 764           tmp_v2.push_back(new_face2face_iter->first);
 765           already_in.insert(new_face2face_iter->first);
 766
 767           TUIDS::const_iterator new_faces_iter = (new_face2face_iter->second).begin();
 768           for(;new_faces_iter!=(new_face2face_iter->second).end();new_faces_iter++){
 769             if(already_in.find(*new_faces_iter) != already_in.end()) continue;
 770             already_in.insert(*new_faces_iter);
 771
 772             already_in_tmp.clear();
 773             already_in_tmp.insert(new_face2face_iter->first);
 774
 775             TUIDS &tmp = output_newid2face[k];
 776             GetAllFacesOnOnePlane(new_face2faces,*new_faces_iter,
 777                                   already_in_tmp,tmp_v2);
 778
 779             for(TUIDS::const_iterator aIter=already_in_tmp.begin();
 780                 aIter!=already_in_tmp.end();
 781                 aIter++)
 782               {
 783                 already_in.insert(*aIter);
 784                 tmp.insert(*aIter);
 785               }
 786           }
 787           k++;
 788         }
 789       }
 790
 791       if(MYDEBUG) {
 792         cout << "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
 793         {
 794           TPTOIDS::const_iterator new_face2face_iter = output_newid2face.begin();
 795           for(;new_face2face_iter!=output_newid2face.end();new_face2face_iter++){
 796             cout << "Group ["<<new_face2face_iter->first<<"] :";
 797             TUIDS::const_iterator new_faces_iter = (new_face2face_iter->second).begin();
 798             for(;new_faces_iter!=(new_face2face_iter->second).end();new_faces_iter++)
 799               cout << " " << *new_faces_iter ;
 800             cout << endl;
 801           }
 802         }
 803         cout << "++++++++++++++++++++++++++ +++++++++++++++++++++++++++++++"<<endl;
 804         cout << "++++++++++++++++++++++++++ +++++++++++++++++++++++++++++++"<<endl;
 805         cout << "+++++++++++++++++++++++ ++ ++ ++++++++++++++++++++++++++++"<<endl;
 806         cout << "+++++++++++++++++++++++++   ++++++++++++++++++++++++++++++"<<endl;
 807         cout << "++++++++++++++++++++++++++ +++++++++++++++++++++++++++++++"<<endl;
 808         {
 809           TCellArray::const_iterator new_face2face_iter = output_newid2face_v2.begin();
 810           for(;new_face2face_iter!=output_newid2face_v2.end();new_face2face_iter++){
 811             cout << "Group ["<<new_face2face_iter->first<<"] :";
 812             TCell::const_iterator new_faces_iter = (new_face2face_iter->second).begin();
 813             for(;new_faces_iter!=(new_face2face_iter->second).end();new_faces_iter++)
 814               cout << " " << *new_faces_iter ;
 815             cout << endl;
 816           }
 817         }
 818       }
 819       TCellArray output_new_face2ids;
 820 //       {
 821 //      TPTOIDS::const_iterator new_face2face_iter = output_newid2face.begin();
 822 //      for(;new_face2face_iter!=output_newid2face.end();new_face2face_iter++){
 823
 824 //        vtkIdType new_faceId = new_face2face_iter->first;
 825 //        TUIDS::const_iterator new_faces_iter = (new_face2face_iter->second).begin();
 826 //        vtkIdType fId0 = *new_faces_iter;
 827 //        TCellArray::const_iterator pIds0_iter = f2points.find(fId0);
 828 //        TCell pIds0 = pIds0_iter->second;
 829 //        TCell &output = output_new_face2ids[new_faceId];
 830 //        new_faces_iter++;
 831 //        if(new_face2face_iter->second.size() > 2 ){}
 832 //        for(;new_faces_iter!=(new_face2face_iter->second).end();new_faces_iter++){
 833
 834 //          vtkIdType faceId = *new_faces_iter;
 835 //          // find how much nodes in face (f2points)
 836 //          TCellArray::const_iterator pIds_iter = f2points.find(faceId);
 837 //          TCell pIds = pIds_iter->second;
 838
 839 //          GetSumm(pIds0,pIds,output);
 840 //          pIds0 = output;
 841
 842 //        } // end new_faces_iter
 843
 844 //      } // new_face2face_iter
 845 //       }
 846
 847       {
 848         TCellArray::const_iterator new_face2face_iter = output_newid2face_v2.begin();
 849         for(;new_face2face_iter!=output_newid2face_v2.end();new_face2face_iter++){
 850
 851           vtkIdType new_faceId = new_face2face_iter->first;
 852           TCell::const_iterator new_faces_iter = (new_face2face_iter->second).begin();
 853           vtkIdType fId0 = *new_faces_iter;
 854           TCellArray::const_iterator pIds0_iter = f2points.find(fId0);
 855           TCell pIds0 = pIds0_iter->second;
 856           TCell &output = output_new_face2ids[new_faceId];
 857           new_faces_iter++;
 858           if(new_face2face_iter->second.size() == 1 ){
 859             TCellArray::const_iterator pIds_iter = f2points.find(fId0);
 860             TCell pIds = pIds_iter->second;
 861             output = pIds;
 862             continue;
 863           }
 864           for(;new_faces_iter!=(new_face2face_iter->second).end();new_faces_iter++){
 865
 866             vtkIdType faceId = *new_faces_iter;
 867             // find how much nodes in face (f2points)
 868             TCellArray::const_iterator pIds_iter = f2points.find(faceId);
 869             TCell pIds = pIds_iter->second;
 870
 871             GetSumm(pIds0,pIds,output);
 872             pIds0 = output;
 873
 874           } // end new_faces_iter
 875
 876         } // new_face2face_iter
 877       }
 878
 879       if(MYDEBUG) {
 880         cout << "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
 881         cout << "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
 882         cout << "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
 883       }
 884       outputCellArray = output_new_face2ids;//f2points;
 885     }
 886   } else {
 887     // not implemented
 888   }
 889 }
 890 }