src/Plugins/ElevationSurface/vtkElevationSurfaceFilter.cxx

   1 // Copyright (C) 2010-2012  CEA/DEN, EDF 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/ or email : webmaster.salome@opencascade.com
  18 //
  19
  20 #include "vtkElevationSurfaceFilter.h"
  21 #include "vtkInformation.h"
  22 #include "vtkInformationVector.h"
  23 #include "vtkObjectFactory.h"
  24 #include "vtkElevationSurfaceFilter.h"
  25 #include "vtkPolyData.h"
  26 #include "vtkIdTypeArray.h"
  27 #include "vtkPolyData.h"
  28 #include "vtkUnstructuredGrid.h"
  29 #include "vtkDemandDrivenPipeline.h"
  30 #include "vtkStreamingDemandDrivenPipeline.h"
  31 #include "vtkGenericCell.h"
  32 #include "vtkSmartPointer.h"
  33 #include "vtkPoints.h"
  34 #include "vtkCellArray.h"
  35 #include "vtkPointData.h"
  36 #include "vtkCellData.h"
  37 //#include "vtkDataSetSurfaceFilter.h"
  38
  39 #include <math.h>
  40
  41 vtkCxxRevisionMacro(vtkElevationSurfaceFilter, "$Revision$");
  42 vtkStandardNewMacro(vtkElevationSurfaceFilter);
  43
  44 vtkElevationSurfaceFilter::vtkElevationSurfaceFilter()
  45 {
  46   this->SetNumberOfInputPorts(1);
  47   this->SetNumberOfOutputPorts(1);
  48
  49   this->ScaleFactor = 0.5;
  50   this->Direction[0] = 0.0;
  51   this->Direction[1] = 0.0;
  52   this->Direction[2] = 1.0;
  53   this->AutoDetectDirection = 1;
  54 }
  55
  56 vtkElevationSurfaceFilter::~vtkElevationSurfaceFilter()
  57 {
  58 }
  59
  60 //----------------------------------------------------------------------------
  61 int vtkElevationSurfaceFilter::ProcessRequest(vtkInformation* request,
  62                                          vtkInformationVector** inputVector,
  63                                          vtkInformationVector* outputVector)
  64 {
  65   // generate the data
  66   if(request->Has(vtkDemandDrivenPipeline::REQUEST_DATA()))
  67     {
  68     return this->RequestData(request, inputVector, outputVector);
  69     }
  70
  71   if(request->Has(vtkStreamingDemandDrivenPipeline::REQUEST_UPDATE_EXTENT()))
  72     {
  73     return this->RequestUpdateExtent(request, inputVector, outputVector);
  74     }
  75
  76   // execute information
  77   if(request->Has(vtkDemandDrivenPipeline::REQUEST_INFORMATION()))
  78     {
  79     return this->RequestInformation(request, inputVector, outputVector);
  80     }
  81
  82   return this->Superclass::ProcessRequest(request, inputVector, outputVector);
  83 }
  84
  85 //----------------------------------------------------------------------------
  86 int vtkElevationSurfaceFilter::FillOutputPortInformation(
  87   int vtkNotUsed(port), vtkInformation* info)
  88 {
  89   // now add our info
  90   info->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkUnstructuredGrid");
  91   return 1;
  92 }
  93
  94 //----------------------------------------------------------------------------
  95 int vtkElevationSurfaceFilter::FillInputPortInformation(
  96   int vtkNotUsed(port), vtkInformation* info)
  97 {
  98   info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkPointSet");
  99   return 1;
 100 }
 101
 102 //----------------------------------------------------------------------------
 103 int vtkElevationSurfaceFilter::RequestUpdateExtent(
 104   vtkInformation* vtkNotUsed(request),
 105   vtkInformationVector** inputVector,
 106   vtkInformationVector* vtkNotUsed(outputVector))
 107 {
 108   int numInputPorts = this->GetNumberOfInputPorts();
 109   for (int i=0; i<numInputPorts; i++)
 110     {
 111     int numInputConnections = this->GetNumberOfInputConnections(i);
 112     for (int j=0; j<numInputConnections; j++)
 113       {
 114       vtkInformation* inputInfo = inputVector[i]->GetInformationObject(j);
 115       inputInfo->Set(vtkStreamingDemandDrivenPipeline::EXACT_EXTENT(), 1);
 116       }
 117     }
 118   return 1;
 119 }
 120
 121 int vtkElevationSurfaceFilter::RequestInformation(vtkInformation *request,
 122     vtkInformationVector **input, vtkInformationVector *output)
 123 {
 124   return 1;
 125 }
 126
 127 int vtkElevationSurfaceFilter::RequestData(vtkInformation *request,
 128     vtkInformationVector **input, vtkInformationVector *output)
 129 {
 130         vtkPointSet *psIn = vtkPointSet::SafeDownCast(
 131       input[0]->GetInformationObject(0)->Get(vtkDataObject::DATA_OBJECT()));
 132
 133   //vtkPolyData *psIn = vtkPolyData::SafeDownCast(
 134   //    input[0]->GetInformationObject(0)->Get(vtkDataObject::DATA_OBJECT()));
 135
 136   /*vtkSmartPointer<vtkDataSetSurfaceFilter> surfaceFilter = vtkSmartPointer<vtkDataSetSurfaceFilter>::New();
 137   surfaceFilter->SetInput(psIn);
 138   vtkPolyData* psIn = vtkPolyData::SafeDownCast(surfaceFilter->GetOutput());*/
 139
 140   vtkUnstructuredGrid *usgOut = vtkUnstructuredGrid::SafeDownCast(
 141       output->GetInformationObject(0)->Get(vtkDataObject::DATA_OBJECT()));
 142
 143   vtkDataArray* array = this->GetInputArrayToProcess(0, input);
 144
 145   if(psIn == NULL || array == NULL || usgOut == NULL
 146      || array->GetNumberOfComponents() != 1)
 147     {
 148     vtkDebugMacro("vtkElevationSurfaceFilter no correctly configured");
 149     return 1;
 150     }
 151
 152   double dir[3];
 153   if(this->AutoDetectDirection)
 154     {
 155     this->ComputeDirection(psIn, dir);
 156     }
 157   else
 158     {
 159     dir[0] = this->Direction[0];
 160     dir[1] = this->Direction[1];
 161     dir[2] = this->Direction[2];
 162     }
 163
 164   double len = dir[0]*dir[0] + dir[1]*dir[1] + dir[2]*dir[2];
 165
 166   if(len == 0)
 167     {
 168     dir[2] = 1;
 169     len = 1.0;
 170     }
 171
 172   len = sqrt(len);
 173
 174   dir[0] /= len;
 175   dir[1] /= len;
 176   dir[2] /= len;
 177
 178   dir[0] *= this->GetScaleFactor();
 179   dir[1] *= this->GetScaleFactor();
 180   dir[2] *= this->GetScaleFactor();
 181
 182   usgOut->Allocate(psIn->GetNumberOfCells());
 183
 184   vtkSmartPointer<vtkPoints> newPts = vtkSmartPointer<vtkPoints>::New();
 185   usgOut->SetPoints(newPts);
 186
 187   usgOut->GetPointData()->CopyAllocate(psIn->GetPointData(),
 188           2*psIn->GetNumberOfPoints());
 189   usgOut->GetCellData()->CopyAllocate(psIn->GetCellData(),
 190                   psIn->GetNumberOfCells());
 191
 192   vtkIdType ncell = psIn->GetNumberOfCells();
 193   vtkSmartPointer<vtkIdList> newIds = vtkSmartPointer<vtkIdList>::New();
 194   vtkSmartPointer<vtkIdList> polyhedronIds = vtkSmartPointer<vtkIdList>::New();
 195   vtkSmartPointer<vtkIdList> neighbors = vtkSmartPointer<vtkIdList>::New();
 196   for(vtkIdType cellId=0; cellId < ncell; cellId++)
 197     {
 198     vtkCell* cell = psIn->GetCell(cellId);
 199     if(cell->GetCellDimension() != 2)
 200       continue;
 201
 202     unsigned char newCellType = VTK_EMPTY_CELL;
 203     unsigned char oldCellType = cell->GetCellType();
 204     switch(oldCellType)
 205       {
 206       case VTK_TRIANGLE :
 207         newCellType = VTK_WEDGE;
 208         break;
 209       case VTK_QUAD :
 210         newCellType = VTK_HEXAHEDRON;
 211         break;
 212       case VTK_POLYGON :
 213         newCellType = VTK_POLYHEDRON;
 214       // default : add new cell types to extrude here
 215       }
 216     if(newCellType == VTK_EMPTY_CELL)
 217       continue;
 218
 219     double cellScalar = array->GetTuple1(cellId);
 220
 221     vtkIdList* oldIds = cell->GetPointIds();
 222     int oldPtsNumber = oldIds->GetNumberOfIds();
 223     int newPtsNumber = oldPtsNumber * 2;
 224     newIds->SetNumberOfIds(newPtsNumber);
 225     double coords[VTK_CELL_SIZE*3];
 226     for(int ptid = 0; ptid < oldPtsNumber; ptid++)
 227       {
 228         psIn->GetPoint(oldIds->GetId(ptid), coords + 3*ptid);
 229       }
 230     for(int ptid = 0; ptid < oldPtsNumber; ptid++)
 231       {
 232       coords[(ptid+oldPtsNumber)*3+0] = coords[ptid*3+0] + cellScalar*dir[0];
 233       coords[(ptid+oldPtsNumber)*3+1] = coords[ptid*3+1] + cellScalar*dir[1];
 234       coords[(ptid+oldPtsNumber)*3+2] = coords[ptid*3+2] + cellScalar*dir[2];
 235       }
 236     double minScalar;
 237     bool minInitialized = false;
 238     for(int ptid = 0; ptid < oldPtsNumber; ptid++)
 239       {
 240       neighbors->Initialize();
 241       psIn->GetPointCells(oldIds->GetId(ptid), neighbors);
 242       for(int neiCellIt = 0; neiCellIt < neighbors->GetNumberOfIds(); neiCellIt++)
 243         {
 244         vtkIdType  neigCellId = neighbors->GetId(neiCellIt);
 245         if(neigCellId == cellId)
 246           continue;
 247         double neighborScalar = array->GetTuple1(neigCellId);
 248         if(neighborScalar != 0.0)
 249                   {
 250           if(!minInitialized)
 251                 minScalar = neighborScalar;
 252           else
 253                     minScalar = (neighborScalar < minScalar ? neighborScalar : minScalar);
 254                   minInitialized = true;
 255                   }
 256         }
 257       if(!minInitialized)
 258         minScalar = 0.0;
 259       }
 260     for(int ptid = 0; ptid < oldPtsNumber; ptid++)
 261       {
 262           if(cellScalar != 0)
 263             {
 264         coords[(ptid)*3+0] = coords[ptid*3+0] + minScalar*dir[0];
 265         coords[(ptid)*3+1] = coords[ptid*3+1] + minScalar*dir[1];
 266         coords[(ptid)*3+2] = coords[ptid*3+2] + minScalar*dir[2];
 267         }
 268       else
 269         {
 270         coords[(ptid+oldPtsNumber)*3+0] = coords[ptid*3+0] + minScalar*dir[0];
 271         coords[(ptid+oldPtsNumber)*3+1] = coords[ptid*3+1] + minScalar*dir[1];
 272         coords[(ptid+oldPtsNumber)*3+2] = coords[ptid*3+2] + minScalar*dir[2];
 273         }
 274       }
 275     for(int ptid=0; ptid<newPtsNumber; ptid++)
 276       {
 277       vtkIdType newId = newPts->InsertNextPoint(coords + 3*ptid);
 278       newIds->SetId(ptid, newId);
 279       usgOut->GetPointData()->CopyData(psIn->GetPointData(),
 280                                        oldIds->GetId(ptid % oldPtsNumber),
 281                                        newIds->GetId(ptid));
 282       }
 283     vtkIdType newCellId;
 284     if(newCellType == VTK_POLYHEDRON)
 285       {
 286       polyhedronIds->Initialize();
 287       // in the polyhedron case, I will generate a quad for each edge
 288       // of the input, and two capping faces
 289       polyhedronIds->InsertNextId(2+oldPtsNumber);
 290       // insert the bottom face
 291       polyhedronIds->InsertNextId(oldPtsNumber);
 292       for(int ptid = 0; ptid < oldPtsNumber; ptid++)
 293         {
 294         polyhedronIds->InsertNextId(newIds->GetId(ptid));
 295         }
 296       // insert the top face
 297       polyhedronIds->InsertNextId(oldPtsNumber);
 298       for(int ptid = oldPtsNumber; ptid < 2*oldPtsNumber; ptid++)
 299         {
 300         polyhedronIds->InsertNextId(newIds->GetId(ptid));
 301         }
 302       // insert the bording quads
 303       for(int qid = 0; qid < oldPtsNumber; qid++)
 304         {
 305         polyhedronIds->InsertNextId(4);
 306         polyhedronIds->InsertNextId(newIds->GetId(qid));
 307         polyhedronIds->InsertNextId(newIds->GetId(qid+oldPtsNumber));
 308         polyhedronIds->InsertNextId(newIds->GetId(qid+((oldPtsNumber+1)%oldPtsNumber)));
 309         polyhedronIds->InsertNextId(newIds->GetId((qid+1)%oldPtsNumber));
 310         }
 311       newIds->Initialize();
 312       for(int jj=0; jj<polyhedronIds->GetNumberOfIds(); jj++)
 313         {
 314         newIds->InsertNextId(polyhedronIds->GetId(jj));
 315         }
 316       }
 317     newCellId = usgOut->InsertNextCell(newCellType, newIds);
 318     usgOut->GetCellData()->CopyData(psIn->GetCellData(),
 319                                    cellId,
 320                                    newCellId);
 321     }
 322
 323   usgOut->GetFieldData()->ShallowCopy(psIn->GetFieldData());
 324
 325   usgOut->Squeeze();
 326
 327   return 1;
 328 }
 329
 330 void  vtkElevationSurfaceFilter::ComputeDirection(vtkPointSet* psIn, double *outDir)
 331 {
 332   double tmp[2][3] = {{0, 0, 0}, {0, 0, 0}};
 333   outDir[0] = outDir[1] = outDir[2] = 0;
 334
 335   vtkPoints* pts = psIn->GetPoints();
 336   vtkSmartPointer<vtkGenericCell> cell = vtkSmartPointer<vtkGenericCell>::New();
 337
 338   for(vtkIdType cellId = 0; cellId < psIn->GetNumberOfCells(); cellId++)
 339     {
 340     psIn->GetCell(cellId, cell);
 341     if(cell->GetCellDimension() != 2)
 342       continue;
 343
 344     vtkIdList* ptIds = cell->GetPointIds();
 345     for(int i=0; i<ptIds->GetNumberOfIds(); i++)
 346       {
 347       vtkIdType firstId = ptIds->GetId(i);
 348       vtkIdType secondId = ptIds->GetId((i+1)%ptIds->GetNumberOfIds());
 349       pts->GetPoint(firstId, tmp[0]);
 350       pts->GetPoint(secondId, tmp[1]);
 351       outDir[0] += tmp[0][1]*tmp[1][2] - tmp[0][2]*tmp[1][1];
 352       outDir[1] += tmp[0][2]*tmp[1][0] - tmp[0][0]*tmp[1][2];
 353       outDir[2] += tmp[0][0]*tmp[1][1] - tmp[0][1]*tmp[1][0];
 354       }
 355     }
 356 }
 357
 358 void vtkElevationSurfaceFilter::PrintSelf(ostream& os, vtkIndent indent)
 359 {
 360   this->Superclass::PrintSelf(os, indent);
 361
 362   os << indent << "ScaleFactor : " << this->ScaleFactor << endl;
 363 }