<html xmlns:v="urn:schemas-microsoft-com:vml" xmlns:o="urn:schemas-microsoft-com:office:office" xmlns:w="urn:schemas-microsoft-com:office:word" xmlns:m="http://schemas.microsoft.com/office/2004/12/omml" xmlns="http://www.w3.org/TR/REC-html40">
<head>
<meta http-equiv=Content-Type content="text/html; charset=iso-8859-1">
<meta name=Generator content="Microsoft Word 12 (filtered medium)">
<style>
<!--
/* Font Definitions */
@font-face
{font-family:Helvetica;
panose-1:2 11 5 4 2 2 2 3 2 4;}
@font-face
{font-family:"Cambria Math";
panose-1:2 4 5 3 5 4 6 3 2 4;}
@font-face
{font-family:Calibri;
panose-1:2 15 5 2 2 2 4 3 2 4;}
@font-face
{font-family:Tahoma;
panose-1:2 11 6 4 3 5 4 4 2 4;}
@font-face
{font-family:"Segoe UI";
panose-1:2 11 5 2 4 2 4 2 2 3;}
/* Style Definitions */
p.MsoNormal, li.MsoNormal, div.MsoNormal
{margin:0cm;
margin-bottom:.0001pt;
font-size:11.0pt;
font-family:"Calibri","sans-serif";}
a:link, span.MsoHyperlink
{mso-style-priority:99;
color:blue;
text-decoration:underline;}
a:visited, span.MsoHyperlinkFollowed
{mso-style-priority:99;
color:purple;
text-decoration:underline;}
p.MsoAcetate, li.MsoAcetate, div.MsoAcetate
{mso-style-priority:99;
mso-style-link:"Texte de bulles Car";
margin:0cm;
margin-bottom:.0001pt;
font-size:11.0pt;
font-family:"Calibri","sans-serif";}
span.TextedebullesCar
{mso-style-name:"Texte de bulles Car";
mso-style-priority:99;
mso-style-link:"Texte de bulles";
font-family:"Tahoma","sans-serif";}
p.msonormal0, li.msonormal0, div.msonormal0
{mso-style-name:msonormal;
mso-margin-top-alt:auto;
margin-right:0cm;
mso-margin-bottom-alt:auto;
margin-left:0cm;
font-size:12.0pt;
font-family:"Times New Roman","serif";}
p.BalloonText, li.BalloonText, div.BalloonText
{mso-style-name:"Balloon Text";
mso-style-link:"Balloon Text Char";
margin:0cm;
margin-bottom:.0001pt;
font-size:11.0pt;
font-family:"Calibri","sans-serif";}
span.BalloonTextChar
{mso-style-name:"Balloon Text Char";
mso-style-priority:99;
mso-style-link:"Balloon Text";
font-family:"Segoe UI","sans-serif";}
span.EmailStyle22
{mso-style-type:personal;
font-family:"Calibri","sans-serif";
color:windowtext;}
span.EmailStyle23
{mso-style-type:personal;
font-family:"Calibri","sans-serif";
color:black;}
span.EmailStyle24
{mso-style-type:personal;
font-family:"Calibri","sans-serif";
color:#1F497D;}
span.EmailStyle25
{mso-style-type:personal-reply;
font-family:"Calibri","sans-serif";
color:black;}
.MsoChpDefault
{mso-style-type:export-only;
font-size:10.0pt;}
@page Section1
{size:612.0pt 792.0pt;
margin:70.85pt 70.85pt 70.85pt 70.85pt;}
div.Section1
{page:Section1;}
-->
</style>
<!--[if gte mso 9]><xml>
<o:shapedefaults v:ext="edit" spidmax="1026" />
</xml><![endif]--><!--[if gte mso 9]><xml>
<o:shapelayout v:ext="edit">
<o:idmap v:ext="edit" data="1" />
</o:shapelayout></xml><![endif]-->
</head>
<body lang=FR link=blue vlink=purple>
<div class=Section1>
<p class=MsoNormal><i><span lang=EN-US style='color:black'>Dear Kenneth, thank
you for your answer. Here is the message from my colleague.<o:p></o:p></span></i></p>
<p class=MsoNormal><i><span lang=EN-US style='color:black'>Bertrand<o:p></o:p></span></i></p>
<p class=MsoNormal><span lang=EN-US style='color:black'><o:p> </o:p></span></p>
<p class=MsoNormal><span lang=EN-US style='color:black'>Hello, I am the
colleague using ParaView 4.3.1!<o:p></o:p></span></p>
<p class=MsoNormal><span lang=EN-US style='color:black'><o:p> </o:p></span></p>
<p class=MsoNormal><span lang=EN-US style='color:black'>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.<o:p></o:p></span></p>
<p class=MsoNormal><span lang=EN-US style='color:black'><o:p> </o:p></span></p>
<p class=MsoNormal><span lang=EN-US style='color:black'>Unfortunately the
physics behind my algorithm forces me to look for all the neighbors inside a
sphere of each point of my unstructured grid.<o:p></o:p></span></p>
<p class=MsoNormal><span lang=EN-US style='color:black'><o:p> </o:p></span></p>
<p class=MsoNormal><span lang=EN-US style='color:black'>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? <o:p></o:p></span></p>
<p class=MsoNormal><span lang=EN-US style='color:black'>I have read a more
straightforward </span><span style='color:black'><a
href="http://floating-point-gui.de/"><span lang=EN-US>website</span></a></span><span
lang=EN-US style='color:black'> about floating point number : </span><span
style='color:black'><a href="http://floating-point-gui.de/" target="_blank"><span
lang=EN-US>http://floating-point-gui.de/</span></a></span><span lang=EN-US
style='color:black'> which I hope offer good information.<o:p></o:p></span></p>
<p class=MsoNormal><span lang=EN-US style='color:black'>The page about </span><span
style='color:black'><a href="http://floating-point-gui.de/errors/comparison/"><span
lang=EN-US>comparison</span></a></span><span lang=EN-US style='color:black'>
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.<o:p></o:p></span></p>
<p class=MsoNormal><span lang=EN-US style='color:black'><o:p> </o:p></span></p>
<p class=MsoNormal><span lang=EN-US style='color:black'>I am not sure if it is
vtk responsibility to find out how to deal with those case but I find the
solution interesting.<o:p></o:p></span></p>
<p class=MsoNormal><span lang=EN-US style='color:black'><o:p> </o:p></span></p>
<p class=MsoNormal><span style='color:black'>Thank you!<o:p></o:p></span></p>
<p class=MsoNormal><span style='color:black'><o:p> </o:p></span></p>
<p class=MsoNormal><span style='color:black'><o:p> </o:p></span></p>
<div>
<div style='border:none;border-top:solid #B5C4DF 1.0pt;padding:3.0pt 0cm 0cm 0cm'>
<p class=MsoNormal><b><span style='font-size:10.0pt;font-family:"Tahoma","sans-serif"'>From:</span></b><span
style='font-size:10.0pt;font-family:"Tahoma","sans-serif"'> Moreland, Kenneth
[mailto:kmorel@sandia.gov] <br>
<b>Sent:</b> Thursday, June 23, 2016 5:11 PM<br>
<b>To:</b> Bertrand Gazanion; paraview@paraview.org<br>
<b>Subject:</b> Re: [Paraview] Find points within radius<o:p></o:p></span></p>
</div>
</div>
<p class=MsoNormal><o:p> </o:p></p>
<p class=MsoNormal><span lang=EN-US style='color:#1F497D'>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: <a
href="https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html">https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html</a>.<o:p></o:p></span></p>
<p class=MsoNormal><span lang=EN-US style='color:#1F497D'><o:p> </o:p></span></p>
<p class=MsoNormal><span lang=EN-US style='color:#1F497D'>Anyway, the simple
solution is to increase the radius in FindPointsWithinRadius by a small amount.
Changing the line to:<o:p></o:p></span></p>
<p class=MsoNormal><span lang=EN-US style='color:#1F497D'><o:p> </o:p></span></p>
<p class=MsoNormal><span lang=EN-US style='font-family:"Courier New";
color:#1F497D'> loc.FindPointsWithinRadius(0.01,xp,id_list)<o:p></o:p></span></p>
<p class=MsoNormal><span lang=EN-US style='color:#1F497D'><o:p> </o:p></span></p>
<p class=MsoNormal><span lang=EN-US style='color:#1F497D'>fixes the problem for
all point locators.<o:p></o:p></span></p>
<p class=MsoNormal><span lang=EN-US style='color:#1F497D'><o:p> </o:p></span></p>
<p class=MsoNormal><span lang=EN-US style='color:#1F497D'>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.<o:p></o:p></span></p>
<p class=MsoNormal><span lang=EN-US style='color:#1F497D'><o:p> </o:p></span></p>
<p class=MsoNormal><span lang=EN-US style='color:#1F497D'>-Ken<o:p></o:p></span></p>
<p class=MsoNormal><span lang=EN-US style='color:#1F497D'><o:p> </o:p></span></p>
<div>
<div style='border:none;border-top:solid #E1E1E1 1.0pt;padding:3.0pt 0cm 0cm 0cm'>
<p class=MsoNormal><b><span lang=EN-US>From:</span></b><span lang=EN-US>
ParaView [mailto:paraview-bounces@paraview.org] <b>On Behalf Of </b>Bertrand
Gazanion<br>
<b>Sent:</b> Thursday, June 23, 2016 8:13 AM<br>
<b>To:</b> paraview@paraview.org<br>
<b>Subject:</b> [EXTERNAL] Re: [Paraview] Find points within radius<o:p></o:p></span></p>
</div>
</div>
<p class=MsoNormal><span lang=EN-US><o:p> </o:p></span></p>
<p class=MsoNormal><span lang=EN-US style='color:black'>Forgot the pictures...
Sorry for the spam.<o:p></o:p></span></p>
<p class=MsoNormal><span lang=EN-US style='color:black'>Bertrand<o:p></o:p></span></p>
<p class=MsoNormal><span lang=EN-US style='color:black'><o:p> </o:p></span></p>
<div>
<div style='border:none;border-top:solid #B5C4DF 1.0pt;padding:3.0pt 0cm 0cm 0cm'>
<p class=MsoNormal><b><span lang=EN-US style='font-size:10.0pt;font-family:
"Tahoma","sans-serif"'>De :</span></b><span lang=EN-US style='font-size:
10.0pt;font-family:"Tahoma","sans-serif"'> Bertrand Gazanion [<a
href="mailto:b.gazanion@gantha.com">mailto:b.gazanion@gantha.com</a>] <br>
<b>Envoyé :</b> jeudi 23 juin 2016 16:11<br>
<b>À :</b> 'paraview@paraview.org</span><span style='font-size:10.0pt;
font-family:"Tahoma","sans-serif"'>'<br>
<b>Objet :</b> Find points within radius<o:p></o:p></span></p>
</div>
</div>
<p class=MsoNormal><o:p> </o:p></p>
<p class=MsoNormal style='background:white'><i><span lang=EN-US
style='font-family:"Helvetica","sans-serif";color:black'>On behalf of a
colleague of mine working with ParaView 4.3.1<o:p></o:p></span></i></p>
<p class=MsoNormal style='background:white'><span lang=EN-US style='font-family:
"Helvetica","sans-serif";color:black'><o:p> </o:p></span></p>
<p class=MsoNormal style='background:white'><span lang=EN-US style='font-family:
"Helvetica","sans-serif";color:black'>Hello,<o:p></o:p></span></p>
<p class=MsoNormal style='background:white'><span lang=EN-US style='font-family:
"Helvetica","sans-serif";color:black'><br>
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.<br
id="yui_3_16_0_ym19_1_1466689264692_2872">
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.<br
id="yui_3_16_0_ym19_1_1466689264692_2873">
<br id="yui_3_16_0_ym19_1_1466689264692_2874">
Here is the algorithm that I use to get the number of neighbors each points has
:<br id="yui_3_16_0_ym19_1_1466689264692_2875">
</span><span lang=EN-US style='font-size:10.0pt;font-family:"Helvetica","sans-serif";
color:black'>import vtk<br id="yui_3_16_0_ym19_1_1466689264692_2876">
<br id="yui_3_16_0_ym19_1_1466689264692_2877">
pdi = self.GetInputDataObject(0,0)<br id="yui_3_16_0_ym19_1_1466689264692_2878">
nb_pts = pdi.GetNumberOfPoints()<br id="yui_3_16_0_ym19_1_1466689264692_2879">
<br id="yui_3_16_0_ym19_1_1466689264692_2880">
# Point locator<br id="yui_3_16_0_ym19_1_1466689264692_2881">
loc = vtk.vtkPointLocator()<br id="yui_3_16_0_ym19_1_1466689264692_2882">
#loc = vtk.vtkOctreePointLocator()<br id="yui_3_16_0_ym19_1_1466689264692_2883">
#loc = vtk. vtkKdTreePointLocator()<br id="yui_3_16_0_ym19_1_1466689264692_2884">
<br id="yui_3_16_0_ym19_1_1466689264692_2885">
loc.SetDataSet(pdi)<br id="yui_3_16_0_ym19_1_1466689264692_2886">
loc.BuildLocator()<br id="yui_3_16_0_ym19_1_1466689264692_2887">
<br id="yui_3_16_0_ym19_1_1466689264692_2888">
pts = pdi.GetPoints()<br id="yui_3_16_0_ym19_1_1466689264692_2889">
<br id="yui_3_16_0_ym19_1_1466689264692_2890">
neighbours = vtk.vtkTypeInt64Array()<br
id="yui_3_16_0_ym19_1_1466689264692_2891">
neighbours.SetNumberOfComponents(1)<br id="yui_3_16_0_ym19_1_1466689264692_2892">
neighbours.SetNumberOfTuples(nb_pts)<br
id="yui_3_16_0_ym19_1_1466689264692_2893">
neighbours.SetName('neighbours')<br id="yui_3_16_0_ym19_1_1466689264692_2894">
<br id="yui_3_16_0_ym19_1_1466689264692_2895">
for k in range(nb_pts):<br id="yui_3_16_0_ym19_1_1466689264692_2896">
<br id="yui_3_16_0_ym19_1_1466689264692_2897">
xp = pdi.GetPoint(k)<br
id="yui_3_16_0_ym19_1_1466689264692_2898">
<br id="yui_3_16_0_ym19_1_1466689264692_2899">
id_list = vtk.vtkIdList()<br
id="yui_3_16_0_ym19_1_1466689264692_2900">
<br id="yui_3_16_0_ym19_1_1466689264692_2901">
loc.FindPointsWithinRadius(0.01,xp,id_list)<br
id="yui_3_16_0_ym19_1_1466689264692_2902">
loc.Update()<br id="yui_3_16_0_ym19_1_1466689264692_2903">
<br id="yui_3_16_0_ym19_1_1466689264692_2904">
N = id_list.GetNumberOfIds()<br
id="yui_3_16_0_ym19_1_1466689264692_2905">
neighbours.InsertTuple1(k,N)<br
id="yui_3_16_0_ym19_1_1466689264692_2906">
<br id="yui_3_16_0_ym19_1_1466689264692_2907">
<br id="yui_3_16_0_ym19_1_1466689264692_2908">
self.GetOutput().GetPointData().AddArray(neighbours)</span><span lang=EN-US
style='font-family:"Helvetica","sans-serif";color:black'><o:p></o:p></span></p>
<p class=MsoNormal style='background:white'><span lang=EN-US style='font-family:
"Helvetica","sans-serif";color:black'><o:p> </o:p></span></p>
<p class=MsoNormal style='background:white'><span lang=EN-US style='font-family:
"Helvetica","sans-serif";color:black'>You will find attached sample picture
describing the result I get with each point locator.<br>
<br id="yui_3_16_0_ym19_1_1466689264692_2911">
thank you very much<o:p></o:p></span></p>
</div>
</body>
</html>