[Paraview] Find points within radius

[Moreland, Kenneth]

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:paraview-bounces at paraview.org] On Behalf Of Bertrand Gazanion
Sent: Thursday, June 23, 2016 8:13 AM
To: paraview at paraview.org
Subject: [EXTERNAL] Re: [Paraview] Find points within radius

Forgot the pictures... Sorry for the spam.
Bertrand

De : Bertrand Gazanion [mailto:b.gazanion at gantha.com]
Envoyé : jeudi 23 juin 2016 16:11
À : 'paraview at paraview.org'
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/paraview/attachments/20160623/dca87e8b/attachment.html>