Ordering 3dpoints ?

Didier Bur · 18 Jan 2011, 22:11

Hi,
Anyone knows a way to order 4 points in a clockwise or ccw manner, given the 4 points and assuming they are coplanar ?
I'm actually generating faces that are "butterfly shaped", and the only way I've found to detect these sort of face is to test if two of its edges intersect.

TIG · 18 Jan 2011, 22:28

@didier bur said:

Hi,
Anyone knows a way to order 4 points in a clockwise or ccw manner, given the 4 points and assuming their are coplanar ?
I'm actually generating faces that are "butterfly shaped", and the only way I've found to detect these sort of face is to test if two of its edges intersect.

Yoy should find that faces.outer_loop.vertices always return ccwise [about face.normal] and all other (face.loops-[face.outer_loop]).each{loop|loop.vertices} [holes etc] are cwise. If they are loop vertices they must be coplanar... Why not quickly make a grouped face using the points and get the vertices back to points ?
I deduce you want to take any 4 coplanar points and make a face from them that's NOT a 'bow-tie' ?
Think of it another way...
Work inside a group's entities.
If the points are coplanar any three can define a plane.
Join three with edges.
Take the last point and add edges from it to the three other points that are the nearest.
Now find faces for all of these 6 edges.
Now two of these edges has two faces - erase them.
You now have a quad that's not 'bow-tied' ?

fredo6 · 18 Jan 2011, 22:43

@didier bur said:

Hi,
Anyone knows a way to order 4 points in a clockwise or ccw manner, given the 4 points and assuming they are coplanar ?
I'm actually generating faces that are "butterfly shaped", and the only way I've found to detect these sort of face is to test if two of its edges intersect.

Didier,

If points P1, P2, P3, P4 are coplanar, then, I think you can

compute the barycenter G of the 4 points
sort the points P2, P3, P4 by the angle between G->P1 and G->Pi

Fredo

Cleverbeans · 18 Jan 2011, 22:44

First it's important to note that such a face may not be unique. Say for example if you take the points (0,0),(1,1),(-1,1) and (0,-1) then every ordering that starts with (0,0) will either be clockwise or counter-clockwise but depending on way you put the other points in will get a different face. If you can always count on the face being convex however I believe this ambiguity vanishes.

To order them appropriately, you can first choose any three points, and just put them in any order, which will either be clockwise or counter-clockwise. Now you want to add the fourth point into the order at the correct spot so that you don't get that intersection issue. You can convince yourself pretty easily that if you draw a triangle, and take a point outside that triangle you basically want to connect it to two points without crossing any of the existing lines on the triangle. From here, it's easy to see that you want to connect it to the two points which are closest to it. Here is some code which will return the points in the desired order under the conditions given.


def order_points(points)
   #Given an array of four coplanar, convex point3d objects,  returns the points in convex order

   pt = points[0]
   triangle = points[1..-1] #the triangle which is already either clockwise or counter-clockwise

   distances = points.map{|x| pt.distance_to(x)}
   i = distances.index(distances.max()) #this finds the index of the point which is furthest away in the points array

   if i == 0 then
       return triangle[0..1] + [pt] + triangle[2]
   elsif i==1 then
       return triangle + [pt]
   else
       return [pt] + triangle
   end
end

I think that should do the trick.

Didier Bur · 18 Jan 2011, 22:53

Me again,
Thanks all for your answers,
@TIG: I'd better use a 'theoric' way, instead of using SketchUp (that slows down the computation too much)
@Fredo6: why didn't I thought of that
Cleverbeans has a clever solution too

Didier Bur · 20 Jan 2011, 20:47

I've implemented CleverBeans solution to order the points but found that it does not work:

elsif i==1 then
           return triangle + [pt]
       else
           return [pt] + triangle

two cases that creates exactly the same quad face is not cool

Here is my slightly different coding:

# Given a set of 4 coplanar 3dPoints, returns the points in convex order
    def quadOrderPoints(pArray)
		pt = pArray[0]
		triangle = pArray[1..-1] #the triangle which is already either clockwise or counter-clockwise
		distances =triangle.map{ |x| pt.distance(x) }
		#Find the index of the point which is furthest away in the points array
		i = distances.index(distances.max())
		case i
			when 0
				return [pArray[1],pArray[2],pArray[0],pArray[3]]
			when 1
				return [pArray[1],pArray[2],pArray[3],pArray[0]]
			when 2
				return [pArray[0],pArray[1],pArray[2],pArray[3]]
		end
	end

Why it doesn't work: here on the image, triangle is p1,p2,p3, and p0 is the first point in the array to order.
Greatest distance from p0 to one of the triangle vertices is shown in red.
Left figure shows, as 'distances.max' index will be 0, that quad will be (according to the algorithm): p1,p2,p0,p3
Right figure shows that 'distances.max index will be 1, so the quad will be (according to the algorithm): p1,p2,p3,P0, or p0,p1,p2,p3.

BUT both figures suggest that the quad has to be p0,p1,p2,p3.
I've tested it and about half of the quads are wrong.

Dan Rathbun · 20 Jan 2011, 22:16

What about about using a temporary BouningBox around the 4 pts.
Each pt would lie on one of the BB's planes. 2 of the BB's planes would not have a point.
Either the vertical planes will have the pts, or the top and bottom and 2 side planes will have the points.

Chris Fullmer · 20 Jan 2011, 22:44

I can only manage to describe this in 2d, so you get to implement it onto any 3d plane. I tihnk it is similar to what Fredo suggested - his is probably faster. Mine requires re-calculating on each point. But mine works to find the convex face that bounds MANY points if needed.

Get the point in the most negative x location. Then get the angle between the vector -1,0,0 and the vector from your point to all other points. The smallest angle is the next point in the convex face. You have to make sure you are always getting the angle going in the correct direction (clockwise or counter-clockwise). Then use the previous vector as the new base vector when you go to test for the next point.

I made a video that shows it. Think of the protractor as doing the testing. It rotates until it finds the next point to use - it is the point with the smallest angle.

There is no sound.
[flash=638,509:26c6u44v]http://www.chrisfullmer.com/forums/convex_face_finder.swf[/flash:26c6u44v]

Jim · 20 Jan 2011, 23:19

Start by looking at the cross product of each vertex's edge vectors.


model = Sketchup.active_model
entities = model.entities
selection = model.selection

pts = [  [0, 0], [0, 1], [1, 1], [1, 0] ]
pts2 = [ pts[0], pts[2], pts[1], pts[3] ]

g1 = entities.add_group
g1.entities.add_face pts
pts.each_with_index{|pt, i| g1.entities.add_text((i+1).to_s, pt)}
g1.name = "Normal"

g2 = entities.add_group
g2.entities.add_face pts2
g2.name = "Swapped"
pts2.each_with_index{|pt, i| g2.entities.add_text((i+1).to_s, pt)}

def show_cross(ents, ary)
  p1, p2, p3, p4 = ary
  n1 = ((p1.vector_to p4) * (p1.vector_to p2)).normalize
  n2 = ((p2.vector_to p1) * (p2.vector_to p3)).normalize
  n3 = ((p3.vector_to p2) * (p3.vector_to p4)).normalize
  n4 = ((p4.vector_to p3) * (p4.vector_to p1)).normalize
  ents.add_line(p1, p1.offset(n1))
  ents.add_line(p2, p2.offset(n2))
  ents.add_line(p3, p3.offset(n3))
  ents.add_line(p4, p4.offset(n4))
end

show_cross(g1.entities, pts)
show_cross(g2.entities, pts2)

TIG · 21 Jan 2011, 09:16

http://en.wikipedia.org/wiki/Graham_scan ?
Ensures 'convex' shapes...

thomthom · 21 Jan 2011, 10:24

@tig said:

http://en.wikipedia.org/wiki/Graham_scan ?
Ensures 'convex' shapes...

I got an implementation of Grahams Scan in TT_Lib (version 1)
From geom2d.rb

<span class="syntaxdefault"><br /></span><span class="syntaxcomment"># Graham's Scan<br /># Computational Geometry Algorythms and Appliations (ISBN-10; 3540779736)<br /># Chapter 1.1<br />#<br /># @param [Array<Geom;;Point3d>] points<br /># @return [Array<Geom;;Point3d>]<br /># @since 1.1.0<br /></span><span class="syntaxdefault">def self</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">convex_hull</span><span class="syntaxkeyword">(</span><span class="syntaxdefault">points</span><span class="syntaxkeyword">)<br /></span><span class="syntaxdefault">  return nil if points</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">length </span><span class="syntaxkeyword"><</span><span class="syntaxdefault"> 2<br />  <br />  points </span><span class="syntaxkeyword">=</span><span class="syntaxdefault"> self</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">sort_points_by_x_y</span><span class="syntaxkeyword">(</span><span class="syntaxdefault">points</span><span class="syntaxkeyword">)<br /></span><span class="syntaxdefault">  <br />  l_upper </span><span class="syntaxkeyword">=</span><span class="syntaxdefault"> </span><span class="syntaxkeyword">[</span><span class="syntaxdefault"> points</span><span class="syntaxkeyword">[</span><span class="syntaxdefault">0</span><span class="syntaxkeyword">],</span><span class="syntaxdefault"> points</span><span class="syntaxkeyword">[</span><span class="syntaxdefault">1</span><span class="syntaxkeyword">]</span><span class="syntaxdefault"> </span><span class="syntaxkeyword">]<br /></span><span class="syntaxdefault">  2</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">upto</span><span class="syntaxkeyword">(</span><span class="syntaxdefault">points</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">length </span><span class="syntaxkeyword">-</span><span class="syntaxdefault"> 1</span><span class="syntaxkeyword">)</span><span class="syntaxdefault"> </span><span class="syntaxkeyword">{</span><span class="syntaxdefault"> </span><span class="syntaxkeyword">|</span><span class="syntaxdefault">i</span><span class="syntaxkeyword">|<br /></span><span class="syntaxdefault">    l_upper </span><span class="syntaxkeyword"><<</span><span class="syntaxdefault"> points</span><span class="syntaxkeyword">[</span><span class="syntaxdefault">i</span><span class="syntaxkeyword">]<br /></span><span class="syntaxdefault">    while l_upper</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">length </span><span class="syntaxkeyword">></span><span class="syntaxdefault"> 2 </span><span class="syntaxkeyword">&&</span><span class="syntaxdefault"> </span><span class="syntaxkeyword">!</span><span class="syntaxdefault">self</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">right_turn</span><span class="syntaxkeyword">?(</span><span class="syntaxdefault">l_upper</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">last</span><span class="syntaxkeyword">(</span><span class="syntaxdefault">3</span><span class="syntaxkeyword">))<br /></span><span class="syntaxdefault">      l_upper</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">delete_at</span><span class="syntaxkeyword">(-</span><span class="syntaxdefault">2</span><span class="syntaxkeyword">)<br /></span><span class="syntaxdefault">    end<br />  </span><span class="syntaxkeyword">}<br /></span><span class="syntaxdefault">  <br />  l_lower </span><span class="syntaxkeyword">=</span><span class="syntaxdefault"> </span><span class="syntaxkeyword">[</span><span class="syntaxdefault"> points</span><span class="syntaxkeyword">[-</span><span class="syntaxdefault">1</span><span class="syntaxkeyword">],</span><span class="syntaxdefault"> points</span><span class="syntaxkeyword">[-</span><span class="syntaxdefault">2</span><span class="syntaxkeyword">]</span><span class="syntaxdefault"> </span><span class="syntaxkeyword">]<br /></span><span class="syntaxdefault">  </span><span class="syntaxkeyword">(</span><span class="syntaxdefault">points</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">length </span><span class="syntaxkeyword">-</span><span class="syntaxdefault"> 3</span><span class="syntaxkeyword">).</span><span class="syntaxdefault">downto</span><span class="syntaxkeyword">(</span><span class="syntaxdefault">0</span><span class="syntaxkeyword">)</span><span class="syntaxdefault"> </span><span class="syntaxkeyword">{</span><span class="syntaxdefault"> </span><span class="syntaxkeyword">|</span><span class="syntaxdefault">i</span><span class="syntaxkeyword">|<br /></span><span class="syntaxdefault">    l_lower </span><span class="syntaxkeyword"><<</span><span class="syntaxdefault"> points</span><span class="syntaxkeyword">[</span><span class="syntaxdefault">i</span><span class="syntaxkeyword">]<br /></span><span class="syntaxdefault">    while l_lower</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">length </span><span class="syntaxkeyword">></span><span class="syntaxdefault"> 2 </span><span class="syntaxkeyword">&&</span><span class="syntaxdefault"> </span><span class="syntaxkeyword">!</span><span class="syntaxdefault">self</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">right_turn</span><span class="syntaxkeyword">?(</span><span class="syntaxdefault">l_lower</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">last</span><span class="syntaxkeyword">(</span><span class="syntaxdefault">3</span><span class="syntaxkeyword">))<br /></span><span class="syntaxdefault">      l_lower</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">delete_at</span><span class="syntaxkeyword">(-</span><span class="syntaxdefault">2</span><span class="syntaxkeyword">)<br /></span><span class="syntaxdefault">    end<br />  </span><span class="syntaxkeyword">}<br /></span><span class="syntaxdefault">  l_lower</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">delete_at</span><span class="syntaxkeyword">(</span><span class="syntaxdefault">0</span><span class="syntaxkeyword">)<br /></span><span class="syntaxdefault">  l_lower</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">delete_at</span><span class="syntaxkeyword">(-</span><span class="syntaxdefault">1</span><span class="syntaxkeyword">)<br /></span><span class="syntaxdefault">  <br />  return l_upper </span><span class="syntaxkeyword">+</span><span class="syntaxdefault"> l_lower<br />end<br /><br /><br /></span><span class="syntaxcomment"># Lexicographic order x, y<br />#<br /># @param [Array<Geom;;Point3d>] points<br /># @return [Array<Geom;;Point3d>]<br /># @since 1.1.0<br /></span><span class="syntaxdefault">def self</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">sort_points_by_x_y</span><span class="syntaxkeyword">(</span><span class="syntaxdefault">points</span><span class="syntaxkeyword">)<br /></span><span class="syntaxdefault">  return points</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">sort </span><span class="syntaxkeyword">{</span><span class="syntaxdefault"> </span><span class="syntaxkeyword">|</span><span class="syntaxdefault">a</span><span class="syntaxkeyword">,</span><span class="syntaxdefault">b</span><span class="syntaxkeyword">|</span><span class="syntaxdefault"> </span><span class="syntaxkeyword">(</span><span class="syntaxdefault">a</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">x</span><span class="syntaxkeyword">==</span><span class="syntaxdefault">b</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">x</span><span class="syntaxkeyword">)</span><span class="syntaxdefault"> </span><span class="syntaxkeyword">?</span><span class="syntaxdefault"> a</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">y </span><span class="syntaxkeyword"><=></span><span class="syntaxdefault"> b</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">y </span><span class="syntaxkeyword">;</span><span class="syntaxdefault"> a</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">x </span><span class="syntaxkeyword"><=></span><span class="syntaxdefault"> b</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">x </span><span class="syntaxkeyword">}<br /></span><span class="syntaxdefault">end<br /></span>

Didier Bur · 21 Jan 2011, 10:45

Thanks TIG and TT, I was just starting coding Graham's scan, this will save me time.
Left turn/right turn and cross products will do the trick I presume.

What's all that for: I'm working on subdivision curves and surfaces: Catmull-Clark, Loop, Doo-Sabin, Peters-Reif, Kobbelt, Chaikin, and others. Peters-Reif seemed too me the simplest algorithm but turned out to be one of the most complex, comparing to Catmull which I have coded in an hour or so.
I'm facing the CW/CCW and ordering problem, and the 'correctly orienting faces' problem. That's fun

@TT: didn't you encounter problems with == method when comparing Floats (in (a.x==b.x) for instance) ?
I've come across this:

plane1=Geom.fit_plane_to_points(pArray[0],pArray[1],pArray[2])
plane2=(Geom.fit_plane_to_points(pArray[0],pArray[1],pArray[3]))
return plane1==plane2

Assuming the 4 3dPoints are coplanar, this always return false. I had to code another method 'isEqual?' for the Float class that compares with a tolerance of 1 to the power of -10 to get correct results when comparing such arrays of floats.

thomthom · 21 Jan 2011, 11:13

@didier bur said:

@TT: didn't you encounter problems with == method when comparing Floats (in (a.x==b.x) for instance) ?

Good point. I haven't tested it that much, I just ported some pseudo code which I read in a book.
One should probably do a tolerance comparison.

[quote="Didier Bur"]


    plane1=Geom.fit_plane_to_points(pArray[0],pArray[1],pArray[2])
    plane2=(Geom.fit_plane_to_points(pArray[0],pArray[1],pArray[3]))
    return plane1==plane2

<span class="syntaxdefault"><br /></span><span class="syntaxcomment">#&nbsp;Check&nbsp;if&nbsp;all&nbsp;points&nbsp;in&nbsp;the&nbsp;array&nbsp;are&nbsp;on&nbsp;the&nbsp;same&nbsp;plane.<br />#<br />#&nbsp;@todo&nbsp;Ensure&nbsp;that&nbsp;the&nbsp;first&nbsp;three&nbsp;points&nbsp;are&nbsp;not&nbsp;co-linear.<br />#<br />#&nbsp;@param&nbsp;[Array<Geom;;Point3d>]&nbsp;points<br />#<br />#&nbsp;@since&nbsp;2.0.0<br /></span><span class="syntaxdefault">def&nbsp;self</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">planar_points</span><span class="syntaxkeyword">?(</span><span class="syntaxdefault">points</span><span class="syntaxkeyword">)<br />&nbsp;&nbsp;return&nbsp;</span><span class="syntaxdefault">false&nbsp;</span><span class="syntaxkeyword">if&nbsp;</span><span class="syntaxdefault">points</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">length&nbsp;</span><span class="syntaxkeyword"><&nbsp;</span><span class="syntaxdefault">3<br />&nbsp;&nbsp;v1&nbsp;</span><span class="syntaxkeyword">=&nbsp;</span><span class="syntaxdefault">points</span><span class="syntaxkeyword">[</span><span class="syntaxdefault">0</span><span class="syntaxkeyword">].</span><span class="syntaxdefault">vector_to</span><span class="syntaxkeyword">(</span><span class="syntaxdefault">points</span><span class="syntaxkeyword">[</span><span class="syntaxdefault">1</span><span class="syntaxkeyword">])<br />&nbsp;&nbsp;</span><span class="syntaxdefault">v2&nbsp;</span><span class="syntaxkeyword">=&nbsp;</span><span class="syntaxdefault">points</span><span class="syntaxkeyword">[</span><span class="syntaxdefault">0</span><span class="syntaxkeyword">].</span><span class="syntaxdefault">vector_to</span><span class="syntaxkeyword">(</span><span class="syntaxdefault">points</span><span class="syntaxkeyword">[</span><span class="syntaxdefault">2</span><span class="syntaxkeyword">])<br />&nbsp;&nbsp;return&nbsp;</span><span class="syntaxdefault">false&nbsp;unless&nbsp;v1</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">valid</span><span class="syntaxkeyword">?&nbsp;&&&nbsp;</span><span class="syntaxdefault">v2</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">valid</span><span class="syntaxkeyword">?<br />&nbsp;&nbsp;</span><span class="syntaxdefault">vector&nbsp;</span><span class="syntaxkeyword">=&nbsp;</span><span class="syntaxdefault">v1&nbsp;</span><span class="syntaxkeyword">*&nbsp;</span><span class="syntaxdefault">v2<br />&nbsp;&nbsp;plane&nbsp;</span><span class="syntaxkeyword">=&nbsp;[&nbsp;</span><span class="syntaxdefault">points</span><span class="syntaxkeyword">[</span><span class="syntaxdefault">0</span><span class="syntaxkeyword">],&nbsp;</span><span class="syntaxdefault">vector&nbsp;</span><span class="syntaxkeyword">]<br />&nbsp;&nbsp;</span><span class="syntaxdefault">3.upto</span><span class="syntaxkeyword">(</span><span class="syntaxdefault">points</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">length</span><span class="syntaxkeyword">-</span><span class="syntaxdefault">1</span><span class="syntaxkeyword">)&nbsp;{&nbsp;|</span><span class="syntaxdefault">i</span><span class="syntaxkeyword">|<br />&nbsp;&nbsp;&nbsp;&nbsp;return&nbsp;</span><span class="syntaxdefault">false&nbsp;unless&nbsp;points</span><span class="syntaxkeyword">[</span><span class="syntaxdefault">i</span><span class="syntaxkeyword">].</span><span class="syntaxdefault">on_plane</span><span class="syntaxkeyword">?(</span><span class="syntaxdefault">plane</span><span class="syntaxkeyword">)<br />&nbsp;&nbsp;}<br />&nbsp;&nbsp;return&nbsp;</span><span class="syntaxdefault">true<br />end<br /></span>

thomthom · 21 Jan 2011, 11:14

@didier bur said:

What's all that for: I'm working on subdivision curves and surfaces: Catmull-Clark, Loop, Doo-Sabin, Peters-Reif, Kobbelt, Chaikin, and others. Peters-Reif seemed too me the simplest algorithm but turned out to be one of the most complex, comparing to Catmull which I have coded in an hour or so.

I wish I could understand algorithm formulas. I always need to search up a pseudo code version, or an implementation in another language to work it out.

thomthom · 21 Jan 2011, 11:17

@didier bur said:

@TT: didn't you encounter problems with == method when comparing Floats (in (a.x==b.x) for instance) ?

Actually, Point3d.x, Point3d.y and Point3d.z return Length objects. They have built in tolerances - same as SU's.
You can easily get in trouble with unpredictable results if you're not 100% sure if you're dealing with Floats or Lengths.

Didier Bur · 21 Jan 2011, 11:28

@unknownuser said:

I wish I could understand algorithm formulas. I always need to search up a pseudo code version, or an implementation in another language to work it out.

Me too. When I speak of 'algorithm' I mean pseudo-code or (even better) natural language

Here is my snippet for coplanarity test of 4 points (can be extended to any number of points):

# Given a set of 4 3dPoints, return true if coplanar, or false
	def quadPointsCoplanar?(pArray)
		return (Geom.fit_plane_to_points(pArray[0],pArray[1],pArray[2]).isEqual?(Geom.fit_plane_to_points(pArray[0],pArray[1],pArray[3])))
	end

Where Array of floats comparison method must be 'aliased' with isEqual?

class Array
	def isEqual?(otherArray)
		self.each_index { |i| return false if self[i].isNotEqual(otherArray[i]) }
		return true
	end
end

and where == method for floats must be aliased with isNotEqual()

class Float
	def isNotEqual(value)
		delta = self - value
		return ((delta > Math;;EPSILON) or (delta < -Math;;EPSILON))
	end
end

and finaly requires a constant

Math.const_set("EPSILON", 1.0e-10)

thomthom · 21 Jan 2011, 11:44

I came across a very good site on the topic of floating precision: http://floating-point-gui.de/
Comparison with tolerance is a bit more complex in order to account for the various pit-falls.

The Floating-Point Guide - Comparison

Explanation of the various pitfalls in comparing floating-point numbers.

(floating-point-gui.de)

thomthom · 21 Jan 2011, 11:52

Ooops! Seems that I forgot to post all required methods:

<span class="syntaxdefault"><br />module TT_Lib<br /> module Geom2D<br />  <br />  </span><span class="syntaxcomment"># Graham's Scan<br /></span><span class="syntaxdefault">  </span><span class="syntaxcomment"># Computational Geometry Algorythms and Appliations (ISBN-10; 3540779736)<br /></span><span class="syntaxdefault">  </span><span class="syntaxcomment"># Chapter 1.1<br /></span><span class="syntaxdefault">  </span><span class="syntaxcomment">#<br /></span><span class="syntaxdefault">  </span><span class="syntaxcomment"># @param [Array<Geom;;Point3d>] points<br /></span><span class="syntaxdefault">  </span><span class="syntaxcomment"># @return [Array<Geom;;Point3d>]<br /></span><span class="syntaxdefault">  </span><span class="syntaxcomment"># @since 1.1.0<br /></span><span class="syntaxdefault">  def self</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">convex_hull</span><span class="syntaxkeyword">(</span><span class="syntaxdefault">points</span><span class="syntaxkeyword">)<br /></span><span class="syntaxdefault">    return nil if points</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">length </span><span class="syntaxkeyword"><</span><span class="syntaxdefault"> 2<br />    <br />    points </span><span class="syntaxkeyword">=</span><span class="syntaxdefault"> self</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">sort_points_by_x_y</span><span class="syntaxkeyword">(</span><span class="syntaxdefault">points</span><span class="syntaxkeyword">)<br /></span><span class="syntaxdefault">    <br />    l_upper </span><span class="syntaxkeyword">=</span><span class="syntaxdefault"> </span><span class="syntaxkeyword">[</span><span class="syntaxdefault"> points</span><span class="syntaxkeyword">[</span><span class="syntaxdefault">0</span><span class="syntaxkeyword">],</span><span class="syntaxdefault"> points</span><span class="syntaxkeyword">[</span><span class="syntaxdefault">1</span><span class="syntaxkeyword">]</span><span class="syntaxdefault"> </span><span class="syntaxkeyword">]<br /></span><span class="syntaxdefault">    2</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">upto</span><span class="syntaxkeyword">(</span><span class="syntaxdefault">points</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">length </span><span class="syntaxkeyword">-</span><span class="syntaxdefault"> 1</span><span class="syntaxkeyword">)</span><span class="syntaxdefault"> </span><span class="syntaxkeyword">{</span><span class="syntaxdefault"> </span><span class="syntaxkeyword">|</span><span class="syntaxdefault">i</span><span class="syntaxkeyword">|<br /></span><span class="syntaxdefault">      l_upper </span><span class="syntaxkeyword"><<</span><span class="syntaxdefault"> points</span><span class="syntaxkeyword">[</span><span class="syntaxdefault">i</span><span class="syntaxkeyword">]<br /></span><span class="syntaxdefault">      while l_upper</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">length </span><span class="syntaxkeyword">></span><span class="syntaxdefault"> 2 </span><span class="syntaxkeyword">&&</span><span class="syntaxdefault"> </span><span class="syntaxkeyword">!</span><span class="syntaxdefault">self</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">right_turn</span><span class="syntaxkeyword">?(</span><span class="syntaxdefault">l_upper</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">last</span><span class="syntaxkeyword">(</span><span class="syntaxdefault">3</span><span class="syntaxkeyword">))<br /></span><span class="syntaxdefault">        l_upper</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">delete_at</span><span class="syntaxkeyword">(-</span><span class="syntaxdefault">2</span><span class="syntaxkeyword">)<br /></span><span class="syntaxdefault">      end<br />    </span><span class="syntaxkeyword">}<br /></span><span class="syntaxdefault">    <br />    l_lower </span><span class="syntaxkeyword">=</span><span class="syntaxdefault"> </span><span class="syntaxkeyword">[</span><span class="syntaxdefault"> points</span><span class="syntaxkeyword">[-</span><span class="syntaxdefault">1</span><span class="syntaxkeyword">],</span><span class="syntaxdefault"> points</span><span class="syntaxkeyword">[-</span><span class="syntaxdefault">2</span><span class="syntaxkeyword">]</span><span class="syntaxdefault"> </span><span class="syntaxkeyword">]<br /></span><span class="syntaxdefault">    </span><span class="syntaxkeyword">(</span><span class="syntaxdefault">points</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">length </span><span class="syntaxkeyword">-</span><span class="syntaxdefault"> 3</span><span class="syntaxkeyword">).</span><span class="syntaxdefault">downto</span><span class="syntaxkeyword">(</span><span class="syntaxdefault">0</span><span class="syntaxkeyword">)</span><span class="syntaxdefault"> </span><span class="syntaxkeyword">{</span><span class="syntaxdefault"> </span><span class="syntaxkeyword">|</span><span class="syntaxdefault">i</span><span class="syntaxkeyword">|<br /></span><span class="syntaxdefault">      l_lower </span><span class="syntaxkeyword"><<</span><span class="syntaxdefault"> points</span><span class="syntaxkeyword">[</span><span class="syntaxdefault">i</span><span class="syntaxkeyword">]<br /></span><span class="syntaxdefault">      while l_lower</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">length </span><span class="syntaxkeyword">></span><span class="syntaxdefault"> 2 </span><span class="syntaxkeyword">&&</span><span class="syntaxdefault"> </span><span class="syntaxkeyword">!</span><span class="syntaxdefault">self</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">right_turn</span><span class="syntaxkeyword">?(</span><span class="syntaxdefault">l_lower</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">last</span><span class="syntaxkeyword">(</span><span class="syntaxdefault">3</span><span class="syntaxkeyword">))<br /></span><span class="syntaxdefault">        l_lower</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">delete_at</span><span class="syntaxkeyword">(-</span><span class="syntaxdefault">2</span><span class="syntaxkeyword">)<br /></span><span class="syntaxdefault">      end<br />    </span><span class="syntaxkeyword">}<br /></span><span class="syntaxdefault">    l_lower</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">delete_at</span><span class="syntaxkeyword">(</span><span class="syntaxdefault">0</span><span class="syntaxkeyword">)<br /></span><span class="syntaxdefault">    l_lower</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">delete_at</span><span class="syntaxkeyword">(-</span><span class="syntaxdefault">1</span><span class="syntaxkeyword">)<br /></span><span class="syntaxdefault">    <br />    return l_upper </span><span class="syntaxkeyword">+</span><span class="syntaxdefault"> l_lower<br />  end<br />  <br />  <br />  </span><span class="syntaxcomment"># @param [Array<Geom;;Point3d>] points<br /></span><span class="syntaxdefault">  </span><span class="syntaxcomment"># @return [Array<Geom;;Point3d>]<br /></span><span class="syntaxdefault">  </span><span class="syntaxcomment"># @since 1.1.0<br /></span><span class="syntaxdefault">  def self</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">sort_points_by_x</span><span class="syntaxkeyword">(</span><span class="syntaxdefault">points</span><span class="syntaxkeyword">)<br /></span><span class="syntaxdefault">    return points</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">sort </span><span class="syntaxkeyword">{</span><span class="syntaxdefault"> </span><span class="syntaxkeyword">|</span><span class="syntaxdefault">a</span><span class="syntaxkeyword">,</span><span class="syntaxdefault">b</span><span class="syntaxkeyword">|</span><span class="syntaxdefault"> a</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">x </span><span class="syntaxkeyword"><=></span><span class="syntaxdefault"> b</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">x </span><span class="syntaxkeyword">}<br /></span><span class="syntaxdefault">  end<br />  <br />  <br />  </span><span class="syntaxcomment"># Lexicographic order x, y<br /></span><span class="syntaxdefault">  </span><span class="syntaxcomment">#<br /></span><span class="syntaxdefault">  </span><span class="syntaxcomment"># @param [Array<Geom;;Point3d>] points<br /></span><span class="syntaxdefault">  </span><span class="syntaxcomment"># @return [Array<Geom;;Point3d>]<br /></span><span class="syntaxdefault">  </span><span class="syntaxcomment"># @since 1.1.0<br /></span><span class="syntaxdefault">  def self</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">sort_points_by_x_y</span><span class="syntaxkeyword">(</span><span class="syntaxdefault">points</span><span class="syntaxkeyword">)<br /></span><span class="syntaxdefault">    return points</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">sort </span><span class="syntaxkeyword">{</span><span class="syntaxdefault"> </span><span class="syntaxkeyword">|</span><span class="syntaxdefault">a</span><span class="syntaxkeyword">,</span><span class="syntaxdefault">b</span><span class="syntaxkeyword">|</span><span class="syntaxdefault"> </span><span class="syntaxkeyword">(</span><span class="syntaxdefault">a</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">x</span><span class="syntaxkeyword">==</span><span class="syntaxdefault">b</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">x</span><span class="syntaxkeyword">)</span><span class="syntaxdefault"> </span><span class="syntaxkeyword">?</span><span class="syntaxdefault"> a</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">y </span><span class="syntaxkeyword"><=></span><span class="syntaxdefault"> b</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">y </span><span class="syntaxkeyword">;</span><span class="syntaxdefault"> a</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">x </span><span class="syntaxkeyword"><=></span><span class="syntaxdefault"> b</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">x </span><span class="syntaxkeyword">}<br /></span><span class="syntaxdefault">  end<br />  <br />  <br />  </span><span class="syntaxcomment"># @param [Array<Geom;;Point3d>] points array of 3 +Geom;;Point3d+ objects.<br /></span><span class="syntaxdefault">  </span><span class="syntaxcomment"># @since 1.1.0 - Renamed from +.right_turn+ to +.right_turn?+ at version 1.2.0<br /></span><span class="syntaxdefault">  def self</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">right_turn</span><span class="syntaxkeyword">?(</span><span class="syntaxdefault">points</span><span class="syntaxkeyword">)<br /></span><span class="syntaxdefault">    p</span><span class="syntaxkeyword">,</span><span class="syntaxdefault"> q</span><span class="syntaxkeyword">,</span><span class="syntaxdefault"> r </span><span class="syntaxkeyword">=</span><span class="syntaxdefault"> points<br />    array </span><span class="syntaxkeyword">=</span><span class="syntaxdefault"> </span><span class="syntaxkeyword">[</span><span class="syntaxdefault"> 1</span><span class="syntaxkeyword">,</span><span class="syntaxdefault"> p</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">x</span><span class="syntaxkeyword">,</span><span class="syntaxdefault"> p</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">y</span><span class="syntaxkeyword">,<br /></span><span class="syntaxdefault">          1</span><span class="syntaxkeyword">,</span><span class="syntaxdefault"> q</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">x</span><span class="syntaxkeyword">,</span><span class="syntaxdefault"> q</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">y</span><span class="syntaxkeyword">,<br /></span><span class="syntaxdefault">          1</span><span class="syntaxkeyword">,</span><span class="syntaxdefault"> r</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">x</span><span class="syntaxkeyword">,</span><span class="syntaxdefault"> r</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">y </span><span class="syntaxkeyword">]<br /></span><span class="syntaxdefault">    return self</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">determinant_3x3</span><span class="syntaxkeyword">(array)</span><span class="syntaxdefault"> </span><span class="syntaxkeyword"><</span><span class="syntaxdefault"> 0.0<br />  end<br />  <br />  </span><span class="syntaxcomment"># @param [Array<Number>] array<br /></span><span class="syntaxdefault">  </span><span class="syntaxcomment"># @return [Number]<br /></span><span class="syntaxdefault">  </span><span class="syntaxcomment"># @since 1.1.0<br /></span><span class="syntaxdefault">  def self</span><span class="syntaxkeyword">.</span><span class="syntaxdefault">determinant_3x3</span><span class="syntaxkeyword">(array)<br /></span><span class="syntaxdefault">    a</span><span class="syntaxkeyword">,</span><span class="syntaxdefault">b</span><span class="syntaxkeyword">,</span><span class="syntaxdefault">c</span><span class="syntaxkeyword">,</span><span class="syntaxdefault">d</span><span class="syntaxkeyword">,</span><span class="syntaxdefault">e</span><span class="syntaxkeyword">,</span><span class="syntaxdefault">f</span><span class="syntaxkeyword">,</span><span class="syntaxdefault">g</span><span class="syntaxkeyword">,</span><span class="syntaxdefault">h</span><span class="syntaxkeyword">,</span><span class="syntaxdefault">i </span><span class="syntaxkeyword">=</span><span class="syntaxdefault"> array<br />    return a</span><span class="syntaxkeyword">*</span><span class="syntaxdefault">e</span><span class="syntaxkeyword">*</span><span class="syntaxdefault">i </span><span class="syntaxkeyword">-</span><span class="syntaxdefault"> a</span><span class="syntaxkeyword">*</span><span class="syntaxdefault">f</span><span class="syntaxkeyword">*</span><span class="syntaxdefault">h </span><span class="syntaxkeyword">+</span><span class="syntaxdefault"> b</span><span class="syntaxkeyword">*</span><span class="syntaxdefault">f</span><span class="syntaxkeyword">*</span><span class="syntaxdefault">g </span><span class="syntaxkeyword">-</span><span class="syntaxdefault"> b</span><span class="syntaxkeyword">*</span><span class="syntaxdefault">d</span><span class="syntaxkeyword">*</span><span class="syntaxdefault">i </span><span class="syntaxkeyword">+</span><span class="syntaxdefault"> c</span><span class="syntaxkeyword">*</span><span class="syntaxdefault">d</span><span class="syntaxkeyword">*</span><span class="syntaxdefault">h </span><span class="syntaxkeyword">-</span><span class="syntaxdefault"> c</span><span class="syntaxkeyword">*</span><span class="syntaxdefault">e</span><span class="syntaxkeyword">*</span><span class="syntaxdefault">g<br />  end<br />  <br /> end </span><span class="syntaxcomment"># module Geom2D<br /></span><span class="syntaxdefault">end </span><span class="syntaxcomment"># module TT_Lib<br /></span><span class="syntaxdefault"> </span>

fredo6 · 22 Jan 2011, 16:49

@didier bur said:

Here is my snippet for coplanarity test of 4 points (can be extended to any number of points):

# Given a set of 4 3dPoints, return true if coplanar, or false
> 	def quadPointsCoplanar?(pArray)
> 		return (Geom.fit_plane_to_points(pArray[0],pArray[1],pArray[2]).isEqual?(Geom.fit_plane_to_points(pArray[0],pArray[1],pArray[3])))
> 	end

I had a method for that in LibFredo6. Here is the code (slightly modified in order to return the plane or nil)


#Check if a set of points are coplanar
#  --> Return the plane or nil
def G6.points_coplanar?(pts)
	n = pts.length
	return nil if n < 3
	return Geom.fit_plane_to_points(pts) if n == 3
	
	n1 = [n / 2, 3].max
	n2 = [n - n1, 3].max
	plane1 = Geom.fit_plane_to_points *(pts[0..n1-1])
	plane2 = Geom.fit_plane_to_points *(pts[n-n2..n-1])
	normal1 = Geom;;Vector3d.new *plane1[0..2]
	normal2 = Geom;;Vector3d.new *plane2[0..2]
	(normal1.parallel?(normal2) && pts[0].on_plane?(plane2)) ? plane1 ; nil
end

Fredo

PS: I don't think it is a good idea to extend the Array and Float ruby classes. Better used a specific private method


def compare_float_arrays(array1, array2, epsilon)
   for i in 0..array1.length
       return false if (array1[i] - array2[i]).abs > epsilon
   end
   true
end

thomthom · 22 Jan 2011, 17:04

@unknownuser said:

PS: I don't think it is a good idea to extend the Array and Float ruby classes. Better used a specific private method

Ditto. Any Ruby and Sketchup base classes. There's been a few cases with clashes over methods like this.