bbtk_creaVtk_PKG/src/bbcreaVtkDelaunay3D.cxx

   1 //=====
   2 // Before editing this file, make sure it's a file of your own (i.e.: it wasn't generated from xml description; if so : your modifications will be lost)
   3 //=====
   4 #include "bbcreaVtkDelaunay3D.h"
   5 #include "bbcreaVtkPackage.h"
   6
   7 #include "vtkPoints.h"
   8 #include "vtkDelaunay2D.h"
   9 #include "vtkDelaunay3D.h"
  10 #include "vtkShrinkFilter.h"
  11 #include "vtkGeometryFilter.h"
  12 #include <vtkUnstructuredGrid.h>
  13 #include <vtkCleanPolyData.h>
  14
  15 #include <vtkPCANormalEstimation.h>
  16 #include <vtkSignedDistance.h>
  17 #include <vtkExtractSurface.h>
  18 #include <vtkPointData.h>
  19
  20 #include <vtkSurfaceReconstructionFilter.h>
  21 #include <vtkContourFilter.h>
  22 #include <vtkReverseSense.h>
  23
  24
  25
  26 namespace bbcreaVtk
  27 {
  28
  29 BBTK_ADD_BLACK_BOX_TO_PACKAGE(creaVtk,Delaunay3D)
  30 BBTK_BLACK_BOX_IMPLEMENTATION(Delaunay3D,bbtk::AtomicBlackBox);
  31 //=====
  32 // Before editing this file, make sure it's a file of your own (i.e.: it wasn't generated from xml description; if so : your modifications will be lost)
  33 //=====
  34 void Delaunay3D::Process()
  35 {
  36
  37 // THE MAIN PROCESSING METHOD BODY
  38 //   Here we simply set the input 'In' value to the output 'Out'
  39 //   And print out the output value
  40 // INPUT/OUTPUT ACCESSORS ARE OF THE FORM :
  41 //    void bbSet{Input|Output}NAME(const TYPE&)
  42 //    const TYPE& bbGet{Input|Output}NAME() const
  43 //    Where :
  44 //    * NAME is the name of the input/output
  45 //      (the one provided in the attribute 'name' of the tag 'input')
  46 //    * TYPE is the C++ type of the input/output
  47 //      (the one provided in the attribute 'type' of the tag 'input')
  48 //    bbSetOutputOut( bbGetInputIn() );
  49 //    std::cout << "Output value = " <<bbGetOutputOut() << std::endl;
  50
  51
  52 printf("EED Delaunay3D::Process Start\n");
  53         std::vector<double> lstX=bbGetInputLstPointsX();
  54         std::vector<double> lstY=bbGetInputLstPointsY();
  55         std::vector<double> lstZ=bbGetInputLstPointsZ();
  56         if (lstX.size()!=0)
  57         {
  58                 vtkPoints *points = vtkPoints::New();
  59                 int i,size=lstX.size();
  60                 for (i=0;i<size;i++)
  61                 {
  62                         points->InsertNextPoint( lstX[i], lstY[i], lstZ[i] );
  63                 } // for i
  64
  65                 vtkPolyData *inputpolydata = vtkPolyData::New();
  66                 inputpolydata->SetPoints( points );
  67
  68 /* Delaunay */
  69                 vtkDelaunay3D* delaunay = vtkDelaunay3D::New();
  70                 delaunay->SetInputData( inputpolydata );
  71                 delaunay->SetTolerance( bbGetInputTolerance() ); //0.01
  72                 delaunay->SetAlpha( bbGetInputAlpha() ); //0.2
  73                 delaunay->BoundingTriangulationOff();
  74                 delaunay->Update();
  75
  76                 vtkGeometryFilter *geometry = vtkGeometryFilter::New();
  77                 geometry->SetInputData( delaunay->GetOutput() );
  78                 geometry->Update();
  79
  80                 bbSetOutputOut( geometry->GetOutput() );
  81
  82
  83 /* vtkExtractSurface
  84   double bounds[6];
  85   inputpolydata->GetBounds(bounds);
  86   double range[3];
  87   for (int i = 0; i < 3; ++i)
  88   {
  89     range[i] = bounds[2*i + 1] - bounds[2*i];
  90   }
  91
  92   int sampleSize = inputpolydata->GetNumberOfPoints() * .00005;
  93   if (sampleSize < 10)
  94   {
  95     sampleSize = 10;
  96   }
  97   std::cout << "Sample size is: " << sampleSize << std::endl;
  98   // Do we need to estimate normals?
  99   vtkSmartPointer<vtkSignedDistance> distance =
 100     vtkSmartPointer<vtkSignedDistance>::New();
 101   if (inputpolydata->GetPointData()->GetNormals())
 102   {
 103     std::cout << "Using normals from input file" << std::endl;
 104     distance->SetInputData (inputpolydata);
 105   }  else{
 106         std::cout << "Estimating normals using PCANormalEstimation" << std::endl;
 107     vtkSmartPointer<vtkPCANormalEstimation> normals =
 108       vtkSmartPointer<vtkPCANormalEstimation>::New();
 109     normals->SetInputData (inputpolydata);
 110     normals->SetSampleSize(sampleSize);
 111     normals->SetNormalOrientationToGraphTraversal();
 112     normals->FlipNormalsOn();
 113     distance->SetInputConnection (normals->GetOutputPort());
 114   }
 115   std::cout << "Range: "
 116             << range[0] << ", "
 117             << range[1] << ", "
 118             << range[2] << std::endl;
 119   int dimension = 256;
 120   double radius;
 121
 122   radius = std::max(std::max(range[0], range[1]), range[2])
 123     / static_cast<double>(dimension) * 4; // ~4 voxels
 124
 125 //EED
 126  radius = bbGetInputTolerance();
 127
 128
 129   std::cout << "Radius: " << radius << std::endl;
 130
 131   distance->SetRadius(radius);
 132   distance->SetDimensions(dimension, dimension, dimension);
 133   distance->SetBounds(
 134     bounds[0] - range[0] * .1,
 135     bounds[1] + range[0] * .1,
 136     bounds[2] - range[1] * .1,
 137     bounds[3] + range[1] * .1,
 138     bounds[4] - range[2] * .1,
 139     bounds[5] + range[2] * .1);
 140
 141
 142                 vtkExtractSurface *surface = vtkExtractSurface::New();
 143                 surface->SetInputConnection (distance->GetOutputPort());
 144                 surface->SetRadius(radius * .99 ); //
 145                 surface->Update();
 146                 bbSetOutputOut( surface->GetOutput() );
 147 */
 148
 149 /*vtkPoissonReconstruction
 150                  vtkPolyData *polyData=inputpolydata;
 151
 152                  vtkSmartPointer<vtkPoissonReconstruction> surface =
 153                         vtkSmartPointer<vtkPoissonReconstruction>::New();
 154                   surface->SetDepth(12);
 155
 156                   int sampleSize = polyData->GetNumberOfPoints() * .00005;
 157                   if (sampleSize < 10)
 158                   {
 159                         sampleSize = 10;
 160                   }
 161                   if (polyData->GetPointData()->GetNormals())
 162                   {
 163                         std::cout << "Using normals from input file" << std::endl;
 164                         surface->SetInputData (polyData);
 165                   }
 166                   else
 167                   {
 168                         std::cout << "Estimating normals using PCANormalEstimation" << std::endl;
 169                         vtkSmartPointer<vtkPCANormalEstimation> normals =
 170                           vtkSmartPointer<vtkPCANormalEstimation>::New();
 171                         normals->SetInputData (polyData);
 172                         normals->SetSampleSize(sampleSize);
 173                         normals->SetNormalOrientationToGraphTraversal();
 174                         normals->FlipNormalsOff();
 175                         surface->SetInputConnection(normals->GetOutputPort());
 176                   }
 177
 178                 surface->Update();
 179                 bbSetOutputOut( surface->GetOutput() );
 180
 181 */
 182
 183
 184
 185 /*SurfaceReconstructionFilter
 186
 187                 vtkPolyData *polydata=inputpolydata;
 188
 189                 // Construct the surface and create isosurface.
 190                   vtkSmartPointer<vtkSurfaceReconstructionFilter> surf =
 191                         vtkSmartPointer<vtkSurfaceReconstructionFilter>::New();
 192                   surf->SetInputData(polydata);
 193
 194                   vtkSmartPointer<vtkContourFilter> cf =
 195                         vtkSmartPointer<vtkContourFilter>::New();
 196                   cf->SetInputConnection(surf->GetOutputPort());
 197                   cf->SetValue(0, 0.0);
 198
 199                   // Sometimes the contouring algorithm can create a volume whose gradient
 200                   // vector and ordering of polygon (using the right hand rule) are
 201                   // inconsistent. vtkReverseSense cures this problem.
 202                   vtkSmartPointer<vtkReverseSense> reverse =
 203                         vtkSmartPointer<vtkReverseSense>::New();
 204                   reverse->SetInputConnection(cf->GetOutputPort());
 205                   reverse->ReverseCellsOn();
 206                   reverse->ReverseNormalsOn();
 207
 208                   bbSetOutputOut( reverse->GetOutput() );
 209
 210 */
 211
 212         } else {
 213                 printf("Warnning! Delaunay3D::Process:  list of points empty. \n");
 214         } // if lstX.size
 215
 216 printf("EED Delaunay3D::Process End\n");
 217
 218 }
 219 //=====
 220 // Before editing this file, make sure it's a file of your own (i.e.: it wasn't generated from xml description; if so : your modifications will be lost)
 221 //=====
 222 void Delaunay3D::bbUserSetDefaultValues()
 223 {
 224
 225 //  SET HERE THE DEFAULT INPUT/OUTPUT VALUES OF YOUR BOX
 226 //    Here we initialize the input 'In' to 0
 227    bbSetInputTolerance( 0.01 );
 228    bbSetInputAlpha( 0.2 );
 229    bbSetInputShrinkFactor( 0.9 );
 230
 231 }
 232 //=====
 233 // Before editing this file, make sure it's a file of your own (i.e.: it wasn't generated from xml description; if so : your modifications will be lost)
 234 //=====
 235 void Delaunay3D::bbUserInitializeProcessing()
 236 {
 237
 238 //  THE INITIALIZATION METHOD BODY :
 239 //    Here does nothing
 240 //    but this is where you should allocate the internal/output pointers
 241 //    if any
 242
 243
 244 }
 245 //=====
 246 // Before editing this file, make sure it's a file of your own (i.e.: it wasn't generated from xml description; if so : your modifications will be lost)
 247 //=====
 248 void Delaunay3D::bbUserFinalizeProcessing()
 249 {
 250
 251 //  THE FINALIZATION METHOD BODY :
 252 //    Here does nothing
 253 //    but this is where you should desallocate the internal/output pointers
 254 //    if any
 255
 256 }
 257 }
 258 // EO namespace bbcreaVtk
 259
 260