Commit 67122e84 authored by Cory Quammen's avatar Cory Quammen

BUG: Fixed error in area-weighted normal

When the closest point on a surface to a test point is on an edge, the
area-weighted normal is the average of the normals of the faces
adjacent to the edge. The contribution from one of the faces was being
added twice, producing the incorrect area-weighted normal. This
resulted in the wrong sign for the signed distance in some cases.
This patch fixes the error.

Changed TestImplicitPolyDataDistance.cxx to fail prior to this change
and pass after the change. The previous test exercised
vtkImplicityPolyDataDistance but did not check the results. The new
test checks against a baseline.
parent 050c1726
......@@ -18,7 +18,7 @@ vtk_add_test_cxx(${vtk-module}CxxTests tests
TestFeatureEdges.cxx,NO_VALID
TestGlyph3D.cxx
TestHedgeHog.cxx,NO_VALID
TestImplicitPolyDataDistance.cxx,NO_VALID
TestImplicitPolyDataDistance.cxx
TestMaskPoints.cxx,NO_VALID
TestNamedComponents.cxx,NO_VALID
TestPolyDataConnectivityFilter.cxx,NO_VALID
......
......@@ -11,77 +11,130 @@
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notice for more information.
=========================================================================*/
#include <vtkSmartPointer.h>
#include <vtkImplicitPolyDataDistance.h>
#include <vtkPlaneSource.h>
#include "vtkActor.h"
#include "vtkCamera.h"
#include "vtkGlyph3D.h"
#include "vtkImplicitPolyDataDistance.h"
#include "vtkNew.h"
#include "vtkPlaneSource.h"
#include "vtkPolyDataMapper.h"
#include "vtkProperty.h"
#include "vtkRegressionTestImage.h"
#include "vtkRenderer.h"
#include "vtkRenderWindow.h"
#include "vtkRenderWindowInteractor.h"
#include "vtkSmartPointer.h"
#include "vtkSphereSource.h"
#include "vtkTesting.h"
#include "vtkXMLPolyDataReader.h"
#include <vector>
int TestImplicitPolyDataDistance(int, char*[])
int TestImplicitPolyDataDistance(int argc, char* argv[])
{
vtkSmartPointer<vtkPlaneSource> plane =
vtkSmartPointer<vtkPlaneSource>::New();
plane->SetXResolution(5);
plane->SetYResolution(5);
plane->Update();
vtkSmartPointer<vtkImplicitPolyDataDistance> distance =
vtkSmartPointer<vtkImplicitPolyDataDistance>::New();
distance->SetInput(plane->GetOutput());
// Evaluate at a few points
std::vector<double *> probes;
{
double *probe = new double[3];
probe[0] = probe[1] = 0.0;
probe[2] = 1.0;
probes.push_back(probe);
}
{
double *probe = new double[3];
probe[0] = probe[1] = 0.0;
probe[2] = -1.0;
probes.push_back(probe);
}
{
double *probe = new double[3];
probe[0] = 1.0;
probe[1] = 1.0;
probe[2] = 1.0;
probes.push_back(probe);
}
// Probe at each of the polydata points
for (int i = 0; i < plane->GetOutput()->GetNumberOfPoints(); ++i)
vtkSmartPointer<vtkTesting> testHelper =
vtkSmartPointer<vtkTesting>::New();
testHelper->AddArguments(argc, argv);
if (!testHelper->IsFlagSpecified("-D"))
{
double *point = new double[3];
plane->GetOutput()->GetPoint(i, point);
probes.push_back(point);
std::cerr << "Error: -D /path/to/data was not specified.";
return EXIT_FAILURE;
}
std::vector<double*>::iterator it = probes.begin();
for ( ; it != probes.end(); ++it)
std::string dataRoot = testHelper->GetDataRoot();
std::string fileName = dataRoot + "/Data/CuspySurface.vtp";
std::cout << fileName << std::endl;
// Set up reader
vtkNew<vtkXMLPolyDataReader> reader;
reader->SetFileName(fileName.c_str());
reader->Update();
// Set up distance calculator
vtkNew<vtkImplicitPolyDataDistance> implicitDistance;
implicitDistance->SetInput(reader->GetOutput());
// Compute distances to test points, saving those within the cuspy surface for display
vtkNew<vtkPoints> points;
double xRange[2] = {-47.6, 46.9};
double yRange[2] = {-18.2, 82.1};
double zRange[2] = {1.63, 102};
const double spacing = 10.0;
for (double z = zRange[0]; z < zRange[1]; z += spacing)
{
double *p = *it;
double result = distance->EvaluateFunction(p);
double gradient[3];
distance->EvaluateGradient(p, gradient);
std::cout << "(" << p[0] << ", " << p[1] << ", " << p[2] << "): "
<< result
<< " gradient(" << gradient[0] << ", "
<< gradient[1] << ", "
<< gradient[2] << ")"
<< std::endl;
for (double y = yRange[0]; y < yRange[1]; y += spacing)
{
for (double x = xRange[0]; x < xRange[1]; x += spacing)
{
double pt[3] = {x, y, z};
double distance = implicitDistance->EvaluateFunction(pt);
if (distance <= 0.0)
{
points->InsertNextPoint(pt);
}
}
}
}
distance->Print(std::cout);
// Set up the point data structure
vtkNew<vtkPolyData> insidePoints;
insidePoints->SetPoints(points.GetPointer());
// Glyph the points
vtkNew<vtkSphereSource> sphere;
sphere->SetRadius(3);
vtkNew<vtkGlyph3D> glypher;
glypher->SetInputData(insidePoints.GetPointer());
glypher->SetSourceConnection(sphere->GetOutputPort());
// Display the glyphs
vtkNew<vtkPolyDataMapper> mapper;
mapper->SetInputConnection(glypher->GetOutputPort());
vtkNew<vtkActor> actor;
actor->SetMapper(mapper.GetPointer());
actor->GetProperty()->SetColor(1.0, 0.0, 0.0);
// Display the bounding surface
vtkNew<vtkPolyDataMapper> surfaceMapper;
surfaceMapper->SetInputConnection(reader->GetOutputPort());
it = probes.begin();
for ( ; it != probes.end(); ++it)
vtkNew<vtkActor> surfaceActor;
surfaceActor->SetMapper(surfaceMapper.GetPointer());
surfaceActor->GetProperty()->FrontfaceCullingOn();
// Standard rendering classes
vtkSmartPointer<vtkRenderer> renderer =
vtkSmartPointer<vtkRenderer>::New();
vtkSmartPointer<vtkRenderWindow> renWin =
vtkSmartPointer<vtkRenderWindow>::New();
renWin->SetMultiSamples(0);
renWin->AddRenderer(renderer);
vtkSmartPointer<vtkRenderWindowInteractor> iren =
vtkSmartPointer<vtkRenderWindowInteractor>::New();
iren->SetRenderWindow(renWin);
renderer->AddActor(actor.GetPointer());
renderer->AddActor(surfaceActor.GetPointer());
// Standard testing code.
renderer->SetBackground(0.0, 0.0, 0.0);
renWin->SetSize(300,300);
vtkCamera *camera = renderer->GetActiveCamera();
renderer->ResetCamera();
camera->Azimuth(30);
camera->Elevation(-20);
iren->Initialize();
renWin->Render();
int retVal = vtkRegressionTestImage(renWin.GetPointer());
if (retVal == vtkRegressionTester::DO_INTERACTOR)
{
delete [] *it;
iren->Start();
}
return EXIT_SUCCESS;
return !retVal;
}
......@@ -164,12 +164,10 @@ double vtkImplicitPolyDataDistance::SharedEvaluate(double x[3], double n[3])
{
count += (fabs(weights[i]) < this->Tolerance ? 1 : 0);
}
// if weights contains no 0s
if ( count == 0 || count == 1 )
// Face case - weights contains no 0s
if ( count == 0 )
{
// Compute face normal.
// For count == 0, this is all we need.
// For count = 1, we'll add in the normals from adjacent faces.
if ( cnorms )
{
cnorms->GetTuple(cellId, awnorm);
......@@ -179,9 +177,8 @@ double vtkImplicitPolyDataDistance::SharedEvaluate(double x[3], double n[3])
vtkPolygon::ComputeNormal(cell->Points, awnorm);
}
}
// if weights contains 1 0s
if ( count == 1 )
// Edge case - weights contain one 0
else if ( count == 1 )
{
// ... edge ... get two adjacent faces, compute average normal
int a = -1, b = -1;
......@@ -202,7 +199,9 @@ double vtkImplicitPolyDataDistance::SharedEvaluate(double x[3], double n[3])
return this->NoValue;
}
this->Input->GetCellEdgeNeighbors(0, a, b, idList);
// The first argument is the cell ID. We pass a bogus cell ID so that
// all face IDs attached to the edge are returned in the idList.
this->Input->GetCellEdgeNeighbors(VTK_ID_MAX, a, b, idList);
for (int i = 0; i < idList->GetNumberOfIds(); i++)
{
double norm[3];
......@@ -221,7 +220,7 @@ double vtkImplicitPolyDataDistance::SharedEvaluate(double x[3], double n[3])
vtkMath::Normalize(awnorm);
}
// If weights contains 2 0s
// Vertex case - weights contain two 0s
else if ( count == 2 )
{
// ... vertex ... this is the expensive case, get all adjacent
......
b78e62003c1960033b7b5e80c4b3e38a
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