tools/vtkEDFOverloads/vtkEDFCutter.cxx

   1 // Copyright (C) 2010-2016  CEA/DEN, EDF R&D, OPEN CASCADE
   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, or (at your option) any later version.
   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/ or email : webmaster.salome@opencascade.com
  18 //
  19
  20 #include "vtkEDFCutter.h"
  21
  22 #include "vtkInformationVector.h"
  23 #include "vtkInformation.h"
  24 #include "vtkSmartPointer.h"
  25 #include "vtkGenericCell.h"
  26 #include "vtkPolyData.h"
  27 #include "vtkObjectFactory.h"
  28 #include "vtkIdTypeArray.h"
  29 #include "vtkCellData.h"
  30 #include "vtkCellArray.h"
  31 #include "vtkIdList.h"
  32
  33 #include <list>
  34 #include <set>
  35 #include <map>
  36 #include <deque>
  37
  38 class vtkEDFEdge
  39 {
  40 public :
  41   vtkIdType v0;
  42   vtkIdType v1;
  43
  44   vtkEDFEdge(vtkIdType a, vtkIdType b) : v0(a), v1(b){}
  45   vtkEDFEdge(){}
  46 };
  47
  48 bool operator == (const vtkEDFEdge& e0, const vtkEDFEdge& e1)
  49 {
  50   return (e0.v0 == e1.v0 && e0.v1 == e1.v1) ||
  51       (e0.v1 == e1.v0 && e0.v0 == e1.v1);
  52 }
  53
  54 bool operator != (const vtkEDFEdge& e0, const vtkEDFEdge& e1)
  55 {
  56   return !(e0==e1);
  57 }
  58
  59 bool operator < (const vtkEDFEdge& e0, const vtkEDFEdge& e1)
  60 {
  61   vtkEDFEdge the_e0;
  62   vtkEDFEdge the_e1;
  63   if(e0.v0 < e0.v1)
  64     {
  65     the_e0.v0 = e0.v0;
  66     the_e0.v1 = e0.v1;
  67     }
  68   else
  69     {
  70     the_e0.v0 = e0.v1;
  71     the_e0.v1 = e0.v0;
  72     }
  73   if(e1.v0 < e1.v1)
  74     {
  75     the_e1.v0 = e1.v0;
  76     the_e1.v1 = e1.v1;
  77     }
  78   else
  79     {
  80     the_e1.v0 = e1.v1;
  81     the_e1.v1 = e1.v0;
  82     }
  83
  84   if(the_e0.v0 == the_e1.v0)
  85     return (the_e0.v1 < the_e1.v1);
  86
  87   return the_e0.v0 < the_e1.v0;
  88 }
  89
  90 vtkStandardNewMacro(vtkEDFCutter);
  91
  92 vtkEDFCutter::vtkEDFCutter()
  93 {
  94   this->OriginalCellIdsName = NULL;
  95 }
  96
  97 vtkEDFCutter::~vtkEDFCutter()
  98 {
  99   this->SetOriginalCellIdsName(NULL);
 100 }
 101
 102 int vtkEDFCutter::RequestData(vtkInformation * request,
 103                         vtkInformationVector ** inVector,
 104                         vtkInformationVector * outVector)
 105 {
 106   // get the info objects
 107   vtkInformation *inInfo = inVector[0]->GetInformationObject(0);
 108   vtkInformation *outInfo = outVector->GetInformationObject(0);
 109
 110   // get the input and output
 111   vtkDataSet *input = vtkDataSet::SafeDownCast(
 112     inInfo->Get(vtkDataObject::DATA_OBJECT()));
 113
 114   vtkSmartPointer<vtkIdTypeArray> cellIdArray =
 115       vtkSmartPointer<vtkIdTypeArray>::New();
 116   cellIdArray->SetName(this->GetOriginalCellIdsName());
 117   cellIdArray->SetNumberOfComponents(1);
 118   cellIdArray->SetNumberOfTuples(input->GetNumberOfCells());
 119   for(vtkIdType id=0; id < cellIdArray->GetNumberOfTuples(); id++)
 120     {
 121     cellIdArray->SetTuple1(id, id);
 122     }
 123   input->GetCellData()->AddArray(cellIdArray);
 124
 125   int ret = this->Superclass::RequestData(request, inVector, outVector);
 126
 127   if(ret == 0)
 128     return 0;
 129
 130   vtkPolyData *output = vtkPolyData::SafeDownCast(
 131     outInfo->Get(vtkDataObject::DATA_OBJECT()));
 132
 133   vtkSmartPointer<vtkPolyData> tmpOutput;
 134   tmpOutput.TakeReference(output->NewInstance());
 135
 136   this->AssembleOutputTriangles(output, tmpOutput);
 137
 138   output->ShallowCopy(tmpOutput);
 139
 140   return ret;
 141 }
 142
 143
 144 void  vtkEDFCutter::AssembleOutputTriangles(vtkPolyData* inpd,
 145                                             vtkPolyData* outpd)
 146 {
 147   outpd->ShallowCopy(inpd);
 148
 149   vtkIdTypeArray* originalCellIds = vtkIdTypeArray::SafeDownCast(
 150       inpd->GetCellData()->GetArray(this->GetOriginalCellIdsName()));
 151
 152   if(originalCellIds == NULL)
 153     {
 154     return;
 155     }
 156
 157   outpd->GetCellData()->Initialize();
 158   outpd->GetCellData()->CopyAllocate(inpd->GetCellData());
 159
 160   vtkSmartPointer<vtkCellArray> verts = vtkSmartPointer<vtkCellArray>::New();
 161   vtkSmartPointer<vtkCellArray> lines = vtkSmartPointer<vtkCellArray>::New();
 162   vtkSmartPointer<vtkCellArray> polys = vtkSmartPointer<vtkCellArray>::New();
 163   vtkSmartPointer<vtkCellArray> strips = vtkSmartPointer<vtkCellArray>::New();
 164   outpd->SetVerts(verts);
 165   outpd->SetLines(lines);
 166   outpd->SetPolys(polys);
 167   outpd->SetStrips(strips);
 168
 169   for(vtkIdType cellId=0; cellId<inpd->GetNumberOfCells(); cellId++)
 170     {
 171     unsigned char type = inpd->GetCellType(cellId);
 172     if(type != VTK_TRIANGLE)
 173       {
 174       vtkIdType npts;
 175       vtkIdType* pts = NULL;
 176       inpd->GetCellPoints(cellId, npts, pts);
 177       vtkIdType newCellId =
 178           outpd->InsertNextCell(type, npts, pts);
 179       outpd->GetCellData()->CopyData(inpd->GetCellData(), cellId, newCellId);
 180       }
 181     else
 182       {
 183       vtkIdType cellIdEnd = cellId+1;
 184       vtkIdType originalCellId = originalCellIds->GetValue(cellId);
 185       while(cellIdEnd < inpd->GetNumberOfCells() &&
 186             inpd->GetCellType(cellIdEnd) == VTK_TRIANGLE &&
 187             originalCellIds->GetValue(cellIdEnd) == originalCellId)
 188         {
 189         cellIdEnd++;
 190         }
 191
 192       // all triangles from cellId to cellIdEnd come from the same
 193       // original cell.
 194
 195       // A batch is composed of triangles which are connected by the edges.
 196       std::map<vtkIdType, std::set<vtkIdType> > connectedTriangles;
 197       for(vtkIdType firstCell = cellId; firstCell < cellIdEnd-1; firstCell++)
 198         {
 199         vtkIdType npts;
 200         vtkIdType* pts = NULL;
 201         inpd->GetCellPoints(firstCell, npts, pts);
 202         vtkEDFEdge fe0 = vtkEDFEdge(pts[0], pts[1]);
 203         vtkEDFEdge fe1 = vtkEDFEdge(pts[1], pts[2]);
 204         vtkEDFEdge fe2 = vtkEDFEdge(pts[2], pts[0]);
 205         for(vtkIdType secondCell = firstCell+1; secondCell < cellIdEnd; secondCell++)
 206           {
 207           vtkIdType snpts;
 208           vtkIdType* spts = NULL;
 209           inpd->GetCellPoints(secondCell, snpts, spts);
 210           vtkEDFEdge se0 = vtkEDFEdge(spts[0], spts[1]);
 211           vtkEDFEdge se1 = vtkEDFEdge(spts[1], spts[2]);
 212           vtkEDFEdge se2 = vtkEDFEdge(spts[2], spts[0]);
 213
 214           if(fe0 == se0 || fe0 == se1 || fe0 == se2 ||
 215              fe1 == se0 || fe1 == se1 || fe1 == se2 ||
 216              fe2 == se0 || fe2 == se1 || fe2 == se2)
 217             {
 218             connectedTriangles[firstCell].insert(secondCell);
 219             connectedTriangles[secondCell].insert(firstCell);
 220             }
 221           }
 222         }
 223
 224       std::set<vtkIdType> visitedCell;
 225       for(vtkIdType id=cellId; id<cellIdEnd; id++)
 226         {
 227         if(visitedCell.find(id) != visitedCell.end())
 228           continue;
 229
 230         // if this cell has not yet been visited, I create a batch of all
 231         // cells connected to this one
 232
 233         visitedCell.insert(id);
 234         std::set<vtkIdType> batch;
 235         std::list<vtkIdType> cellList;
 236         cellList.push_back(id);
 237         while(cellList.size() > 0)
 238           {
 239           vtkIdType currentId = *(cellList.begin());
 240           batch.insert(currentId);
 241           cellList.pop_front();
 242           if(connectedTriangles.find(currentId) != connectedTriangles.end())
 243             {
 244             const std::set<vtkIdType>& adj = connectedTriangles[currentId];
 245             std::set<vtkIdType>::const_iterator it = adj.begin();
 246             while(it != adj.end())
 247               {
 248               vtkIdType other = *it;
 249               if(visitedCell.find(other) == visitedCell.end())
 250                 {
 251                 cellList.push_back(other);
 252                 visitedCell.insert(other);
 253                 }
 254               it++;
 255               }
 256             }
 257           }
 258
 259
 260
 261         // then I add this batch to the output,
 262         // creating a unique cell for the whole batch.
 263
 264         if(batch.size() == 1)
 265           {
 266           vtkIdType tid = *(batch.begin());
 267           vtkIdType npts;
 268           vtkIdType* pts = NULL;
 269           inpd->GetCellPoints(tid, npts, pts);
 270           vtkIdType newCellId =
 271               outpd->InsertNextCell(VTK_TRIANGLE, npts, pts);
 272           outpd->GetCellData()->CopyData(inpd->GetCellData(), cellId, newCellId);
 273           }
 274         else if(batch.size() == 2)
 275           { // two triangles connected together --> create a VTK_QUAD
 276           vtkIdType fid = *(batch.begin());
 277           vtkIdType sid = *(batch.rbegin());
 278           vtkIdType fnpts;
 279           vtkIdType* fpts = NULL;
 280           inpd->GetCellPoints(fid, fnpts, fpts);
 281           vtkIdType snpts;
 282           vtkIdType* spts = NULL;
 283           inpd->GetCellPoints(sid, snpts, spts);
 284
 285           int findex = 0;
 286           vtkIdType fv = fpts[findex];
 287           while(((fv == spts[0]) ||
 288                  (fv == spts[1]) ||
 289                  (fv == spts[2])) && findex < 3)
 290             {
 291             findex++;
 292             fv = fpts[findex];
 293             }
 294           if(findex == 3)
 295             {// this is a degenerate case : one of the triangles use
 296             // only 2 vertices
 297             findex = 0;
 298             }
 299           int sindex = 0;
 300           vtkIdType sv = spts[sindex];
 301           while(((sv == fpts[0]) ||
 302                  (sv == fpts[1]) ||
 303                  (sv == fpts[2])) && sindex < 3)
 304             {
 305             sindex++;
 306             sv = spts[sindex];
 307             }
 308           if(sindex == 3)
 309             {// this is a degenerate case : one of the triangles use
 310             // only 2 vertices
 311             sindex = 0;
 312             }
 313
 314           vtkIdType pts[4];
 315           pts[0] = fpts[findex];
 316           pts[1] = fpts[(findex+1)%3];
 317           pts[2] = spts[sindex];
 318           pts[3] = fpts[(findex+2)%3];
 319
 320           vtkIdType newCellId = outpd->InsertNextCell(VTK_QUAD, 4, pts);
 321           outpd->GetCellData()->CopyData(inpd->GetCellData(), cellId, newCellId);
 322           }
 323         else
 324           {
 325           std::deque<vtkEDFEdge> contour;
 326
 327           std::list<vtkIdType> toVisit;
 328           std::set<vtkIdType> visited;
 329
 330           toVisit.push_back(*(batch.begin()));
 331
 332           std::set<vtkIdType> triedAgain;
 333
 334           while(toVisit.size() > 0)
 335             {
 336             vtkIdType currentId = *(toVisit.begin());
 337             toVisit.pop_front();
 338             if(visited.find(currentId) != visited.end())
 339               continue;
 340
 341             visited.insert(currentId);
 342             const std::set<vtkIdType>& adj = connectedTriangles[currentId];
 343             std::set<vtkIdType>::const_iterator it = adj.begin();
 344             while(it != adj.end())
 345               {
 346               vtkIdType adjid = *it;
 347               it++;
 348               if(visited.find(adjid) != visited.end())
 349                 continue;
 350
 351               toVisit.push_back(adjid);
 352               }
 353
 354             vtkIdType npts;
 355             vtkIdType* pts = NULL;
 356             inpd->GetCellPoints(currentId, npts, pts);
 357             vtkEDFEdge edge[3] = {vtkEDFEdge(pts[0], pts[1]),
 358                                   vtkEDFEdge(pts[1], pts[2]),
 359                                   vtkEDFEdge(pts[2], pts[0])};
 360
 361             // special case : initialization of the contour
 362             if(contour.size() == 0)
 363               {
 364               contour.push_back(edge[0]);
 365               contour.push_back(edge[1]);
 366               contour.push_back(edge[2]);
 367               continue;
 368               }
 369
 370             // Find which edge of the contour
 371             // is connected to the current triangle
 372             int orient = 0;
 373             std::deque<vtkEDFEdge>::iterator contourit = contour.begin();
 374             bool found = false;
 375             while(contourit != contour.end())
 376               {
 377               vtkEDFEdge& e = *contourit;
 378               for(orient = 0; orient<3; orient++)
 379                 {
 380                 if(e == edge[orient])
 381                   {
 382                   found = true;
 383                   break;
 384                   }
 385                 }
 386               if(found)
 387                 break;
 388
 389               contourit++;
 390               }
 391             if(contourit == contour.end())
 392               {// this triangle is not connected to the current contour.
 393               // put it back in the queue for later processing
 394               if(triedAgain.find(currentId) == triedAgain.end())
 395                 {
 396                 triedAgain.insert(currentId);
 397                 toVisit.push_back(currentId);
 398                 visited.erase(currentId);
 399                 continue;
 400                 }
 401               else
 402                 {
 403                 vtkDebugMacro( << "triangle " << currentId
 404                   << "could not be added to the contour of the current batch");
 405                 continue;
 406                 }
 407               }
 408             // if I reach this point, I will add the triangle to the contour
 409             // --> the contour will be modified and I can try again
 410             // to add the previously rejected triangles
 411             triedAgain.clear();
 412
 413             // Now, merge the edges of the current triangle with
 414             // the contour
 415             vtkEDFEdge& tbeforeedge = edge[(orient+1)%3];
 416             vtkEDFEdge& tafteredge = edge[(orient+2)%3];
 417
 418             std::deque<vtkEDFEdge>::iterator beforeit;
 419             if(contourit == contour.begin())
 420               beforeit = contour.end()-1;
 421             else
 422               beforeit = contourit - 1;
 423             vtkEDFEdge& beforeedge = *beforeit;
 424
 425             std::deque<vtkEDFEdge>::iterator afterit;
 426             if(contourit == contour.end()-1)
 427               afterit = contour.begin();
 428             else
 429               afterit = contourit + 1;
 430             vtkEDFEdge& afteredge = *afterit;
 431
 432             if(beforeedge == tbeforeedge)
 433               {
 434               if(afteredge == tafteredge)
 435                 {// in this case, I am adding a triangle that is fully inside
 436                 // the contour. I need to remove the three edges from the
 437                 // contour.
 438                 if(contour.size() == 3)
 439                   {
 440                   contour.clear();
 441                   }
 442                 else
 443                   {
 444                   std::deque<vtkEDFEdge>::iterator lastit;
 445                   if(afterit == contour.end()-1)
 446                     lastit = contour.begin();
 447                   else
 448                     lastit = afterit + 1;
 449
 450                   if(lastit < beforeit)
 451                     {
 452                     contour.erase(beforeit, contour.end());
 453                     contour.erase(contour.begin(), lastit);
 454                     }
 455                   else
 456                     {
 457                     contour.erase(beforeit, lastit);
 458                     }
 459                   }
 460                 }
 461               else
 462                 {// the edge before is the glued, remove the two adjacent edges
 463                 // and add the edge after
 464                 if(beforeit == contour.end()-1)
 465                   {
 466                   contour.erase(beforeit, contour.end());
 467                   contour.erase(contour.begin(), contour.begin()+1);
 468                   contour.push_back(tafteredge);
 469                   }
 470                 else
 471                   {
 472                   int index = beforeit - contour.begin();
 473                   contour.erase(beforeit, contourit+1);
 474                   contour.insert(contour.begin()+index, tafteredge);
 475                   }
 476                 }
 477               }
 478             else if(afteredge == tafteredge)
 479               {// the edge after is glued, removed the two glued edges and add
 480               // the edge new edge
 481               if(contourit == contour.end() -1)
 482                 {
 483                 contour.erase(contour.end() - 1);
 484                 contour.erase(contour.begin());
 485                 contour.push_back(tbeforeedge);
 486                 }
 487               else
 488                 {
 489                 int index = contourit - contour.begin();
 490                 contour.erase(contourit, afterit+1);
 491                 contour.insert(contour.begin()+index, tbeforeedge);
 492                 }
 493               }
 494             else
 495               {
 496               int index = contourit - contour.begin();
 497               contour.erase(contourit);
 498               contour.insert(contour.begin()+index, tbeforeedge);
 499               contour.insert(contour.begin()+index+1, tafteredge);
 500               }
 501             }
 502           vtkSmartPointer<vtkIdList> ids = vtkSmartPointer<vtkIdList>::New();
 503           std::deque<vtkEDFEdge>::iterator cit = contour.begin();
 504           while(cit != contour.end())
 505             {
 506             vtkEDFEdge& e = *cit;
 507             cit++;
 508             ids->InsertNextId(e.v0);
 509             }
 510
 511           vtkIdType newCellId = outpd->InsertNextCell(VTK_POLYGON, ids);
 512           outpd->GetCellData()->CopyData(inpd->GetCellData(), cellId, newCellId);
 513           }
 514         }
 515       cellId = cellIdEnd - 1;
 516       }
 517     }
 518 }
 519
 520 void  vtkEDFCutter::PrintSelf(ostream& os, vtkIndent indent)
 521 {
 522   this->Superclass::PrintSelf(os, indent);
 523 }
 524
 525