src/PIPELINE/VISU_LabelPointsFilter.cxx

   1 //  Copyright (C) 2007-2008  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 //  Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 //  CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 //  This library is free software; you can redistribute it and/or
   7 //  modify it under the terms of the GNU Lesser General Public
   8 //  License as published by the Free Software Foundation; either
   9 //  version 2.1 of the License.
  10 //
  11 //  This library is distributed in the hope that it will be useful,
  12 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 //  Lesser General Public License for more details.
  15 //
  16 //  You should have received a copy of the GNU Lesser General Public
  17 //  License along with this library; if not, write to the Free Software
  18 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 //  See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22 //  VISU OBJECT : interactive object for VISU entities implementation
  23 // File:    VISU_LabelPointsFilter.cxx
  24 // Author:  Vitaly Smetannikov
  25 // Module : VISU
  26 //
  27 #include "VISU_LabelPointsFilter.hxx"
  28
  29 #include <vtkPolyData.h>
  30 #include <vtkObjectFactory.h>
  31 #include <vtkInformation.h>
  32 #include <vtkInformationVector.h>
  33 #include <vtkFeatureEdges.h>
  34 #include <vtkCellArray.h>
  35 #include <vtkPointData.h>
  36 #include <vtkCellData.h>
  37 #include <vtkDataArray.h>
  38 #include <vtkDoubleArray.h>
  39 #include <vtkGeometryFilter.h>
  40 #include <vtkPolyDataConnectivityFilter.h>
  41 #include <vtkMath.h>
  42
  43 #include <set>
  44 #include <vector>
  45
  46
  47 #define CONTAINS(SET, PT) (SET.find(PT) != SET.end())
  48
  49 struct ltIdType
  50 {
  51   bool operator()(const vtkIdType a1, const vtkIdType a2) const
  52   {
  53     return a1 < a2;
  54   }
  55 };
  56
  57
  58
  59 //----------------------------------------------------------------------------
  60 vtkStandardNewMacro(VISU_LabelPointsFilter);
  61
  62
  63 //----------------------------------------------------------------------------
  64 void VISU_LabelPointsFilter::SetPointsNb(int theNb)
  65 {
  66   if (myPointsNb == theNb) return;
  67   myPointsNb = (theNb < 1)? 1:theNb;
  68   Modified();
  69 }
  70
  71 //----------------------------------------------------------------------------
  72 VISU_LabelPointsFilter::VISU_LabelPointsFilter():
  73   vtkPolyDataAlgorithm(),
  74   myPointsNb(3)
  75 {
  76 }
  77
  78 //----------------------------------------------------------------------------
  79 VISU_LabelPointsFilter::~VISU_LabelPointsFilter()
  80 {}
  81
  82
  83
  84 //----------------------------------------------------------------------------
  85 int VISU_LabelPointsFilter::RequestData(vtkInformation* vtkNotUsed(request),
  86                                         vtkInformationVector** inputVector,
  87                                         vtkInformationVector* outputVector)
  88 {
  89   // get the info objects
  90   vtkInformation* inInfo = inputVector[0]->GetInformationObject(0);
  91   vtkInformation* outInfo = outputVector->GetInformationObject(0);
  92
  93   // get the input and ouptut
  94   vtkPolyData* input = vtkPolyData::SafeDownCast(inInfo->Get(vtkDataObject::DATA_OBJECT()));
  95   vtkPolyData* output = vtkPolyData::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT()));
  96
  97   GetRegions(input, output);
  98
  99   return 1;
 100 }
 101
 102 int VISU_LabelPointsFilter::RequestUpdateExtent(vtkInformation* vtkNotUsed(request),
 103                                                 vtkInformationVector** inputVector,
 104                                                 vtkInformationVector* outputVector)
 105 {
 106   // get the info objects
 107   vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
 108   vtkInformation *outInfo = outputVector->GetInformationObject(0);
 109
 110   vtkPolyData* input = vtkPolyData::SafeDownCast(inInfo->Get(vtkDataObject::DATA_OBJECT()));
 111   vtkPolyData* output = vtkPolyData::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT()));
 112
 113   GetRegions(input, output);
 114   return 1;
 115 }
 116
 117
 118
 119 int VISU_LabelPointsFilter::GetRegions(vtkPolyData* theInput,
 120                                        vtkPolyData* theOutput)
 121 {
 122   vtkIdType cellId, i;
 123   vtkIdType numPts, numCells;
 124   vtkPoints *inPts;
 125
 126   vtkPointData *aInpPD=theInput->GetPointData(), *aOutPD=theOutput->GetPointData();
 127   vtkCellData *aInpCD=theInput->GetCellData(), *aOutCD=theOutput->GetCellData();
 128
 129   //  Check input/allocate storage
 130   //
 131   inPts = theInput->GetPoints();
 132
 133   if (inPts == NULL)
 134     return 1;
 135
 136   numPts = inPts->GetNumberOfPoints();
 137   numCells = theInput->GetNumberOfCells();
 138
 139   if ( numPts < 1 || numCells < 1 )
 140     return 1;
 141
 142   // Build cell structure
 143   //
 144   vtkPolyData* aMesh = vtkPolyData::New();
 145   aMesh->CopyStructure(theInput);
 146   aMesh->BuildLinks();
 147
 148   // Initialize.  Keep track of points and cells visited.
 149   //
 150   vtkIdTypeArray* aRegionSizes = vtkIdTypeArray::New();
 151   int* aVisited = new int[numCells];
 152   for ( i=0; i < numCells; i++ )
 153      aVisited[i] = -1;
 154
 155   vtkIdType* aPointMap = new vtkIdType[numPts];
 156   for ( i=0; i < numPts; i++ )
 157     aPointMap[i] = -1;
 158
 159   vtkPoints* newPts = vtkPoints::New();
 160   newPts->Allocate(numPts);
 161
 162   // Traverse all cells marking those visited.  Each new search
 163   // starts a new connected region. Connected region grows
 164   // using a connected wave propagation.
 165   //
 166   vtkIdList* aWave = vtkIdList::New();
 167   aWave->Allocate(numPts/4+1,numPts);
 168   vtkIdList* aWave2 = vtkIdList::New();
 169   aWave2->Allocate(numPts/4+1,numPts);
 170
 171   vtkIdType aPointNumber = 0;
 172   int aRegionNumber = 0;
 173
 174   vtkIdList* aCellIds = vtkIdList::New();
 175   aCellIds->Allocate(8, VTK_CELL_SIZE);
 176   vtkIdList* aPointIds = vtkIdList::New();
 177   aPointIds->Allocate(8, VTK_CELL_SIZE);
 178
 179   //  vtkIdType aNumCellsInRegion;
 180
 181   aOutPD->CopyAllocate(aInpPD);
 182   aOutCD->CopyAllocate(aInpCD);
 183
 184   //visit all cells marking with region number
 185   for (cellId=0; cellId < numCells; cellId++) {
 186     if ( aVisited[cellId] < 0 ) {
 187       aWave->InsertNextId(cellId);
 188       aPointNumber = 0;
 189       TraverseAndMark(aWave, aWave2, aVisited, aPointMap,
 190                       aRegionNumber, aPointNumber, aMesh);
 191
 192       if (aPointNumber >= myPointsNb) {
 193         std::set<vtkIdType, ltIdType> aIdxSet;
 194         for (i=0; i < numPts; i++) {
 195           if ( aPointMap[i] > -1 ) {
 196             aIdxSet.insert(i);
 197             aPointMap[i] = -1;
 198           }
 199         }
 200         std::vector<vtkIdType> aIdx(aIdxSet.begin(), aIdxSet.end());
 201         int aActualPts = aIdx.size();
 202         int aNewId;
 203         if (myPointsNb > 2) {
 204           int k = aActualPts/(myPointsNb - 1);
 205           int count;
 206           for (i=0, count = 0; i < aActualPts; i+=k, count++) {
 207             aNewId = newPts->InsertNextPoint(inPts->GetPoint(aIdx[i]));
 208             aOutPD->CopyData(aInpPD, aIdx[i], aNewId);
 209           }
 210           if (count < myPointsNb) {
 211             aNewId = newPts->InsertNextPoint(inPts->GetPoint(aIdx[aActualPts - 1]));
 212             aOutPD->CopyData(aInpPD, aIdx[aActualPts - 1], aNewId);
 213           }
 214         } else {
 215           aNewId = newPts->InsertNextPoint(inPts->GetPoint(aIdx[0]));
 216           aOutPD->CopyData(aInpPD, aIdx[0], aNewId);
 217           if (myPointsNb == 2) {
 218             aNewId = newPts->InsertNextPoint(inPts->GetPoint(aIdx[aActualPts - 1]));
 219             aOutPD->CopyData(aInpPD, aIdx[aActualPts - 1], aNewId);
 220           }
 221         }
 222       }
 223       aWave->Reset();
 224       aWave2->Reset();
 225     }
 226   }
 227
 228   aWave->Delete();
 229   aWave2->Delete();
 230
 231   theOutput->SetPoints(newPts);
 232   newPts->Delete();
 233
 234
 235   delete [] aVisited;
 236   delete [] aPointMap;
 237   aMesh->Delete();
 238   theOutput->Squeeze();
 239   aCellIds->Delete();
 240   aPointIds->Delete();
 241
 242   return aRegionSizes->GetMaxId() + 1;
 243 }
 244
 245
 246 // Mark current cell as visited and assign region number.  Note:
 247 // traversal occurs across shared vertices.
 248 //
 249 void VISU_LabelPointsFilter::TraverseAndMark (vtkIdList* theWave,
 250                                               vtkIdList* theWave2,
 251                                               int* theVisited,
 252                                               vtkIdType* thePointMap,
 253                                               int& theRegionNumber,
 254                                               vtkIdType& thePointNumber,
 255                                               vtkPolyData* theMesh)
 256 {
 257   vtkIdType cellId, ptId, numIds, i;
 258   int j, k;
 259   vtkIdType *pts, *cells, npts;
 260   vtkIdList *tmpWave;
 261   unsigned short ncells;
 262   vtkIdList* aNeighborCellPointIds = vtkIdList::New();
 263
 264
 265   while ( (numIds=theWave->GetNumberOfIds()) > 0 ) {
 266     for ( i=0; i < numIds; i++ ) {
 267       cellId = theWave->GetId(i);
 268       if ( theVisited[cellId] < 0 ) {
 269         theVisited[cellId] = theRegionNumber;
 270         theMesh->GetCellPoints(cellId, npts, pts);
 271
 272         for (j=0; j < npts; j++) {
 273           if ( thePointMap[ptId=pts[j]] < 0 ) {
 274             thePointMap[ptId] = thePointNumber++;
 275           }
 276           theMesh->GetPointCells(ptId,ncells,cells);
 277
 278           // check connectivity criterion (geometric + scalar)
 279           for (k=0; k < ncells; k++) {
 280             cellId = cells[k];
 281             theWave2->InsertNextId(cellId);
 282             //              }
 283           }//for all cells using this point
 284         }//for all points of this cell
 285       }//if cell not yet visited
 286     }//for all cells in this wave
 287
 288     tmpWave = theWave;
 289     theWave = theWave2;
 290     theWave2 = tmpWave;
 291     tmpWave->Reset();
 292   } //while wave is not empty
 293 }