tools/vtkEDFOverloads/vtkEDFCutter.cxx

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