The problem you are facing is with the precision of floating point numbers. The
way you are setting up your problem is that the points are exactly 0.01 units
away from each other. When you ask the point locator for all points within a
radius of 0.01, then what happens when a point is exactly 0.01 units away from
the center? That depends on what happens in the limited precision of your
computer's floating point values. Sometimes it will be 0.01 + epsilon, in which
case the point will be considered inside the radius. Sometimes it will be 0.01
- epsilon, in which case the point will be considered outside the radius. This
is a well known issue with floating point numbers. There are lots of resources
describing it. A quick Google search gave this article:
https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html.
Anyway, the simple solution is to increase the radius in FindPointsWithinRadius
by a small amount. Changing the line to:
loc.FindPointsWithinRadius(0.01,xp,id_list)
fixes the problem for all point locators.
BTW, although perhaps convenient, using a point locator to find neighbors is
not the most robust way to do it. It would be safer to use the find neighbor
facilities of the different vtkDataSet classes.
-Ken
From: ParaView [mailto:[email protected]] On Behalf Of Bertrand
Gazanion
Sent: Thursday, June 23, 2016 8:13 AM
To: [email protected]
Subject: [EXTERNAL] Re: [Paraview] Find points within radius
Forgot the pictures... Sorry for the spam.
Bertrand
De : Bertrand Gazanion [mailto:[email protected]]
Envoyé : jeudi 23 juin 2016 16:11
À : '[email protected]'
Objet : Find points within radius
On behalf of a colleague of mine working with ParaView 4.3.1
Hello,
I am using an algorithm doing calculation on points and their neighbours and I
find the results of those points locators to be strange. To test it I made a
simple algorithm counting how many neighbours each point has.
When I create a plane and set its x and y resolution to 100 each, all points
should have the same number of neighbours except on the sides. However
vtkPointLocator, vtkKdTreePointLocator and vtkOctreePointLocator do not find
the same number of neighbours for each point.
Here is the algorithm that I use to get the number of neighbors each points has
:
import vtk
pdi = self.GetInputDataObject(0,0)
nb_pts = pdi.GetNumberOfPoints()
# Point locator
loc = vtk.vtkPointLocator()
#loc = vtk.vtkOctreePointLocator()
#loc = vtk. vtkKdTreePointLocator()
loc.SetDataSet(pdi)
loc.BuildLocator()
pts = pdi.GetPoints()
neighbours = vtk.vtkTypeInt64Array()
neighbours.SetNumberOfComponents(1)
neighbours.SetNumberOfTuples(nb_pts)
neighbours.SetName('neighbours')
for k in range(nb_pts):
xp = pdi.GetPoint(k)
id_list = vtk.vtkIdList()
loc.FindPointsWithinRadius(0.01,xp,id_list)
loc.Update()
N = id_list.GetNumberOfIds()
neighbours.InsertTuple1(k,N)
self.GetOutput().GetPointData().AddArray(neighbours)
You will find attached sample picture describing the result I get with each
point locator.
thank you very much
_______________________________________________
Powered by www.kitware.com
Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html
Please keep messages on-topic and check the ParaView Wiki at:
http://paraview.org/Wiki/ParaView
Search the list archives at: http://markmail.org/search/?q=ParaView
Follow this link to subscribe/unsubscribe:
http://public.kitware.com/mailman/listinfo/paraview
