src/libengrid/vtkImplicitPolyData.cpp

   1 //
   2 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   3 // +                                                                      +
   4 // + This file is part of enGrid.                                         +
   5 // +                                                                      +
   6 // + Copyright 2008-2012 enGits GmbH                                     +
   7 // +                                                                      +
   8 // + enGrid is free software: you can redistribute it and/or modify       +
   9 // + it under the terms of the GNU General Public License as published by +
  10 // + the Free Software Foundation, either version 3 of the License, or    +
  11 // + (at your option) any later version.                                  +
  12 // +                                                                      +
  13 // + enGrid is distributed in the hope that it will be useful,            +
  14 // + but WITHOUT ANY WARRANTY; without even the implied warranty of       +
  15 // + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the        +
  16 // + GNU General Public License for more details.                         +
  17 // +                                                                      +
  18 // + You should have received a copy of the GNU General Public License    +
  19 // + along with enGrid. If not, see <http://www.gnu.org/licenses/>.       +
  20 // +                                                                      +
  21 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  22 //
  23 /*=========================================================================
  24
  25   Program:   Boolean
  26   Module:    $RCSfile: vtkImplicitPolyData.cpp,v $
  27   Language:  C++
  28   Date:      $Date: 2003/07/12 11:22:38 $
  29   Version:   $Revision: 1.1.1.1 $
  30
  31   Copyright (c) 2002 Denis Shamonin
  32   Section Computational Science
  33   University of Amsterdam
  34   Kruislaan 403, 1098 SJ Amsterdam
  35   the Netherlands
  36
  37   E-mail        : dshamoni@science.uva.nl
  38   URL           : http://www.science.uva.nl/~dshamoni/
  39
  40   All rights reserved.
  41
  42      This software is distributed WITHOUT ANY WARRANTY; without even
  43      the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  44      PURPOSE.  See the above copyright notice for more information.
  45
  46 =========================================================================
  47 ToDo:
  48 -   Needs to be restructured in line with vtk principles of on demand
  49     execution. Eg BuildLinks and BuildLocator on first call to Evaluate*?.
  50     Handle input modification correctly etc.
  51 -   Drop internal triangle filter and just check cell type as encountered?.
  52 -   Get working with CellLocator. Currently crashes every time, although
  53     a large number of cells are located successfully. Too hard to debug,
  54     maybe try a tiny dataset.
  55
  56     PLEASE SEND ME YOUR UPDATES, BUG FIXES! at dshamoni@science.uva.nl
  57 =========================================================================*/
  58 //#include <stdafx.h>
  59 #include "vtkImplicitPolyData.h"
  60 #include "vtkPolygon.h"
  61
  62 // Constructor
  63 vtkImplicitPolyData::vtkImplicitPolyData()
  64 {
  65   this->NoGradient[0] = 0.0;
  66   this->NoGradient[1] = 0.0;
  67   this->NoGradient[2] = 1.0;
  68
  69   this->tri = NULL;
  70   this->input = NULL;
  71   this->locator = NULL;
  72   this->poly = vtkPolygon::New();
  73   this->cells = NULL;
  74   this->EvaluateBoundsSet=0;
  75   this->Tolerance = 10.0;
  76 }
  77
  78 void vtkImplicitPolyData::SetInput(vtkPolyData *input)
  79 {
  80   if( this->input != input )
  81     {
  82     vtkDebugMacro( <<" setting Input to " << (void *)input );
  83
  84     // use a tringle filter on the polydata input
  85     // this is done to filter out lines and vertices to leave only
  86     // polygons which are required by this algorithm for cell normals
  87     if( this->tri == NULL )
  88       {
  89       this->tri = vtkTriangleFilter::New();
  90       this->tri->PassVertsOff();
  91       this->tri->PassLinesOff();
  92       }
  93     this->tri->SetInput( input );
  94     this->tri->Update();
  95
  96     this->input = this->tri->GetOutput();
  97     this->input->BuildLinks();  // to enable calls to GetPointCells
  98     this->NoValue = this->input->GetLength();
  99
 100     if( this->locator != NULL ) this->locator->Delete();
 101     locator = PointLocator::New();
 102     this->locator->SetDataSet( this->input );
 103 //  if( this->EvaluateBoundsSet )     // experimental
 104 //  this->locator->SetEvaluateBounds( this->EvaluateBounds );
 105     this->locator->BuildLocator();
 106
 107     if( this->cells   != NULL ) this->cells->Delete();
 108     this->cells = vtkIdList::New();
 109     this->cells->SetNumberOfIds( this->input->GetNumberOfCells() );
 110
 111     this->poly = vtkPolygon::New();
 112     this->Modified();
 113     }
 114 }
 115
 116 // used to extend bounds of point locator
 117 void vtkImplicitPolyData::SetEvaluateBounds( double eBounds[6] )
 118 {
 119   int i;
 120   for( i=0; i<6; i++ )
 121     {
 122     this->EvaluateBounds[i] = eBounds[i];
 123     }
 124
 125   this->EvaluateBoundsSet = 1;
 126 }
 127
 128 unsigned long vtkImplicitPolyData::GetMTime()
 129 {
 130   unsigned long mTime=this->vtkImplicitFunction::GetMTime();
 131   unsigned long inputMTime;
 132
 133   if ( this->input != NULL )
 134     {
 135     this->input->Update ();
 136     inputMTime = this->input->GetMTime();
 137     mTime = ( inputMTime > mTime ? inputMTime : mTime );
 138     }
 139
 140   return mTime;
 141 }
 142
 143 // Destructor
 144 vtkImplicitPolyData::~vtkImplicitPolyData()
 145 {
 146   if( this->tri     != NULL ) this->tri->Delete();
 147   if( this->locator != NULL ) this->locator->Delete();
 148   if( this->poly    != NULL ) this->poly->Delete();
 149   if( this->cells   != NULL ) this->cells->Delete();
 150 }
 151
 152 // Evaluate for point x[3].
 153 // Method using combination of implicit planes
 154 double vtkImplicitPolyData::EvaluateFunction(double x[3])
 155 {
 156   // See if data set with polygons has been specified
 157   if( this->input == NULL || input->GetNumberOfCells() == 0 )
 158     {
 159     vtkErrorMacro(<<"No polygons to evaluate function!");
 160     return this->NoValue;
 161     }
 162
 163   int cellNum, pid;
 164   double dist;
 165   vtkCell *cell;
 166   double dot, ret=-VTK_LARGE_FLOAT, cNormal[3], closestPoint[3];
 167
 168   // get point id of closest point in data set according Tolerance
 169   pid = this->locator->FindClosestPointWithinRadius(Tolerance,x,dist);
 170   if(( pid != -1 ) && (dist < Tolerance))
 171     {
 172     this->input->GetPoint( pid, closestPoint );
 173
 174     // get cells it belongs to
 175     this->input->GetPointCells( pid, cells );
 176     // for each cell
 177     for( cellNum=0; cellNum<cells->GetNumberOfIds(); cellNum++ )
 178       {
 179       cell = this->input->GetCell( cells->GetId( cellNum ) );
 180
 181       // get cell normal
 182       poly->ComputeNormal( cell->GetPoints(), cNormal );
 183
 184       dot = ( cNormal[0]*(x[0]-closestPoint[0]) +
 185               cNormal[1]*(x[1]-closestPoint[1]) +
 186               cNormal[2]*(x[2]-closestPoint[2]) );
 187
 188       if( dot > ret ) ret = dot;
 189       }
 190     }
 191
 192   if( ret == -VTK_LARGE_FLOAT ) ret = NoValue;
 193   return ret;
 194 }
 195
 196 // Evaluate function gradient at point x[3].
 197 void vtkImplicitPolyData::EvaluateGradient( double x[3], double n[3] )
 198 {
 199   int i;
 200   // See if data set with polygons has been specified
 201   if( this->input == NULL || input->GetNumberOfCells() == 0 )
 202     {
 203     vtkErrorMacro(<<"No polygons to evaluate gradient!");
 204     for( i=0; i<3; i++ ) n[i] = this->NoGradient[i];
 205     return;
 206     }
 207
 208   int cellNum, pid;
 209   double dist;
 210   vtkCell *cell;
 211   double dot, ret=-VTK_LARGE_FLOAT, cNormal[3], closestPoint[3];
 212
 213   // get point id of closest point in data set according Tolerance
 214   pid = this->locator->FindClosestPointWithinRadius(Tolerance,x,dist);
 215   if(( pid != -1 ) && (dist < Tolerance))
 216     {
 217     this->input->GetPoint( pid, closestPoint );
 218
 219     // get cells it belongs to
 220     this->input->GetPointCells( pid, cells );
 221     // for each cell
 222     for( cellNum=0; cellNum<cells->GetNumberOfIds(); cellNum++ )
 223        {
 224        cell = this->input->GetCell( cells->GetId( cellNum ) );
 225
 226        // get cell normal
 227        poly->ComputeNormal( cell->GetPoints(), cNormal );
 228
 229        dot = ( cNormal[0]*(x[0]-closestPoint[0]) +
 230                cNormal[1]*(x[1]-closestPoint[1]) +
 231                cNormal[2]*(x[2]-closestPoint[2]) );
 232
 233        if( dot > ret )
 234          {
 235          for( i=0; i<3; i++ ) n[i] = cNormal[i];
 236          }
 237
 238         }
 239     }
 240
 241   if( ret == -VTK_LARGE_FLOAT )
 242     {
 243     for( i=0; i<3; i++ ) n[i] = this->NoGradient[i];
 244     }
 245 }
 246
 247 void vtkImplicitPolyData::PrintSelf(ostream& os, vtkIndent indent)
 248 {
 249   vtkImplicitFunction::PrintSelf(os,indent);
 250
 251   os << indent << "No polydata Value: " << this->NoValue << "\n";
 252   os << indent << "No polydata Gradient: (" << this->NoGradient[0] << ", "
 253      << this->NoGradient[1] << ", " << this->NoGradient[2] << ")\n";
 254
 255   if ( this->input )
 256     {
 257     os << indent << "Input : " << this->input << "\n";
 258     }
 259   else
 260     {
 261     os << indent << "Input : (none)\n";
 262     }
 263
 264   os << indent << "Tolerance: " << this->Tolerance << "\n";
 265 }