vtkEdgePoints.cxx 7.03 KB
Newer Older
Will Schroeder's avatar
Will Schroeder committed
1 2
/*=========================================================================

Ken Martin's avatar
Ken Martin committed
3
  Program:   Visualization Toolkit
Ken Martin's avatar
Ken Martin committed
4
  Module:    vtkEdgePoints.cxx
Will Schroeder's avatar
Will Schroeder committed
5

6
  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
7 8
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
9

10 11
     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
12
     PURPOSE.  See the above copyright notice for more information.
Will Schroeder's avatar
Will Schroeder committed
13 14

=========================================================================*/
Ken Martin's avatar
Ken Martin committed
15
#include "vtkEdgePoints.h"
16

17 18 19
#include "vtkCell.h"
#include "vtkCellArray.h"
#include "vtkCellData.h"
20
#include "vtkDataSet.h"
21
#include "vtkGenericCell.h"
22 23
#include "vtkInformation.h"
#include "vtkInformationVector.h"
24
#include "vtkMergePoints.h"
25
#include "vtkObjectFactory.h"
26
#include "vtkPointData.h"
27
#include "vtkPolyData.h"
28

Brad King's avatar
Brad King committed
29
vtkStandardNewMacro(vtkEdgePoints);
Will Schroeder's avatar
Will Schroeder committed
30 31

// Construct object with contour value of 0.0.
Ken Martin's avatar
Ken Martin committed
32
vtkEdgePoints::vtkEdgePoints()
Will Schroeder's avatar
Will Schroeder committed
33 34
{
  this->Value = 0.0;
35 36 37 38 39 40
  this->Locator = vtkMergePoints::New();
}

vtkEdgePoints::~vtkEdgePoints()
{
  this->Locator->Delete();
41
  this->Locator = nullptr;
Will Schroeder's avatar
Will Schroeder committed
42 43 44 45 46
}

//
// General filter: handles arbitrary input.
//
47 48 49 50
int vtkEdgePoints::RequestData(
  vtkInformation *vtkNotUsed(request),
  vtkInformationVector **inputVector,
  vtkInformationVector *outputVector)
Will Schroeder's avatar
Will Schroeder committed
51
{
52 53 54 55
  // get the info objects
  vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
  vtkInformation *outInfo = outputVector->GetInformationObject(0);

56
  // get the input and output
57 58 59 60 61
  vtkDataSet *input = vtkDataSet::SafeDownCast(
    inInfo->Get(vtkDataObject::DATA_OBJECT()));
  vtkPolyData *output = vtkPolyData::SafeDownCast(
    outInfo->Get(vtkDataObject::DATA_OBJECT()));

62
  vtkDataArray *inScalars;
63
  vtkPoints *newPts;
Ken Martin's avatar
Ken Martin committed
64
  vtkCellArray *newVerts;
Amy Squillacote's avatar
Amy Squillacote committed
65 66
  vtkIdType cellId, ptId, edgeId, newCellId;
  int above, below, i, numEdges;
67
  vtkCell *edge;
Ken Martin's avatar
Ken Martin committed
68
  double range[2];
69 70
  double s0, s1, x0[3], x1[3], x[3], t;
  double e0Scalar, deltaScalar;
71
  int e0, e1;
Amy Squillacote's avatar
Amy Squillacote committed
72
  vtkIdType p1, p2;
73
  vtkIdType pts[1];
Amy Squillacote's avatar
Amy Squillacote committed
74
  vtkIdType numCells, estimatedSize;
75
  vtkDataArray *cellScalars;
76 77
  vtkPointData *inPd=input->GetPointData(), *outPd=output->GetPointData();
  vtkCellData *inCd=input->GetCellData(), *outCd=output->GetCellData();
Will Schroeder's avatar
Will Schroeder committed
78

Ken Martin's avatar
Ken Martin committed
79
  vtkDebugMacro(<< "Generating edge points");
80

Ken Martin's avatar
Ken Martin committed
81 82
  // Initialize and check input
  //
83
  if ( ! (inScalars = input->GetPointData()->GetScalars()) )
84
  {
Ken Martin's avatar
Ken Martin committed
85
    vtkErrorMacro(<<"No scalar data to contour");
86
    return 1;
87
  }
Will Schroeder's avatar
Will Schroeder committed
88

89
  inScalars->GetRange(range,0);
90
  if ( this->Value < range[0] || this->Value > range[1] )
91
  {
Ken Martin's avatar
Ken Martin committed
92
    vtkWarningMacro(<<"Value lies outside of scalar range");
93
    return 1;
94
  }
95

96
  numCells = input->GetNumberOfCells();
97
  estimatedSize = static_cast<vtkIdType>(numCells * .75);
98
  estimatedSize = estimatedSize / 1024 * 1024; //multiple of 1024
Bill Lorensen's avatar
Bill Lorensen committed
99
  if (estimatedSize < 1024)
100
  {
Bill Lorensen's avatar
Bill Lorensen committed
101
    estimatedSize = 1024;
102
  }
103

104
  newPts = vtkPoints::New();
105
  newPts->Allocate(estimatedSize, estimatedSize/2);
Will Schroeder's avatar
Will Schroeder committed
106
  newVerts = vtkCellArray::New();
107
  newVerts->Allocate(estimatedSize, estimatedSize/2);
108 109
  cellScalars = inScalars->NewInstance();
  cellScalars->SetNumberOfComponents(inScalars->GetNumberOfComponents());
110
  cellScalars->Allocate(VTK_CELL_SIZE*inScalars->GetNumberOfComponents());
111

112
  this->Locator->InitPointInsertion (newPts, input->GetBounds());
113

114
  // interpolate data along edge; copy cell data
115
  outPd->InterpolateAllocate(inPd,5000,10000);
116
  outCd->CopyAllocate(inCd,5000,10000);
117

118 119
  // Traverse all edges. Since edges are not explicitly represented, use a
  // trick: traverse all cells and obtain cell edges and then cell edge
luz.paz's avatar
luz.paz committed
120
  // neighbors. If cell id < all edge neighbors ids, then this edge has not
121 122
  // yet been visited and is processed.
  //
123
  int abort=0;
Amy Squillacote's avatar
Amy Squillacote committed
124
  vtkIdType progressInterval = numCells/20 + 1;
125 126
  vtkGenericCell *cell = vtkGenericCell::New();
  for (cellId=0; cellId < numCells && !abort; cellId++)
127
  {
128
    if ( ! (cellId % progressInterval) )
129
    {
130
      vtkDebugMacro(<<"Processing #" << cellId);
131
      this->UpdateProgress (static_cast<double>(cellId)/numCells);
132
      abort = this->GetAbortExecute();
133
    }
134 135

    input->GetCell(cellId,cell);
136
    inScalars->GetTuples(cell->PointIds, cellScalars);
Will Schroeder's avatar
Will Schroeder committed
137

138
    // loop over cell points to check if cell straddles isosurface value
Will Schroeder's avatar
Will Schroeder committed
139
    for ( above=below=0, ptId=0; ptId < cell->GetNumberOfPoints(); ptId++ )
140
    {
141
      if ( cellScalars->GetComponent(ptId,0) >= this->Value )
142
      {
Will Schroeder's avatar
Will Schroeder committed
143
        above = 1;
144
      }
145
      else
146
      {
Will Schroeder's avatar
Will Schroeder committed
147 148
        below = 1;
      }
149
    }
Will Schroeder's avatar
Will Schroeder committed
150 151

    if ( above && below ) //contour passes through cell
152
    {
Will Schroeder's avatar
Will Schroeder committed
153
      if ( cell->GetCellDimension() < 2 ) //only points can be generated
154
      {
155
        cell->Contour(this->Value, cellScalars, this->Locator, newVerts,
156
                      nullptr, nullptr, inPd, outPd, inCd, cellId, outCd);
157
      }
Will Schroeder's avatar
Will Schroeder committed
158 159

      else //
160
      {
Will Schroeder's avatar
Will Schroeder committed
161 162
        numEdges = cell->GetNumberOfEdges();
        for (edgeId=0; edgeId < numEdges; edgeId++)
163
        {
Will Schroeder's avatar
Will Schroeder committed
164
          edge = cell->GetEdge(edgeId);
165
          inScalars->GetTuples(edge->PointIds, cellScalars);
Will Schroeder's avatar
Will Schroeder committed
166

167 168
          s0 = cellScalars->GetComponent(0,0);
          s1 = cellScalars->GetComponent(1,0);
169
          if ( (s0 < this->Value && s1 >= this->Value) ||
Will Schroeder's avatar
Will Schroeder committed
170
          (s0 >= this->Value && s1 < this->Value) )
171
          {
172
            //ordering intersection direction avoids numerical problems
173
            deltaScalar = s1 - s0;
174
            if (deltaScalar > 0)
175
            {
176 177
              e0 = 0; e1 = 1;
              e0Scalar = s0;
178
            }
179
            else
180
            {
181 182 183
              e0 = 1; e1 = 0;
              e0Scalar = s1;
              deltaScalar = -deltaScalar;
184
            }
185 186

            t = (this->Value - e0Scalar) / deltaScalar;
187

188 189
            edge->Points->GetPoint(e0,x0);
            edge->Points->GetPoint(e1,x1);
190

Bill Lorensen's avatar
Bill Lorensen committed
191
            for (i=0; i<3; i++)
192
            {
193
              x[i] = x0[i] + t * (x1[i] - x0[i]);
194
            }
195
            if ( this->Locator->InsertUniquePoint(x, pts[0]) )
196
            {
197
              newCellId = newVerts->InsertNextCell(1,pts);
198
              outCd->CopyData(inCd,cellId,newCellId);
199 200
              p1 = edge->PointIds->GetId(e0);
              p2 = edge->PointIds->GetId(e1);
201
              outPd->InterpolateEdge(inPd,pts[0],p1,p2,t);
202 203 204 205 206 207
            } //if point not created before
          } //if edge straddles contour value
        } //for each edge
      } //dimension 2 and higher
    } //above and below
  } //for all cells
208
  cell->Delete();
Will Schroeder's avatar
Will Schroeder committed
209

Ken Martin's avatar
Ken Martin committed
210
  vtkDebugMacro(<<"Created: " << newPts->GetNumberOfPoints() << " points");
211

212 213
  // Update ourselves.  Because we don't know up front how many verts we've
  // created, take care to reclaim memory.
214
  //
Ken Martin's avatar
Ken Martin committed
215
  output->SetPoints(newPts);
216 217
  newPts->Delete();

Ken Martin's avatar
Ken Martin committed
218
  output->SetVerts(newVerts);
219 220
  newVerts->Delete();

221
  this->Locator->Initialize();//free up any extra memory
Ken Martin's avatar
Ken Martin committed
222
  output->Squeeze();
223

224
  cellScalars->Delete();
225 226 227 228

  return 1;
}

229 230 231 232 233 234
int vtkEdgePoints::FillInputPortInformation(int, vtkInformation *info)
{
  info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkDataSet");
  return 1;
}

235
void vtkEdgePoints::PrintSelf(ostream& os, vtkIndent indent)
Will Schroeder's avatar
Will Schroeder committed
236
{
Brad King's avatar
Brad King committed
237
  this->Superclass::PrintSelf(os,indent);
Will Schroeder's avatar
Will Schroeder committed
238

239
  os << indent << "Contour Value: " << this->Value << "\n";
Will Schroeder's avatar
Will Schroeder committed
240
}