[Paraview] Find points within radius

[Bertrand Gazanion]

Dear Kenneth, thank you for your answer. Here is the message from my
colleague.

Bertrand

 

Hello, I am the colleague using ParaView 4.3.1!

 

First of all, thank you for your response, it has explained a lot of results
since I tend to forget about computer and floating point.

 

Unfortunately the physics behind my algorithm forces me to look for all the
neighbors inside a sphere of each point of my unstructured grid.

 

My curiosity has pushed me to look how the FindPointsWithinRadius function
works and I have seen that it does an inferior or equal comparison. Isn't
the equal part the source of inconsistent return values? 

I have read a more straightforward  <http://floating-point-gui.de/> website
about floating point number :  <http://floating-point-gui.de/>
http://floating-point-gui.de/  which I hope offer good information.

The page about  <http://floating-point-gui.de/errors/comparison/> comparison
offer a way to check if two numbers are equal that seems to be accurate, I
have tried to implement it in the FindPointsWithinRadius function of
vtkPointLocator but I am not proficient at c++, vtk and your coding
guidelines, see attachment.

 

I am not sure if it is vtk responsibility to find out how to deal with those
case but I find the solution interesting.

 

Thank you!

 

 

From: Moreland, Kenneth [mailto:kmorel at sandia.gov] 
Sent: Thursday, June 23, 2016 5:11 PM
To: Bertrand Gazanion; paraview at paraview.org
Subject: Re: [Paraview] Find points within radius

 

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/20160624/cdd74b01/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: FindPointsWithinRadius.cxx
Type: application/octet-stream
Size: 2066 bytes
Desc: not available
URL: <http://public.kitware.com/pipermail/paraview/attachments/20160624/cdd74b01/attachment.obj>