You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
1264 lines
36 KiB
1264 lines
36 KiB
/*=========================================================================
|
|
|
|
Program: Visualization Toolkit
|
|
Module: $RCSfile: vtkGreedyTerrainDecimation.cxx,v $
|
|
|
|
Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
|
|
All rights reserved.
|
|
See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
|
|
|
|
This software is distributed WITHOUT ANY WARRANTY; without even
|
|
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
|
|
PURPOSE. See the above copyright notice for more information.
|
|
|
|
=========================================================================*/
|
|
#include "vtkGreedyTerrainDecimation.h"
|
|
#include "vtkInformation.h"
|
|
#include "vtkInformationVector.h"
|
|
#include "vtkObjectFactory.h"
|
|
#include "vtkPriorityQueue.h"
|
|
#include "vtkImageData.h"
|
|
#include "vtkPolyData.h"
|
|
#include "vtkPointData.h"
|
|
#include "vtkPoints.h"
|
|
#include "vtkCellArray.h"
|
|
#include "vtkTriangle.h"
|
|
#include "vtkDoubleArray.h"
|
|
#include "vtkFloatArray.h"
|
|
#include "vtkMath.h"
|
|
|
|
#include <vtkstd/vector>
|
|
|
|
vtkCxxRevisionMacro(vtkGreedyTerrainDecimation, "$Revision: 1.22.12.1 $");
|
|
vtkStandardNewMacro(vtkGreedyTerrainDecimation);
|
|
|
|
// Define some constants describing vertices
|
|
//
|
|
#define VTK_VERTEX_NO_TRIANGLE -1
|
|
#define VTK_VERTEX_INSERTED -2
|
|
|
|
#define VTK_IN_TRIANGLE 0
|
|
#define VTK_INTERIOR_EDGE 1
|
|
#define VTK_BOUNDARY_EDGE 2
|
|
|
|
//Supporting classes for points
|
|
class vtkTerrainInfo
|
|
{
|
|
public:
|
|
vtkTerrainInfo():TriangleId(VTK_VERTEX_NO_TRIANGLE) {}
|
|
vtkIdType TriangleId;
|
|
};
|
|
|
|
//PIMPL STL encapsulation
|
|
//
|
|
// Maps input point ids to owning mesh triangle
|
|
class vtkGreedyTerrainDecimationTerrainInfoType : public vtkstd::vector<vtkTerrainInfo>
|
|
{
|
|
public:
|
|
typedef vtkstd::vector<vtkTerrainInfo> Superclass;
|
|
typedef Superclass::size_type size_type;
|
|
vtkGreedyTerrainDecimationTerrainInfoType(size_type n, const vtkTerrainInfo& value):
|
|
vtkstd::vector<vtkTerrainInfo>(n,value) {}
|
|
};
|
|
|
|
// Maps mesh point id to input point id
|
|
class vtkGreedyTerrainDecimationPointInfoType : public vtkstd::vector<vtkIdType> {};
|
|
|
|
|
|
// Begin vtkGreedyTerrainDecimation class implementation-----------------------------------------
|
|
//
|
|
vtkGreedyTerrainDecimation::vtkGreedyTerrainDecimation()
|
|
{
|
|
this->ErrorMeasure = VTK_ERROR_SPECIFIED_REDUCTION;
|
|
this->NumberOfTriangles = 1000;
|
|
this->Reduction = 0.90;
|
|
this->AbsoluteError = 1;
|
|
this->RelativeError = 0.01;
|
|
this->BoundaryVertexDeletion = 1;
|
|
|
|
this->ComputeNormals = 0;
|
|
this->Normals = 0;
|
|
}
|
|
|
|
vtkGreedyTerrainDecimation::~vtkGreedyTerrainDecimation()
|
|
{
|
|
}
|
|
|
|
inline void vtkGreedyTerrainDecimation::GetTerrainPoint(int i, int j, double x[3])
|
|
{
|
|
x[0] = this->Origin[0] + i*this->Spacing[0];
|
|
x[1] = this->Origin[1] + j*this->Spacing[1];
|
|
}
|
|
|
|
inline void vtkGreedyTerrainDecimation::ComputeImageCoordinates(vtkIdType inputPtId, int ij[2])
|
|
{
|
|
ij[0] = inputPtId % this->Dimensions[0];
|
|
ij[1] = inputPtId / this->Dimensions[0];
|
|
}
|
|
|
|
inline vtkIdType vtkGreedyTerrainDecimation::InsertNextPoint(vtkIdType inputPtId,
|
|
double x[3])
|
|
{
|
|
if ( (this->CurrentPointId+1) >= (vtkIdType)this->PointInfo->size() )
|
|
{
|
|
this->PointInfo->resize(2*this->PointInfo->size());
|
|
}
|
|
|
|
double *ptr = this->Points->WritePointer(3*this->CurrentPointId,3);
|
|
*ptr++ = *x++;
|
|
*ptr++ = *x++;
|
|
*ptr = *x;
|
|
|
|
this->OutputPD->CopyData(this->InputPD,inputPtId,this->CurrentPointId);
|
|
(*this->PointInfo)[this->CurrentPointId] = inputPtId;
|
|
|
|
return this->CurrentPointId++;
|
|
}
|
|
|
|
inline double *vtkGreedyTerrainDecimation::GetPoint(vtkIdType id)
|
|
{
|
|
return this->Points->GetPointer(3*id);
|
|
}
|
|
|
|
inline void vtkGreedyTerrainDecimation::GetPoint(vtkIdType id, double x[3])
|
|
{
|
|
double *ptr = this->Points->GetPointer(3*id);
|
|
x[0] = *ptr++;
|
|
x[1] = *ptr++;
|
|
x[2] = *ptr;
|
|
}
|
|
|
|
void vtkGreedyTerrainDecimation::EstimateOutputSize(const vtkIdType numInputPts,
|
|
vtkIdType &numPts, vtkIdType &numTris)
|
|
{
|
|
switch (this->ErrorMeasure)
|
|
{
|
|
case VTK_ERROR_NUMBER_OF_TRIANGLES:
|
|
numTris = this->NumberOfTriangles;
|
|
break;
|
|
case VTK_ERROR_SPECIFIED_REDUCTION:
|
|
numTris = static_cast<int>(2*numInputPts*(1.0-this->Reduction));
|
|
break;
|
|
default:
|
|
numTris = numInputPts;
|
|
}
|
|
|
|
numPts = numTris/2 + 1;
|
|
numPts = (numPts < 4 ? 4 : numPts); //insure enough storage for initial four corner points
|
|
return;
|
|
}
|
|
|
|
int vtkGreedyTerrainDecimation::SatisfiesErrorMeasure(double error)
|
|
{
|
|
switch (this->ErrorMeasure)
|
|
{
|
|
case VTK_ERROR_NUMBER_OF_TRIANGLES:
|
|
if ( this->Mesh->GetNumberOfPolys() >= this->NumberOfTriangles ) return 1;
|
|
break;
|
|
|
|
case VTK_ERROR_SPECIFIED_REDUCTION:
|
|
{
|
|
double reduction = (double)this->Mesh->GetNumberOfPolys()/this->MaximumNumberOfTriangles;
|
|
if ( (1.0 - reduction) <= this->Reduction ) return 1;
|
|
}
|
|
break;
|
|
|
|
case VTK_ERROR_ABSOLUTE:
|
|
if ( error <= this->AbsoluteError ) return 1;
|
|
break;
|
|
|
|
case VTK_ERROR_RELATIVE:
|
|
if ( (error/this->Length) <= this->RelativeError ) return 1;
|
|
break;
|
|
}
|
|
|
|
return 0;
|
|
}
|
|
|
|
//Update all triangles connected to this mesh point
|
|
void vtkGreedyTerrainDecimation::UpdateTriangles(vtkIdType ptId)
|
|
{
|
|
unsigned short ncells;
|
|
vtkIdType *cells, npts, *pts;
|
|
|
|
this->Mesh->GetPointCells(ptId,ncells,cells);
|
|
for (unsigned short i=0; i<ncells; i++)
|
|
{
|
|
this->Mesh->GetCellPoints(cells[i], npts, pts);
|
|
this->UpdateTriangle(cells[i], (*this->PointInfo)[pts[0]],
|
|
(*this->PointInfo)[pts[1]], (*this->PointInfo)[pts[2]]);
|
|
}
|
|
}
|
|
|
|
//Update all points as to which triangle they lie in. Basically a scanline algorithm.
|
|
void vtkGreedyTerrainDecimation::UpdateTriangle(vtkIdType triId,
|
|
vtkIdType p1, vtkIdType p2, vtkIdType p3)
|
|
{
|
|
// Scan convert triangle / update points as to which triangle contains each point
|
|
int ij1[2], ij2[2], ij3[2];
|
|
this->ComputeImageCoordinates(p1, ij1);
|
|
this->ComputeImageCoordinates(p2, ij2);
|
|
this->ComputeImageCoordinates(p3, ij3);
|
|
|
|
double h[4]; //extra entry added for interpolated value
|
|
h[0] = (double) this->Heights->GetTuple1(p1);
|
|
h[1] = (double) this->Heights->GetTuple1(p2);
|
|
h[2] = (double) this->Heights->GetTuple1(p3);
|
|
|
|
this->UpdateTriangle(triId, ij1, ij2, ij3, h);
|
|
}
|
|
|
|
void vtkGreedyTerrainDecimation::InsertBoundaryVertices()
|
|
{
|
|
int i, j;
|
|
vtkIdType inputPtId, offset;
|
|
|
|
// Insert vertices around boundary of image
|
|
// Note that the four corner vertices are already inserted.
|
|
// Along x-axis at y=0.
|
|
for (i=1; i<(this->Dimensions[0]-1); i++)
|
|
{
|
|
inputPtId = i;
|
|
this->AddPointToTriangulation(inputPtId);
|
|
}
|
|
|
|
// Along x-axis at y=dim[1].
|
|
offset = this->Dimensions[0]*(this->Dimensions[1]-1);
|
|
for (i=1; i<(this->Dimensions[0]-1); i++)
|
|
{
|
|
inputPtId = offset + i;
|
|
this->AddPointToTriangulation(inputPtId);
|
|
}
|
|
|
|
// Along y-axis at x=0.
|
|
for (j=1; j<(this->Dimensions[1]-1); j++)
|
|
{
|
|
inputPtId = j*this->Dimensions[0];
|
|
this->AddPointToTriangulation(inputPtId);
|
|
}
|
|
|
|
// Along y-axis at x=dims[0].
|
|
offset = this->Dimensions[0]-1;
|
|
for (j=1; j<(this->Dimensions[1]-1); j++)
|
|
{
|
|
inputPtId = offset + j*this->Dimensions[0];
|
|
this->AddPointToTriangulation(inputPtId);
|
|
}
|
|
}
|
|
|
|
// Determine whether point x is inside of circumcircle of triangle
|
|
// defined by points (x1, x2, x3). Returns non-zero if inside circle.
|
|
// (Note that z-component is ignored.)
|
|
int vtkGreedyTerrainDecimation::InCircle (double x[3], double x1[3], double x2[3],
|
|
double x3[3])
|
|
{
|
|
double radius2, center[2], dist2;
|
|
|
|
radius2 = vtkTriangle::Circumcircle(x1,x2,x3,center);
|
|
|
|
// check if inside/outside circumcircle
|
|
dist2 = (x[0]-center[0]) * (x[0]-center[0]) +
|
|
(x[1]-center[1]) * (x[1]-center[1]);
|
|
|
|
if ( dist2 < (0.999999999999*radius2) )
|
|
{
|
|
return 1;
|
|
}
|
|
else
|
|
{
|
|
return 0;
|
|
}
|
|
}
|
|
|
|
#define VTK_DEL2D_TOLERANCE 1.0e-014
|
|
|
|
// Recursive method to locate triangle containing point. Starts with arbitrary
|
|
// triangle (tri) and "walks" towards it. Influenced by some of Guibas and
|
|
// Stolfi's work. Returns id of enclosing triangle, or -1 if no triangle
|
|
// found. Also, the array nei[3] is used to communicate info about points
|
|
// that lie on triangle edges: nei[0] is neighboring triangle id, and nei[1]
|
|
// and nei[2] are the vertices defining the edge.
|
|
vtkIdType vtkGreedyTerrainDecimation::FindTriangle(double x[3], vtkIdType ptIds[3],
|
|
vtkIdType tri, double tol,
|
|
vtkIdType nei[3], vtkIdList *neighbors,
|
|
int& status)
|
|
{
|
|
int i, j, ir, ic, inside, i2, i3;
|
|
vtkIdType *pts, npts, newNei;
|
|
double p[3][3], n[2], vp[2], vx[2], dp, minProj;
|
|
|
|
// get local triangle info
|
|
this->Mesh->GetCellPoints(tri,npts,pts);
|
|
for (i=0; i<3; i++)
|
|
{
|
|
ptIds[i] = pts[i];
|
|
this->GetPoint(ptIds[i], p[i]);
|
|
}
|
|
|
|
// Randomization (of find edge neighbors) avoids walking in
|
|
// circles in certain weird cases
|
|
srand(tri);
|
|
ir = rand() % 3;
|
|
|
|
// evaluate in/out of each edge
|
|
for (inside=1, minProj=VTK_DEL2D_TOLERANCE, ic=0; ic<3; ic++)
|
|
{
|
|
i = (ir+ic) % 3;
|
|
i2 = (i+1) % 3;
|
|
i3 = (i+2) % 3;
|
|
|
|
// create a 2D edge normal to define a "half-space"; evaluate points (i.e.,
|
|
// candiate point and other triangle vertex not on this edge).
|
|
n[0] = -(p[i2][1] - p[i][1]);
|
|
n[1] = p[i2][0] - p[i][0];
|
|
vtkMath::Normalize2D(n);
|
|
|
|
// compute local vectors
|
|
for (j=0; j<2; j++)
|
|
{
|
|
vp[j] = p[i3][j] - p[i][j];
|
|
vx[j] = x[j] - p[i][j];
|
|
}
|
|
|
|
//check for duplicate point
|
|
vtkMath::Normalize2D(vp);
|
|
if ( vtkMath::Normalize2D(vx) <= tol )
|
|
{
|
|
vtkErrorMacro("Duplicate point");
|
|
return -1;
|
|
}
|
|
|
|
// see if two points are in opposite half spaces
|
|
dp = vtkMath::Dot2D(n,vx) * (vtkMath::Dot2D(n,vp) < 0 ? -1.0 : 1.0);
|
|
if ( dp < VTK_DEL2D_TOLERANCE )
|
|
{
|
|
if ( dp < minProj ) //track edge most orthogonal to point direction
|
|
{
|
|
inside = 0;
|
|
nei[1] = ptIds[i];
|
|
nei[2] = ptIds[i2];
|
|
minProj = dp;
|
|
}
|
|
}//outside this edge
|
|
}//for each edge
|
|
|
|
if ( inside ) // all edges have tested positive
|
|
{
|
|
nei[0] = (-1);
|
|
status = VTK_IN_TRIANGLE;
|
|
return tri;
|
|
}
|
|
|
|
else if ( !inside && (fabs(minProj) < VTK_DEL2D_TOLERANCE) ) // on edge
|
|
{
|
|
this->Mesh->GetCellEdgeNeighbors(tri,nei[1],nei[2],neighbors);
|
|
if ( neighbors->GetNumberOfIds() < 1 )
|
|
{
|
|
nei[0] = (-1);
|
|
status = VTK_BOUNDARY_EDGE;
|
|
}
|
|
else
|
|
{
|
|
nei[0] = neighbors->GetId(0);
|
|
status = VTK_INTERIOR_EDGE;
|
|
}
|
|
return tri;
|
|
}
|
|
|
|
else //walk towards point
|
|
{
|
|
this->Mesh->GetCellEdgeNeighbors(tri,nei[1],nei[2],neighbors);
|
|
if ( (newNei=neighbors->GetId(0)) == nei[0] )
|
|
{
|
|
vtkErrorMacro("Degeneracy");
|
|
return -1;
|
|
}
|
|
else
|
|
{
|
|
nei[0] = tri;
|
|
return this->FindTriangle(x,ptIds,newNei,tol,nei,neighbors,status);
|
|
}
|
|
}
|
|
}
|
|
|
|
#undef VTK_DEL2D_TOLERANCE
|
|
|
|
// Recursive method checks whether edge is Delaunay, and if not, swaps edge.
|
|
// Continues until all edges are Delaunay. Points p1 and p2 form the edge in
|
|
// question; x is the coordinates of the inserted point; tri is the current
|
|
// triangle id.
|
|
void vtkGreedyTerrainDecimation::CheckEdge(vtkIdType ptId, double x[3], vtkIdType p1,
|
|
vtkIdType p2, vtkIdType tri, int depth)
|
|
{
|
|
if ( depth > 15 )
|
|
{
|
|
return;
|
|
}
|
|
|
|
int i;
|
|
vtkIdType *pts, npts, numNei, nei, p3;
|
|
double x1[3], x2[3], x3[3];
|
|
vtkIdList *neighbors;
|
|
vtkIdType swapTri[3];
|
|
|
|
this->GetPoint(p1,x1);
|
|
this->GetPoint(p2,x2);
|
|
|
|
neighbors = vtkIdList::New();
|
|
neighbors->Allocate(2);
|
|
|
|
this->Mesh->GetCellEdgeNeighbors(tri,p1,p2,neighbors);
|
|
numNei = neighbors->GetNumberOfIds();
|
|
|
|
if ( numNei > 0 ) //i.e., not a boundary edge
|
|
{
|
|
// get neighbor info including opposite point
|
|
nei = neighbors->GetId(0);
|
|
this->Mesh->GetCellPoints(nei, npts, pts);
|
|
for (i=0; i<2; i++)
|
|
{
|
|
if ( pts[i] != p1 && pts[i] != p2 )
|
|
{
|
|
break;
|
|
}
|
|
}
|
|
p3 = pts[i];
|
|
this->GetPoint(p3,x3);
|
|
|
|
// see whether point is in circumcircle
|
|
if ( this->InCircle (x3, x, x1, x2) )
|
|
{// swap diagonal
|
|
this->Mesh->RemoveReferenceToCell(p1,tri);
|
|
this->Mesh->RemoveReferenceToCell(p2,nei);
|
|
this->Mesh->ResizeCellList(ptId,1);
|
|
this->Mesh->AddReferenceToCell(ptId,nei);
|
|
this->Mesh->ResizeCellList(p3,1);
|
|
this->Mesh->AddReferenceToCell(p3,tri);
|
|
|
|
swapTri[0] = ptId; swapTri[1] = p3; swapTri[2] = p2;
|
|
this->Mesh->ReplaceCell(tri,3,swapTri);
|
|
|
|
swapTri[0] = ptId; swapTri[1] = p1; swapTri[2] = p3;
|
|
this->Mesh->ReplaceCell(nei,3,swapTri);
|
|
|
|
// two new edges become suspect
|
|
this->CheckEdge(ptId, x, p3, p2, tri, ++depth);
|
|
this->CheckEdge(ptId, x, p1, p3, nei, depth);
|
|
|
|
}//in circle
|
|
}//interior edge
|
|
|
|
neighbors->Delete();
|
|
}
|
|
|
|
vtkIdType vtkGreedyTerrainDecimation::AddPointToTriangulation(vtkIdType inputPtId)
|
|
{
|
|
vtkIdType ptId, nei[3], tri[4];
|
|
vtkIdType nodes[4][3], pts[3], numNeiPts, *neiPts;
|
|
vtkIdType i, p1=0, p2=0;
|
|
int ij[2];
|
|
double x[3];
|
|
int status;
|
|
|
|
//Make sure the point has not been previously inserted
|
|
if ( (*this->TerrainInfo)[inputPtId].TriangleId == VTK_VERTEX_INSERTED )
|
|
{
|
|
return -1;
|
|
}
|
|
|
|
//Start off by determining the image coordinates and the position
|
|
this->ComputeImageCoordinates(inputPtId, ij);
|
|
this->GetTerrainPoint(ij[0], ij[1], x);
|
|
x[2] = (double) this->Heights->GetTuple1(inputPtId);
|
|
|
|
//Seed the search
|
|
nei[0] = (*this->TerrainInfo)[inputPtId].TriangleId;
|
|
tri[0] = (nei[0] < 0 ? 0 : nei[0]);
|
|
|
|
tri[0] = this->FindTriangle(x,pts,tri[0],this->Tolerance,nei,this->Neighbors,status);
|
|
if ( tri[0] >= 0 ) //found a triangle
|
|
{
|
|
// Insert the point into the output
|
|
ptId = this->InsertNextPoint(inputPtId, x);
|
|
|
|
if (this->Normals)
|
|
{
|
|
float n[3];
|
|
this->ComputePointNormal(ij[0], ij[1], n);
|
|
this->Normals->InsertNextTuple(n);
|
|
}
|
|
if ( status == VTK_IN_TRIANGLE ) //in triangle
|
|
{
|
|
//delete this triangle; create three new triangles
|
|
//first triangle is replaced with one of the new ones
|
|
nodes[0][0] = ptId; nodes[0][1] = pts[0]; nodes[0][2] = pts[1];
|
|
this->Mesh->RemoveReferenceToCell(pts[2], tri[0]);
|
|
this->Mesh->ReplaceCell(tri[0], 3, nodes[0]);
|
|
this->Mesh->InsertNextLinkedPoint(3);
|
|
this->Mesh->AddReferenceToCell(ptId,tri[0]);
|
|
|
|
//create two new triangles
|
|
nodes[1][0] = ptId; nodes[1][1] = pts[1]; nodes[1][2] = pts[2];
|
|
tri[1] = this->Mesh->InsertNextLinkedCell(VTK_TRIANGLE, 3, nodes[1]);
|
|
|
|
nodes[2][0] = ptId; nodes[2][1] = pts[2]; nodes[2][2] = pts[0];
|
|
tri[2] = this->Mesh->InsertNextLinkedCell(VTK_TRIANGLE, 3, nodes[2]);
|
|
|
|
// Check edge neighbors for Delaunay criterion. If not satisfied, flip
|
|
// edge diagonal. (This is done recursively.)
|
|
this->CheckEdge(ptId, x, pts[0], pts[1], tri[0], 0);
|
|
this->CheckEdge(ptId, x, pts[1], pts[2], tri[1], 0);
|
|
this->CheckEdge(ptId, x, pts[2], pts[0], tri[2], 0);
|
|
}
|
|
|
|
else if ( status == VTK_INTERIOR_EDGE ) // on interior triangle edge; has a neighbor
|
|
{
|
|
//update cell list
|
|
this->Mesh->GetCellPoints(nei[0],numNeiPts,neiPts);
|
|
for (i=0; i<3; i++)
|
|
{
|
|
if ( neiPts[i] != nei[1] && neiPts[i] != nei[2] )
|
|
{
|
|
p1 = neiPts[i];
|
|
}
|
|
if ( pts[i] != nei[1] && pts[i] != nei[2] )
|
|
{
|
|
p2 = pts[i];
|
|
}
|
|
}
|
|
this->Mesh->ResizeCellList(p1,1);
|
|
this->Mesh->ResizeCellList(p2,1);
|
|
|
|
//replace two triangles
|
|
this->Mesh->RemoveReferenceToCell(nei[2],tri[0]);
|
|
this->Mesh->RemoveReferenceToCell(nei[2],nei[0]);
|
|
|
|
nodes[0][0] = ptId; nodes[0][1] = p2; nodes[0][2] = nei[1];
|
|
this->Mesh->ReplaceCell(tri[0], 3, nodes[0]);
|
|
|
|
nodes[1][0] = ptId; nodes[1][1] = nei[1]; nodes[1][2] = p1;
|
|
this->Mesh->ReplaceCell(nei[0], 3, nodes[1]);
|
|
|
|
this->Mesh->InsertNextLinkedPoint(4);
|
|
this->Mesh->AddReferenceToCell(ptId,tri[0]);
|
|
this->Mesh->AddReferenceToCell(ptId,nei[0]);
|
|
|
|
tri[1] = nei[0];
|
|
|
|
//create two new triangles
|
|
nodes[2][0] = ptId; nodes[2][1] = nei[2]; nodes[2][2] = p2;
|
|
tri[2] = this->Mesh->InsertNextLinkedCell(VTK_TRIANGLE, 3, nodes[2]);
|
|
|
|
nodes[3][0] = ptId; nodes[3][1] = p1; nodes[3][2] = nei[2];
|
|
tri[3] = this->Mesh->InsertNextLinkedCell(VTK_TRIANGLE, 3, nodes[3]);
|
|
|
|
// Check edge neighbors for Delaunay criterion.
|
|
for ( i=0; i<4; i++ )
|
|
{
|
|
this->CheckEdge (ptId, x, nodes[i][1], nodes[i][2], tri[i], 0);
|
|
}
|
|
}
|
|
|
|
else //if ( status == VTK_BOUNDARY_EDGE ) // on boundary triangle edge; no neighbor
|
|
{
|
|
//update cell list
|
|
for (i=0; i<3; i++)
|
|
{
|
|
if ( pts[i] != nei[1] && pts[i] != nei[2] )
|
|
{
|
|
p1 = pts[i];
|
|
}
|
|
}
|
|
this->Mesh->ResizeCellList(p1,1);
|
|
|
|
//replace one triangle
|
|
this->Mesh->RemoveReferenceToCell(nei[2],tri[0]);
|
|
|
|
nodes[0][0] = ptId; nodes[0][1] = p1; nodes[0][2] = nei[1];
|
|
this->Mesh->ReplaceCell(tri[0], 3, nodes[0]);
|
|
|
|
this->Mesh->InsertNextLinkedPoint(2);
|
|
this->Mesh->AddReferenceToCell(ptId,tri[0]);
|
|
|
|
//create one new triangles
|
|
nodes[1][0] = ptId; nodes[1][1] = nei[2]; nodes[1][2] = p1;
|
|
tri[1] = this->Mesh->InsertNextLinkedCell(VTK_TRIANGLE, 3, nodes[1]);
|
|
|
|
// Check edge neighbors for Delaunay criterion.
|
|
for ( i=0; i<2; i++ )
|
|
{
|
|
this->CheckEdge (ptId, x, nodes[i][1], nodes[i][2], tri[i], 0);
|
|
}
|
|
}
|
|
|
|
//Indicate that it is now inserted and reinsert errors
|
|
(*this->TerrainInfo)[inputPtId].TriangleId = VTK_VERTEX_INSERTED;
|
|
this->UpdateTriangles(ptId);
|
|
|
|
}//if triangle containing point found
|
|
|
|
return 0;
|
|
}
|
|
|
|
void vtkGreedyTerrainDecimation::ComputePointNormal(int i, int j, float n[3])
|
|
{
|
|
vtkDataArray* scalars;
|
|
double x0, x1, y0, y1, dx, dy;
|
|
float vx[3], vy[3];
|
|
|
|
scalars = this->InputPD->GetScalars();
|
|
|
|
dx = dy = 0;
|
|
// X
|
|
if (i > 0)
|
|
{
|
|
x0 = scalars->GetTuple1(j*this->Dimensions[0] + i - 1);
|
|
dx += this->Spacing[0];
|
|
}
|
|
else
|
|
{
|
|
x0 = scalars->GetTuple1(j*this->Dimensions[0] + i);
|
|
}
|
|
if (i < this->Dimensions[0]-1)
|
|
{
|
|
x1 = scalars->GetTuple1(j*this->Dimensions[0] + i + 1);
|
|
dx += this->Spacing[0];
|
|
}
|
|
else
|
|
{
|
|
x1 = scalars->GetTuple1(j*this->Dimensions[0] + i);
|
|
}
|
|
// Y
|
|
if (j > 0)
|
|
{
|
|
y0 = scalars->GetTuple1((j-1)*this->Dimensions[0] + i);
|
|
dy += this->Spacing[1];
|
|
}
|
|
else
|
|
{
|
|
y0 = scalars->GetTuple1(j*this->Dimensions[0] + i);
|
|
}
|
|
if (j < this->Dimensions[1]-1)
|
|
{
|
|
y1 = scalars->GetTuple1((j+1)*this->Dimensions[0] + i);
|
|
dy += this->Spacing[1];
|
|
}
|
|
else
|
|
{
|
|
y1 = scalars->GetTuple1(j*this->Dimensions[0] + i);
|
|
}
|
|
|
|
if (dx == 0.0 || dy == 0.0)
|
|
{
|
|
// This would only happen if the input was not an XY image.
|
|
vtkErrorMacro("Could not compute normal.");
|
|
return;
|
|
}
|
|
vx[0] = (float)(dx);
|
|
vx[1] = 0.0;
|
|
vx[2] = (float)(x1-x0);
|
|
vy[0] = 0.0;
|
|
vy[1] = (float)(dy);
|
|
vy[2] = (float)(y1-y0);
|
|
vtkMath::Cross(vx, vy, n);
|
|
vtkMath::Normalize(n);
|
|
}
|
|
|
|
int vtkGreedyTerrainDecimation::RequestData(
|
|
vtkInformation *vtkNotUsed(request),
|
|
vtkInformationVector **inputVector,
|
|
vtkInformationVector *outputVector)
|
|
{
|
|
// get the info objects
|
|
vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
|
|
vtkInformation *outInfo = outputVector->GetInformationObject(0);
|
|
|
|
// get the input and ouptut
|
|
vtkImageData *input = vtkImageData::SafeDownCast(
|
|
inInfo->Get(vtkDataObject::DATA_OBJECT()));
|
|
vtkPolyData *output = vtkPolyData::SafeDownCast(
|
|
outInfo->Get(vtkDataObject::DATA_OBJECT()));
|
|
|
|
vtkIdType numInputPts=input->GetNumberOfPoints(), numPts, numTris;
|
|
vtkIdType inputPtId;
|
|
double error, bounds[6], center[3];
|
|
vtkCellArray *triangles;
|
|
this->Mesh = output;
|
|
this->InputPD = input->GetPointData();
|
|
this->OutputPD = this->Mesh->GetPointData();
|
|
|
|
// Check input and initialize the algorithm.
|
|
//
|
|
vtkDebugMacro(<<"Decimating terrain...");
|
|
|
|
if ( input->GetDataDimension() != 2 )
|
|
{
|
|
vtkWarningMacro(<<"This class treats 2D height fields only");
|
|
return 1;
|
|
}
|
|
if ( (this->Heights = this->InputPD->GetScalars()) == NULL )
|
|
{
|
|
vtkWarningMacro(<<"This class requires height scalars");
|
|
return 1;
|
|
}
|
|
|
|
input->GetBounds(bounds);
|
|
input->GetCenter(center);
|
|
input->GetDimensions(this->Dimensions);
|
|
double *origin = input->GetOrigin();
|
|
double *spacing = input->GetSpacing();
|
|
for (int ii=0; ii<3; ii++)
|
|
{
|
|
this->Origin[ii] = (double)origin[ii];
|
|
this->Spacing[ii] = (double)spacing[ii];
|
|
}
|
|
this->Length = input->GetLength();
|
|
this->MaximumNumberOfTriangles = 2 * (this->Dimensions[0]-1) * (this->Dimensions[1]-1);
|
|
this->NumberOfTriangles = (this->NumberOfTriangles < this->MaximumNumberOfTriangles ?
|
|
this->NumberOfTriangles : this->MaximumNumberOfTriangles);
|
|
|
|
// Points within this tolerance are considered coincident...should not happen
|
|
this->Tolerance = 0.01 * this->Spacing[0];
|
|
|
|
// Scratch data structures
|
|
this->Neighbors = vtkIdList::New(); this->Neighbors->Allocate(2);
|
|
|
|
// Top element of VTK's priority queue returns the minimum error value. Since we want the
|
|
// maximum error, we use 1/error relationship to insert errors.
|
|
this->TerrainError = vtkPriorityQueue::New();
|
|
this->TerrainError->Allocate(numInputPts, (vtkIdType)((double)0.25*numInputPts));
|
|
|
|
// Initialize the triangle mesh data structures. Double precision point coordinates
|
|
// are required because of the numerical requirements on the Delaunay algorithm.
|
|
//
|
|
this->EstimateOutputSize(numInputPts, numPts, numTris);
|
|
|
|
vtkPoints *newPts = vtkPoints::New();
|
|
newPts->SetDataTypeToDouble();
|
|
newPts->Allocate(numPts);
|
|
this->Points = static_cast<vtkDoubleArray *>(newPts->GetData());
|
|
|
|
// initailize the normals
|
|
if (this->ComputeNormals)
|
|
{
|
|
this->Normals = vtkFloatArray::New();
|
|
this->Normals->SetNumberOfComponents(3);
|
|
this->Normals->Allocate(numPts*3);
|
|
this->Normals->SetName("Normals");
|
|
}
|
|
|
|
// Supplemental arrays used to accelerate the algorithm.
|
|
// TerrainInfo contains the "containing" triangle for each point. PointInfo maps the
|
|
// triangle mesh point id to the input image point id.
|
|
this->TerrainInfo = new vtkGreedyTerrainDecimationTerrainInfoType(numInputPts,vtkTerrainInfo());
|
|
|
|
this->PointInfo = new vtkGreedyTerrainDecimationPointInfoType;
|
|
this->PointInfo->resize(numPts);
|
|
|
|
// Setup the point attributes
|
|
this->OutputPD->CopyAllocate(this->InputPD,numPts);
|
|
|
|
// Begin the algorithm proper. The image is initially triangulated with two triangles whose
|
|
// four vertices are located at the corners of the input image.
|
|
//
|
|
newPts->Allocate(numPts);
|
|
|
|
inputPtId = 0;
|
|
newPts->InsertPoint(0, bounds[0],bounds[2],
|
|
(double)this->Heights->GetTuple1(inputPtId)); //ptId=0
|
|
this->OutputPD->CopyData(this->InputPD,inputPtId,0);
|
|
(*this->PointInfo)[0] = inputPtId;
|
|
|
|
inputPtId = this->Dimensions[0] - 1;
|
|
newPts->InsertPoint(1, bounds[1],bounds[2],
|
|
(double)this->Heights->GetTuple1(inputPtId)); //ptId=1
|
|
this->OutputPD->CopyData(this->InputPD,inputPtId,1);
|
|
(*this->PointInfo)[1] = inputPtId;
|
|
|
|
inputPtId = this->Dimensions[0]*this->Dimensions[1] - 1;
|
|
newPts->InsertPoint(2, bounds[1],bounds[3],
|
|
(double)this->Heights->GetTuple1(inputPtId)); //ptId=2
|
|
this->OutputPD->CopyData(this->InputPD,inputPtId,2);
|
|
(*this->PointInfo)[2] = inputPtId;
|
|
|
|
inputPtId = this->Dimensions[0]*(this->Dimensions[1]-1);
|
|
newPts->InsertPoint(3, bounds[0],bounds[3],
|
|
(double)this->Heights->GetTuple1(inputPtId)); //ptId=3
|
|
this->OutputPD->CopyData(this->InputPD,inputPtId,3);
|
|
(*this->PointInfo)[3] = inputPtId;
|
|
this->CurrentPointId = 4;
|
|
|
|
// Handle normals of the four corners.
|
|
if (this->Normals)
|
|
{
|
|
float n[3];
|
|
this->ComputePointNormal(0, 0, n);
|
|
this->Normals->InsertNextTuple(n);
|
|
this->ComputePointNormal(this->Dimensions[0]-1, 0, n);
|
|
this->Normals->InsertNextTuple(n);
|
|
this->ComputePointNormal(this->Dimensions[0]-1, this->Dimensions[1]-1, n);
|
|
this->Normals->InsertNextTuple(n);
|
|
this->ComputePointNormal(0, this->Dimensions[1]-1, n);
|
|
this->Normals->InsertNextTuple(n);
|
|
}
|
|
|
|
// Insert initial triangles into output mesh
|
|
triangles = vtkCellArray::New();
|
|
triangles->Allocate(numTris,3);
|
|
|
|
triangles->InsertNextCell(3);
|
|
triangles->InsertCellPoint(0); triangles->InsertCellPoint(1); triangles->InsertCellPoint(3);
|
|
|
|
triangles->InsertNextCell(3);
|
|
triangles->InsertCellPoint(1); triangles->InsertCellPoint(2); triangles->InsertCellPoint(3);
|
|
|
|
// Construct the topological hierarchy for the output mesh. The alternative BuildLinks(num)
|
|
// call reallocates the links from the points to the using triangles.
|
|
this->Mesh->SetPoints(newPts);
|
|
this->Mesh->SetPolys(triangles);
|
|
this->Mesh->BuildLinks(numPts); //build cell structure; give it initial size
|
|
|
|
// Update all (two) triangles connected to this mesh point. The single point
|
|
// in each triangle with maximum error are inserted into the error queue.
|
|
//
|
|
this->UpdateTriangles(3);
|
|
|
|
// If boundary vertex deletion is not allowed, insert the boundary
|
|
// points first.
|
|
if ( ! this->BoundaryVertexDeletion )
|
|
{
|
|
this->InsertBoundaryVertices();
|
|
}
|
|
|
|
// While the error metric is not satisfied, add point with greatest error.
|
|
// Note that this algorithm can terminate "prematurely" (e.g. compared to
|
|
// the number of triangles) if the maximum error in the queue becomes zero.
|
|
//
|
|
int abortExecute=0;
|
|
vtkIdType numInsertedPoints=0;
|
|
int tenth=numPts/10+1;
|
|
|
|
while ( !abortExecute && (inputPtId = this->TerrainError->Pop(0, error)) >= 0 )
|
|
{
|
|
if ( this->SatisfiesErrorMeasure((1.0/error)) )
|
|
{
|
|
break;
|
|
}
|
|
else
|
|
{
|
|
this->AddPointToTriangulation(inputPtId);
|
|
if ( ! (++numInsertedPoints % tenth) )
|
|
{
|
|
this->UpdateProgress( (double)(numInsertedPoints>numPts?numPts:numInsertedPoints)/numPts);
|
|
abortExecute = this->GetAbortExecute();
|
|
}
|
|
}
|
|
}
|
|
|
|
if (this->Normals)
|
|
{
|
|
this->OutputPD->SetNormals(this->Normals);
|
|
this->Normals->Delete();
|
|
this->Normals = 0;
|
|
}
|
|
|
|
vtkDebugMacro(<<"Output TIN contains: " << this->Mesh->GetNumberOfPoints() << " points"
|
|
<<"and " << this->Mesh->GetNumberOfPolys() << " triangles");
|
|
|
|
// The output triangle data was created incrementally by the Delaunay algorithm.
|
|
// Here we just clean up the data structures.
|
|
//
|
|
this->Neighbors->Delete();
|
|
this->TerrainError->Delete();
|
|
delete this->TerrainInfo;
|
|
delete this->PointInfo;
|
|
|
|
newPts->Delete();
|
|
triangles->Delete();
|
|
|
|
return 1;
|
|
}
|
|
|
|
/*----------------------------------------------------------------------
|
|
"Scan conversion" routines to update all points lying in a triangle.
|
|
|
|
Divide a triangle into two subtriangles as shown.
|
|
|
|
o max
|
|
/ \
|
|
| \
|
|
/ \
|
|
| \
|
|
midL o..........o midR
|
|
| _/
|
|
/ _/
|
|
| _/
|
|
/ _/
|
|
| _/
|
|
/_/
|
|
o min
|
|
|
|
This way we can scan the two subtriangles independently without worrying about
|
|
the transistion in interpolation that occurs at the vertices.
|
|
|
|
A triangle may be characterized in one of four ways:
|
|
VTK_TWO_TRIANGLES: We can create a two triangle representation
|
|
VTK_BOTTOM_TRIANGLE: We should only scan the lower triangle
|
|
VTK_TOP_TRIANGLE: We should only scan the upper triangle
|
|
VTK_DEGENERATE: The points are colinear (not scan converted)
|
|
|
|
Configuration of the two triangles
|
|
--------------------------------------------------------------------------*/
|
|
#define VTK_TWO_TRIANGLES 0 //most often
|
|
#define VTK_BOTTOM_TRIANGLE 1
|
|
#define VTK_TOP_TRIANGLE 2
|
|
#define VTK_DEGENERATE 3 //should never happen in this application
|
|
|
|
//---------------------------------------------------------------------------
|
|
// Update all points lying in the given triangle. This means indicating the triangle
|
|
// that the point is in, plus computing the error in the height field.
|
|
//
|
|
void vtkGreedyTerrainDecimation::UpdateTriangle(vtkIdType tri, int ij1[2], int ij2[2], int ij3[2],
|
|
double h[3])
|
|
{
|
|
int *min, *max, *midL, *midR, *mid, mid2[2];
|
|
double t, tt;
|
|
int i, j, xL, xR;
|
|
double hMin, hMax, hMidL, hMidR, hL, hR;
|
|
vtkIdType idx, inputPtId, maxInputPtId=0;
|
|
double error, maxError=0.0;
|
|
|
|
int type = this->CharacterizeTriangle(ij1, ij2, ij3, min, max, midL, midR, mid, mid2,
|
|
h, hMin, hMax, hMidL, hMidR);
|
|
|
|
switch(type)
|
|
{
|
|
case VTK_BOTTOM_TRIANGLE:
|
|
case VTK_TWO_TRIANGLES:
|
|
for (j=min[1]+1; j<midL[1]; j++) //for all scan lines; skip vertices
|
|
{
|
|
idx = j*this->Dimensions[0];
|
|
t = (double)(j - min[1]) / (midL[1] - min[1]);
|
|
xL = (int)((1.0-t)*min[0] + t*midL[0]);
|
|
xR = (int)((1.0-t)*min[0] + t*midR[0]);
|
|
hL = (1.0-t)*hMin + t*hMidL;
|
|
hR = (1.0-t)*hMin + t*hMidR;
|
|
for (i=xL; i<=xR; i++)
|
|
{
|
|
inputPtId = i + idx;
|
|
if ( (*this->TerrainInfo)[inputPtId].TriangleId != VTK_VERTEX_INSERTED )
|
|
{
|
|
(*this->TerrainInfo)[inputPtId].TriangleId = tri;
|
|
if ( (xR-xL) > 0 )
|
|
{
|
|
tt = (double)(i-xL) / (xR-xL);
|
|
error = (1.0-tt)*hL + tt*hR;
|
|
}
|
|
else
|
|
{
|
|
error = hL;
|
|
}
|
|
error = fabs( (double)this->Heights->GetTuple1(inputPtId) - error );
|
|
if ( error > maxError )
|
|
{
|
|
maxError = error;
|
|
maxInputPtId = inputPtId;
|
|
}
|
|
} //if vertex not inserted
|
|
} //for this scanline
|
|
} //for all scanlines in this triangle
|
|
|
|
if ( type == VTK_BOTTOM_TRIANGLE )
|
|
{
|
|
break;
|
|
}
|
|
|
|
case VTK_TOP_TRIANGLE:
|
|
//Start scanning the upper triangle
|
|
for (j=max[1]-1; j>midL[1]; j--) //for all scan lines; skip vertices
|
|
{
|
|
idx = j*this->Dimensions[0];
|
|
t = (double)(j - midL[1]) / (max[1] - midL[1]);
|
|
xL = (int) (t*max[0] + (1.0-t)*midL[0]);
|
|
xR = (int) (t*max[0] + (1.0-t)*midR[0]);
|
|
hL = t*hMax + (1.0-t)*hMidL;
|
|
hR = t*hMax + (1.0-t)*hMidR;
|
|
for (i=xL; i<=xR; i++)
|
|
{
|
|
inputPtId = i + idx;
|
|
if ( (*this->TerrainInfo)[inputPtId].TriangleId != VTK_VERTEX_INSERTED )
|
|
{
|
|
(*this->TerrainInfo)[inputPtId].TriangleId = tri;
|
|
if ( (xR-xL) > 0 )
|
|
{
|
|
tt = (double)(i-xL) / (xR-xL);
|
|
error = (1.0-tt)*hL + tt*hR;
|
|
}
|
|
else
|
|
{
|
|
error = hL;
|
|
}
|
|
error = fabs( (double)this->Heights->GetTuple1(inputPtId) - error );
|
|
if ( error > maxError )
|
|
{
|
|
maxError = error;
|
|
maxInputPtId = inputPtId;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
break;
|
|
|
|
default:
|
|
return;
|
|
}
|
|
|
|
//The maximum error in the triangle has been found. Insert it into the queue.
|
|
if ( maxError > 0.0 )
|
|
{
|
|
this->TerrainError->DeleteId(maxInputPtId); //if previously inserted
|
|
this->TerrainError->Insert((1.0/maxError),maxInputPtId);
|
|
}
|
|
}
|
|
|
|
|
|
// Characterize the configuration of the triangle based on image coordinates
|
|
// (All points in triangulation are from an image).
|
|
//
|
|
int vtkGreedyTerrainDecimation::CharacterizeTriangle(int ij1[2], int ij2[2], int ij3[3],
|
|
int* &min, int* &max, int* &midL, int* &midR,
|
|
int* &mid, int mid2[2], double h[3],
|
|
double &hMin, double &hMax, double &hL,
|
|
double &hR)
|
|
{
|
|
// Check for situations where one edge of triangle is horizontal
|
|
//
|
|
if ( ij1[1] == ij2[1] )
|
|
{
|
|
if ( ij1[0] < ij2[0] )
|
|
{
|
|
midL = ij1;
|
|
midR = ij2;
|
|
hL = h[0];
|
|
hR = h[1];
|
|
}
|
|
else
|
|
{
|
|
midL = ij2;
|
|
midR = ij1;
|
|
hL = h[1];
|
|
hR = h[0];
|
|
}
|
|
if( ij3[1] < ij1[1])
|
|
{
|
|
min = ij3;
|
|
hMin = h[2];
|
|
return VTK_BOTTOM_TRIANGLE;
|
|
}
|
|
else
|
|
{
|
|
max = ij3;
|
|
hMax = h[2];
|
|
return VTK_TOP_TRIANGLE;
|
|
}
|
|
}
|
|
|
|
else if ( ij2[1] == ij3[1] )
|
|
{
|
|
if ( ij2[0] < ij3[0] )
|
|
{
|
|
midL = ij2;
|
|
midR = ij3;
|
|
hL = h[1];
|
|
hR = h[2];
|
|
}
|
|
else
|
|
{
|
|
midL = ij3;
|
|
midR = ij2;
|
|
hL = h[2];
|
|
hR = h[1];
|
|
}
|
|
if( ij1[1] < ij2[1])
|
|
{
|
|
min = ij1;
|
|
hMin = h[0];
|
|
return VTK_BOTTOM_TRIANGLE;
|
|
}
|
|
else
|
|
{
|
|
max = ij1;
|
|
hMax = h[0];
|
|
return VTK_TOP_TRIANGLE;
|
|
}
|
|
}
|
|
|
|
else if ( ij3[1] == ij1[1] )
|
|
{
|
|
if ( ij3[0] < ij1[0] )
|
|
{
|
|
midL = ij3;
|
|
midR = ij1;
|
|
hL = h[2];
|
|
hR = h[0];
|
|
}
|
|
else
|
|
{
|
|
midL = ij1;
|
|
midR = ij3;
|
|
hL = h[0];
|
|
hR = h[2];
|
|
}
|
|
if( ij2[1] < ij3[1])
|
|
{
|
|
min = ij2;
|
|
hMin = h[1];
|
|
return VTK_BOTTOM_TRIANGLE;
|
|
}
|
|
else
|
|
{
|
|
max = ij2;
|
|
hMax = h[1];
|
|
return VTK_TOP_TRIANGLE;
|
|
}
|
|
}
|
|
|
|
// Default situation (two triangles with no horizontal edges).
|
|
// Determine max, min and mid vertices.
|
|
//
|
|
// Find minimum
|
|
if ( ij1[1] < ij2[1] )
|
|
{
|
|
if ( ij1[1] < ij3[1] )
|
|
{
|
|
min = ij1;
|
|
hMin = h[0];
|
|
}
|
|
else
|
|
{
|
|
min = ij3;
|
|
hMin = h[2];
|
|
}
|
|
}
|
|
else
|
|
{
|
|
if ( ij2[1] < ij3[1] )
|
|
{
|
|
min = ij2;
|
|
hMin = h[1];
|
|
}
|
|
else
|
|
{
|
|
min = ij3;
|
|
hMin = h[2];
|
|
}
|
|
}
|
|
|
|
// Find maximum
|
|
if ( ij1[1] > ij2[1] )
|
|
{
|
|
if ( ij1[1] > ij3[1] )
|
|
{
|
|
max = ij1;
|
|
hMax = h[0];
|
|
}
|
|
else
|
|
{
|
|
max = ij3;
|
|
hMax = h[2];
|
|
}
|
|
}
|
|
else
|
|
{
|
|
if ( ij2[1] > ij3[1] )
|
|
{
|
|
max = ij2;
|
|
hMax = h[1];
|
|
}
|
|
else
|
|
{
|
|
max = ij3;
|
|
hMax = h[2];
|
|
}
|
|
}
|
|
|
|
// Find the midL and midR
|
|
double hMid, hMid2;
|
|
if ( ij1 != min && ij1 != max)
|
|
{
|
|
mid = ij1;
|
|
hMid = h[0];
|
|
}
|
|
else if ( ij2 != min && ij2 != max)
|
|
{
|
|
mid = ij2;
|
|
hMid = h[1];
|
|
}
|
|
else //if ( ij3 != min && ij2 != max)
|
|
{
|
|
mid = ij3;
|
|
hMid = h[2];
|
|
}
|
|
|
|
// Computation of the intersection
|
|
//
|
|
mid2[1] = mid[1];
|
|
double t = (double) (mid2[1] - min[1]) / (max[1] - min[1]);
|
|
mid2[0] = (int) ((1.0-t)*min[0] + t*max[0] + 0.5); //rounding
|
|
hMid2 = (1.0-t)*hMin + t*hMax;
|
|
|
|
if ( mid[0] < mid2[0] )
|
|
{
|
|
midL = mid;
|
|
midR = mid2;
|
|
hL = hMid;
|
|
hR = hMid2;
|
|
}
|
|
else
|
|
{
|
|
midL = mid2;
|
|
midR = mid;
|
|
hL = hMid2;
|
|
hR = hMid;
|
|
}
|
|
|
|
return VTK_TWO_TRIANGLES;
|
|
}
|
|
|
|
int vtkGreedyTerrainDecimation::FillInputPortInformation(int, vtkInformation *info)
|
|
{
|
|
info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkImageData");
|
|
return 1;
|
|
}
|
|
|
|
void vtkGreedyTerrainDecimation::PrintSelf(ostream& os, vtkIndent indent)
|
|
{
|
|
this->Superclass::PrintSelf(os,indent);
|
|
|
|
os << indent << "Error Measure: ";
|
|
if ( this->ErrorMeasure == VTK_ERROR_NUMBER_OF_TRIANGLES )
|
|
{
|
|
os << "Number of triangles\n";
|
|
os << indent << "Number of triangles: " << this->NumberOfTriangles << "\n";
|
|
}
|
|
else if ( this->ErrorMeasure == VTK_ERROR_SPECIFIED_REDUCTION )
|
|
{
|
|
os << "Specified reduction\n";
|
|
os << indent << "Reduction: " << this->Reduction << "\n";
|
|
}
|
|
else if ( this->ErrorMeasure == VTK_ERROR_ABSOLUTE )
|
|
{
|
|
os << "Absolute\n";
|
|
os << indent << "Absolute Error: " << this->AbsoluteError << "\n";
|
|
}
|
|
else // this->ErrorMeasure == VTK_ERROR_RELATIVE
|
|
{
|
|
os << "Relative\n";
|
|
os << indent << "Relative Error: " << this->RelativeError << "\n";
|
|
}
|
|
|
|
os << indent << "BoundaryVertexDeletion: "
|
|
<< (this->BoundaryVertexDeletion ? "On\n" : "Off\n");
|
|
os << indent << "ComputeNormals: "
|
|
<< (this->ComputeNormals ? "On\n" : "Off\n");
|
|
}
|
|
|