Commit 0e06f86c authored by Will Schroeder's avatar Will Schroeder Committed by Kitware Robot
Browse files

Merge topic 'LouisBergmann/vtk-Quadric-Decimation-Vol-constraint'

189a6b8d Better test image thanks Cory
b93218aa Added test
dbd1dc19

 Added Volume Constraint to QuadricDecimation
Acked-by: Kitware Robot's avatarKitware Robot <kwrobot@kitware.com>
Reviewed-by: Cory Quammen's avatarCory Quammen <cory.quammen@kitware.com>
Merge-request: !1954
parents 82c941ce 189a6b8d
Pipeline #26874 canceled with stage
in 0 seconds
......@@ -14,6 +14,7 @@ vtk_add_test_python(
Delaunay3D.py
Delaunay3DAlphaTest.py
QuadricDecimation.py
QuadricDecimation2.py
StreamPolyData.py
TestElevationFilter.py
TestFlyingEdges2D.py
......
#!/usr/bin/env python
import vtk
from vtk.test import Testing
from vtk.util.misc import vtkGetDataRoot
VTK_DATA_ROOT = vtkGetDataRoot()
res = 100
# Test the volume constraint. Create "planar" data with jittery
# z-direction coordinates. Decimate the data with and without
# planar constraints and compare.
#
ps = vtk.vtkPlaneSource()
ps.SetResolution(res,res)
tf = vtk.vtkTriangleFilter()
tf.SetInputConnection(ps.GetOutputPort())
attr = vtk.vtkRandomAttributeGenerator()
attr.SetInputConnection(tf.GetOutputPort())
attr.GeneratePointScalarsOn()
attr.SetMinimumComponentValue(-0.05)
attr.SetMaximumComponentValue(0.05)
# This jitters the geometry
warp = vtk.vtkWarpScalar()
warp.SetInputConnection(attr.GetOutputPort())
warp.SetScaleFactor(0.02)
# Decimator without volume constraint
deci = vtk.vtkQuadricDecimation()
deci.SetInputConnection(warp.GetOutputPort())
deci.SetTargetReduction(.95)
deci.AttributeErrorMetricOn()
deci.VolumePreservationOff()
mapper = vtk.vtkPolyDataMapper()
mapper.SetInputConnection(warp.GetOutputPort())
actor = vtk.vtkActor()
actor.SetMapper(mapper)
# Decimator with volume constraint
deci2 = vtk.vtkQuadricDecimation()
deci2.SetInputConnection(warp.GetOutputPort())
deci2.SetTargetReduction(.95)
deci2.AttributeErrorMetricOn()
deci2.VolumePreservationOn()
mapper2 = vtk.vtkPolyDataMapper()
mapper2.SetInputConnection(deci.GetOutputPort())
actor2 = vtk.vtkActor()
actor2.SetMapper(mapper2)
# Original noisy surface
#
mapper3 = vtk.vtkPolyDataMapper()
mapper3.SetInputConnection(deci2.GetOutputPort())
actor3 = vtk.vtkActor()
actor3.SetMapper(mapper3)
# Create rendering instances
#
ren0 = vtk.vtkRenderer()
ren0.SetViewport(0,0,.33,1)
ren1 = vtk.vtkRenderer()
ren1.SetViewport(0.33,0,0.66,1)
ren2 = vtk.vtkRenderer()
ren2.SetViewport(0.66,0,1,1)
renWin = vtk.vtkRenderWindow()
renWin.AddRenderer(ren0)
renWin.AddRenderer(ren1)
renWin.AddRenderer(ren2)
iren = vtk.vtkRenderWindowInteractor()
iren.SetRenderWindow(renWin)
# Set up the camera parameters
#
camera = vtk.vtkCamera()
camera.SetFocalPoint(0,0,0)
camera.SetPosition(0,0,1)
camera.Elevation(45.0)
ren0.SetActiveCamera(camera)
ren1.SetActiveCamera(camera)
ren2.SetActiveCamera(camera)
# Add the actors to the renderer, set the background and size
#
ren0.AddActor(actor)
ren1.AddActor(actor2)
ren2.AddActor(actor3)
ren0.SetBackground(0,0,0)
ren1.SetBackground(0,0,0)
ren2.SetBackground(0,0,0)
renWin.SetSize(600,300)
ren0.ResetCamera()
renWin.Render()
iren.Initialize()
# gather some information
print ( "Bounds (volume preserve off): ({0})".format( deci.GetOutput().GetPoints().GetBounds() ) )
print ( "Bounds (volume preserve on): ({0})".format( deci2.GetOutput().GetPoints().GetBounds() ) )
#iren.Start()
# --- end of script --
......@@ -66,6 +66,7 @@ vtkQuadricDecimation::vtkQuadricDecimation()
this->EndPoint1List = vtkIdList::New();
this->EndPoint2List = vtkIdList::New();
this->ErrorQuadrics = NULL;
this->VolumeConstraints = NULL;
this->TargetPoints = vtkDoubleArray::New();
this->TargetReduction = 0.9;
......@@ -73,6 +74,7 @@ vtkQuadricDecimation::vtkQuadricDecimation()
this->NumberOfComponents = 0;
this->AttributeErrorMetric = 0;
this->VolumePreservation = 0;
this->ScalarsAttribute = 1;
this->VectorsAttribute = 1;
this->NormalsAttribute = 1;
......@@ -200,7 +202,7 @@ int vtkQuadricDecimation::RequestData(
vtkIdType npts, *pts;
vtkIdType numDeletedTris=0;
// check some assuptiona about the data
// check some assumptions about the data
if (input->GetPolys() == NULL || input->GetPoints() == NULL ||
input->GetPointData() == NULL || input->GetFieldData() == NULL)
{
......@@ -238,6 +240,16 @@ int vtkQuadricDecimation::RequestData(
this->ErrorQuadrics =
new vtkQuadricDecimation::ErrorQuadric[numPts];
if (this->VolumePreservation)
{
this->VolumeConstraints =
new double[numPts * 4];
for (i = 0; i < numPts * 4; i++)
{
this->VolumeConstraints[i] = 0.0;
}
}
vtkDebugMacro(<<"Computing Edges");
this->Edges->InitEdgeInsertion(numPts, 1); // storing edge id as attribute
......@@ -268,19 +280,19 @@ int vtkQuadricDecimation::RequestData(
{
this->ComputeNumberOfComponents();
}
x = new double [3+this->NumberOfComponents];
x = new double [3+this->NumberOfComponents+this->VolumePreservation];
this->CollapseCellIds = vtkIdList::New();
this->TempX = new double [3+this->NumberOfComponents];
this->TempQuad = new double[11 + 4 * this->NumberOfComponents];
this->TempX = new double [3+this->NumberOfComponents+this->VolumePreservation];
this->TempQuad = new double[11 + 4 * this->NumberOfComponents+this->VolumePreservation];
this->TempB = new double [3 + this->NumberOfComponents];
this->TempA = new double*[3 + this->NumberOfComponents];
this->TempData = new double [(3 + this->NumberOfComponents)*(3 + this->NumberOfComponents)];
for (i = 0; i < 3 + this->NumberOfComponents; i++)
this->TempB = new double [3 + this->NumberOfComponents+this->VolumePreservation];
this->TempA = new double*[3 + this->NumberOfComponents+this->VolumePreservation];
this->TempData = new double [(3 + this->NumberOfComponents+this->VolumePreservation)*(3 + this->NumberOfComponents+VolumePreservation)];
for (i = 0; i < 3 + this->NumberOfComponents+this->VolumePreservation; i++)
{
this->TempA[i] = this->TempData+i*(3 + this->NumberOfComponents);
this->TempA[i] = this->TempData+i*(3 + this->NumberOfComponents+this->VolumePreservation);
}
this->TargetPoints->SetNumberOfComponents(3+this->NumberOfComponents);
this->TargetPoints->SetNumberOfComponents(3+this->NumberOfComponents+this->VolumePreservation);
vtkDebugMacro(<<"Computing Quadrics");
this->InitializeQuadrics(numPts);
......@@ -363,6 +375,9 @@ int vtkQuadricDecimation::RequestData(
delete [] this->ErrorQuadrics[i].Quadric;
}
delete [] this->ErrorQuadrics;
if (this->VolumePreservation)
delete[] this->VolumeConstraints;
delete [] x;
this->CollapseCellIds->Delete();
delete [] this->TempX;
......@@ -425,7 +440,7 @@ void vtkQuadricDecimation::InitializeQuadrics(vtkIdType numPts)
A[2] = data+8;
A[3] = data+12;
// allocate local QEM sparce matrix
// allocate local QEM sparse matrix
QEM = new double[11 + 4 * this->NumberOfComponents];
// clear and allocate global QEM array
......@@ -552,15 +567,27 @@ void vtkQuadricDecimation::InitializeQuadrics(vtkIdType numPts)
}
}
// add the QEM to all point of the face
// add the QEM to all points of the face
for (i = 0; i < 3; i++)
{
for (j = 0; j < 11 + 4 * this->NumberOfComponents; j++)
{
this->ErrorQuadrics[pts[i]].Quadric[j] += QEM[j] * triArea2;
}
}
}//for all triangles
// Set volume constraint values g_vol and d_vol
if (this->VolumePreservation)
{
// Vector g_vol
for (j = 0; j < 3; j++)
{
this->VolumeConstraints[pts[i] * 4 + j] += n[j] * triArea2 * 2.0; // triangle normal with length triArea * 2
}
// Scalar d_vol
this->VolumeConstraints[pts[i] * 4 + 3] += -d * triArea2 * 2.0; // (triangle normal with length triArea * 2) * (pts[0] position)
}
}
}//for all triangles
delete [] QEM;
}
......@@ -660,6 +687,14 @@ void vtkQuadricDecimation::AddQuadric(vtkIdType oldPtId, vtkIdType newPtId)
this->ErrorQuadrics[newPtId].Quadric[i] +=
this->ErrorQuadrics[oldPtId].Quadric[i];
}
if (this->VolumePreservation)
{
for (i = 0; i < 4; i++)
{
this->VolumeConstraints[newPtId * 4 + i] += this->VolumeConstraints[oldPtId * 4 + i];
}
}
}
//----------------------------------------------------------------------------
......@@ -890,7 +925,7 @@ double vtkQuadricDecimation::ComputeCost(vtkIdType edgeId, double *x)
double vtkQuadricDecimation::ComputeCost2(vtkIdType edgeId, double *x)
{
// this function is so ugly because the functionality of converting an QEM
// into a dence matrix was not extracted into a separate function and
// into a dense matrix was not extracted into a separate function and
// neither was multiplication and some other matrix and vector primitives
static const double errorNumber = 1e-10;
vtkIdType pointIds[2];
......@@ -908,7 +943,7 @@ double vtkQuadricDecimation::ComputeCost2(vtkIdType edgeId, double *x)
}
// copy the temp quad into TempA
// converting from the sparce matrix format into a dence
// converting from the sparse matrix format into a dense
this->TempA[0][0] = this->TempQuad[0];
this->TempA[0][1] = this->TempA[1][0] = this->TempQuad[1];
this->TempA[0][2] = this->TempA[2][0] = this->TempQuad[2];
......@@ -928,6 +963,8 @@ double vtkQuadricDecimation::ComputeCost2(vtkIdType edgeId, double *x)
this->TempB[i] = -this->TempQuad[11+4*(i-3)+3];
}
// Set zero to all components of the submatrix a[3:n;3:n] and al to its diagonal
for (i = 3; i < 3 + this->NumberOfComponents; i++)
{
for (j = 3; j < 3 + this->NumberOfComponents; j++)
......@@ -942,8 +979,30 @@ double vtkQuadricDecimation::ComputeCost2(vtkIdType edgeId, double *x)
}
}
}
if (this->VolumePreservation)
{
// Add row/col for volume constraint
for (i = 0; i < 3 + this->NumberOfComponents + 1; i++)
{
if (i >= 3)
{
this->TempA[i][3 + this->NumberOfComponents] = 0;
this->TempA[3 + this->NumberOfComponents][i] = 0;
}
else
{
this->TempA[i][3 + this->NumberOfComponents] = this->VolumeConstraints[pointIds[0] * 4 + i];
this->TempA[3 + this->NumberOfComponents][i] = this->VolumeConstraints[pointIds[0] * 4 + i];
this->TempA[i][3 + this->NumberOfComponents] += this->VolumeConstraints[pointIds[1] * 4 + i];
this->TempA[3 + this->NumberOfComponents][i] += this->VolumeConstraints[pointIds[1] * 4 + i];
}
}
// Add constraint to b
this->TempB[3 + this->NumberOfComponents] = this->VolumeConstraints[pointIds[0] * 4 + 3];
this->TempB[3 + this->NumberOfComponents] += this->VolumeConstraints[pointIds[1] * 4 + 3];
}
for (i = 0; i < 3 + this->NumberOfComponents; i++)
for (i = 0; i < 3 + this->NumberOfComponents + this->VolumePreservation; i++)
{
x[i] = this->TempB[i];
}
......@@ -951,7 +1010,7 @@ double vtkQuadricDecimation::ComputeCost2(vtkIdType edgeId, double *x)
// solve A*x = b
// this clobers A
// need to develop a quality of the solution test??
solveOk = vtkMath::SolveLinearSystem(this->TempA, x, 3 + this->NumberOfComponents);
solveOk = vtkMath::SolveLinearSystem(this->TempA, x, 3 + this->NumberOfComponents + this->VolumePreservation);
// need to copy back into A
this->TempA[0][0] = this->TempQuad[0];
......@@ -982,6 +1041,25 @@ double vtkQuadricDecimation::ComputeCost2(vtkIdType edgeId, double *x)
}
}
}
if (this->VolumePreservation)
{
// Add row/col for volume constraint
for (i = 0; i < 3 + this->NumberOfComponents + 1; i++)
{
if (i >= 3)
{
this->TempA[i][3 + this->NumberOfComponents] = 0;
this->TempA[3 + this->NumberOfComponents][i] = 0;
}
else
{
this->TempA[i][3 + this->NumberOfComponents] = this->VolumeConstraints[pointIds[0] * 4 + i];
this->TempA[3 + this->NumberOfComponents][i] = this->VolumeConstraints[pointIds[0] * 4 + i];
this->TempA[i][3 + this->NumberOfComponents] += this->VolumeConstraints[pointIds[1] * 4 + i];
this->TempA[3 + this->NumberOfComponents][i] += this->VolumeConstraints[pointIds[1] * 4 + i];
}
}
}
// check for failure to solve the system
if (!solveOk)
......@@ -1067,15 +1145,15 @@ double vtkQuadricDecimation::ComputeCost2(vtkIdType edgeId, double *x)
// Compute the cost
// x'*A*x - 2*b*x + d
for (i = 0; i < 3+this->NumberOfComponents; i++)
for (i = 0; i < 3+this->NumberOfComponents + this->VolumePreservation; i++)
{
cost += this->TempA[i][i]*x[i]*x[i];
for (j = i+1; j < 3+this->NumberOfComponents; j++)
for (j = i+1; j < 3+this->NumberOfComponents + this->VolumePreservation; j++)
{
cost += 2.0*this->TempA[i][j]*x[i]*x[j];
}
}
for (i = 0; i < 3+this->NumberOfComponents; i++)
for (i = 0; i < 3+this->NumberOfComponents + this->VolumePreservation; i++)
{
cost -= 2.0 * this->TempB[i]*x[i];
}
......@@ -1357,6 +1435,8 @@ void vtkQuadricDecimation::PrintSelf(ostream& os, vtkIndent indent)
os << indent << "Attribute Error Metric: "
<< (this->AttributeErrorMetric ? "On\n" : "Off\n");
os << indent << "Volume Preservation: "
<< (this->VolumePreservation ? "On\n" : "Off\n");
os << indent << "Scalars Attribute: "
<< (this->ScalarsAttribute ? "On\n" : "Off\n");
os << indent << "Vectors Attribute: "
......
......@@ -84,6 +84,16 @@ public:
vtkGetMacro(AttributeErrorMetric, int);
vtkBooleanMacro(AttributeErrorMetric, int);
// Description:
// Decide whether to activate volume preservation which greatly reduces errors
// in triangle normal direction. If off, volume preservation is disabled and
// if AttributeErrorMetric is active, these errors can be large.
// By default VolumePreservation is off
// the attribute errors are off.
vtkSetMacro(VolumePreservation, int);
vtkGetMacro(VolumePreservation, int);
vtkBooleanMacro(VolumePreservation, int);
// Description:
// If attribute errors are to be included in the metric (i.e.,
// AttributeErrorMetric is on), then the following flags control which
......@@ -188,6 +198,7 @@ protected:
double TargetReduction;
double ActualReduction;
int AttributeErrorMetric;
int VolumePreservation;
int ScalarsAttribute;
int VectorsAttribute;
......@@ -215,8 +226,13 @@ protected:
double *Quadric;
};
// One ErrorQuadric per point
ErrorQuadric *ErrorQuadrics;
int AttributeComponents[6];
// Contains 4 doubles per point. Length = nPoints * 4
double *VolumeConstraints;
int AttributeComponents[6];
double AttributeScale[6];
// Temporary variables for performance
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment