src/PIPELINE/VISU_LabelPointsFilter.cxx

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