Alternative to angle_between?

[Didier Bur]

Hi,
If this can help:

def angle2Pi(p1,p2)
	twoPI = 2.0*Math;;PI
	return (Math.atan2(p2.y-p1.y, p2.x-p1.x)+twoPI).modulo(twoPI)
end

This returns the angle of a vector that starts at point p1 and ends at point p2, no matter wether p1 or p2 is the starting vertex. this is recommanded you use this only in 2D of course.
Regards,

[thomthom]

@unknownuser said:

@unknownuser said:

You have to choose a Z direction (perpendicular to the plane of the two vectors, in order ot orient your angles.
Otherwise it is normal that you get the same value (just imagine you look at your picture from the back!).

Then you make a difference by the sign of (vec1 * vec2) % vecZ

Fredo

Just to make clearer what I said in an earlier post. Referring to the picture, you have 3 vectors: red, blue and black.

If you wish to measure the 'oriented' angles between 'blue and black' and 'blue and red', you must decide on a reference rotation direction (in plane geometry, we would count in the trigonometric sense, or anti-clockwise).

The reference direction is what you want, but is perpendicular to the other vectors. So, for instance, if you are fine with the angle returned by vec_black.angle_between(vec_blue), then the reference direction would be

> vec_rotation_ref = vec_black * vec_blue
> 

Now, for the red vector, you compute vec_black * vec_red
and then compare its orientation with your reference rotation vector, by making the scalar product (which is negative if the vectors have opposite directions, and positive otherwise).

> angle_red = vec_black.angle_between(vec_red)
> if vec_black * vec_red) % vec_rotation_ref < 0
>    angle_red = 2 * Math;PI - angle_red
> end
> 

Note that if you work in the horizontal plane, the you can take Z_AXIS as your reference rotation vector (and then Bill'sformula works).

Fredo

Ok, I've been trying to get this to work no-planar. But I'm having problems with the reference vector - because I don't have both the vectors to compare against at that point of the code.
I only have a base vector - that one that I compare against.

I've been trying to take the cross of the base vector (black) and the vector I compare against (red/blue), but it seem to be 50/50 chance that the cross vector points in opposite directions. Not sure how to orient it predictably....

[bentleykfrog]

Sorry for digging up another old topic, here's my not-so-rough guess at a long way around.

angle1 = black_vector.angle_between red_vector
vector_cross = black_vector.cross red_vector
vector_90 = black_vector.cross vector_cross
angle2 = vector_90.angle_between red_vector
angle1 = (Math;;PI-angle1) + Math;;PI if angle2 >= (Math;;PI/2)

[bentleykfrog]

Also, as a bit of a bodge to fix the 50/50 direction problem

vector_90 = Geom;;Vector3d.new(vector_90[0].abs,vector_90[1].abs,vector_90[2].abs)

Its probably not efficient, or even tested

[bentleykfrog]

I've been working on this problem quite alot lately, and I've realised the code I posted is wildly incorrect, as you say ThomThom:
@thomthom said:

it seem to be 50/50 chance that the cross vector points in opposite directions.

so we should probably compare the cross vector against something relatively static, like one of the model axes?

@unknownuser said:

Note that if you work in the horizontal plane, the you can take Z_AXIS as your reference rotation vector (and then Bill'sformula works).

could this work if we reverse the cross vector based on its angle to the Z_AXIS (ie. greater than 90 degrees then reverse cross vector).

[bentleykfrog]

@bentleykfrog said:

I've been working on this problem quite alot lately

I've finally got my 2d manifold detection working (albeit quite slowly as I have to calculate an infinite line [planar normal of an edge] intersecting a complex polygon [face] for each edge to determine the side of the edge its associated face is on). Anyway I thought I'd post the code here for some feedback and optimization advice. Here's the code snippet: feel free to chop/change/reuse if you like.

		def intersect_line_line_segment_2d(line,segment)
			x1	= line[0].x
			y1	= line[0].y
			x2	= line[1].x
			y2	= line[1].y
			x3	= segment[0].x
			y3	= segment[0].y
			x4	= segment[1].x
			y4	= segment[1].y
			
			u_line_n = ((x4 - x3)*(y1 - y3)) - ((y4 - y3)*(x1 - x3))	#line equation numerator
			denom	 = ((y4 - y3)*(x2 - x1)) - ((x4 - x3)*(y2 - y1))	#line  & segment equation denominator
			
			u_segment_n	= ((x2 - x1)*(y1 - y3)) - ((y2 - y1)*(x1 - x3))	#segmnet equation numerator
			
			if denom == 0
				#line and segment are parallel
				return false
			end #if
			
			u_segment = u_segment_n.quo(denom)
			
			if u_segment >= 0 || u_segment <=1
				#return [x1 + (u_segment*(x2-x1)),y1 + (u_segment*(y2-y1)),1]	#use to return the point of intersection
				u_line = u_line_n.quo(denom)									#use to return a relative distance value for ordering
				return u_line													#intersections on a complex polygon
			else
				#intersection is outside the line segment
				return false
			end
		end #def
		
		def find_manifold_face(origin_face_id,fold_edge_id,possible_face_ids,check_orientation=false,correct_orientation=false)		
			origin_face		= @@face_objects[origin_face_id]
			origin_edges	= @@faces_to_edges[origin_face_id]
			manifold_face	= false
			
			#	PT1; ORIGIN FACE->FOLD EDGE ORIENTATION
			#	we need to find on which side of the edge the face is, and then
			#	determine the direction of rotation. This is more difficult than
			#	it appears. The most error-free way to do it is to create a line
			#	that represents the planar normal for the edge (ie perpendicular
			#	to the fold edge and on the plane of the origin face) and find the
			#	line segments that intersect with the face using an Intersect Segment
			#	with Polygon algorithm. [see; http://softsurfer.com/Archive/algorithm_0111/algorithm_0111.htm]
			#	A COUPLE OF NOTES FIRST;
			# 1;We only need a 2d coordinate system for this so we can use the UVHelper
			#	to generate this and hopefully make the code more efficient
			# 2;We only need to find the first intersected line segment that shares
			#	one of its points with the fold edge
			fold_verts 		= @@edges_to_verts[fold_edge_id]					#construct points on the fold edge that will help with
			fold_arr		= fold_verts.values									#orientation; pt1; start, pt2; middle, pt3; end
			fold_pt1		= fold_arr[0].position
			fold_pt3		= fold_arr[1].position
			fold_pt2		= fold_pt1 + Geom;;Vector3d.new((fold_pt3.x-fold_pt1.x).quo(2),(fold_pt3.y-fold_pt1.y).quo(2),(fold_pt3.z-fold_pt1.z).quo(2))
			fold_vec		= fold_pt1.vector_to fold_pt3
			trans 			= Geom;;Transformation.rotation fold_pt2, origin_face.normal, -90.degrees
			o_vec1 			= fold_vec.transform trans
			
			origin_uvhelp 	= origin_face.get_UVHelper true, false, @@tw		#UVHelper can give us 2d coordinates for a 3d face ;)
			fold_uv_pt1		= origin_uvhelp.get_front_UVQ(fold_pt1)				#construct the corresponding uv points that are uv
			fold_uv_pt2		= origin_uvhelp.get_front_UVQ(fold_pt2)				#representations of the 3d coordinates
			fold_uv_pt3		= origin_uvhelp.get_front_UVQ(fold_pt3)
			fold_uv_vec		= fold_uv_pt1.vector_to fold_uv_pt3
			
			trans_uv_vec	= Geom;;Vector3d.new(fold_uv_vec[1],-fold_uv_vec[0],1)		#construct a fold edge normal line to
			trans_uv_pt		= fold_uv_pt2 + trans_uv_vec								#calculate the line segments of intersection
			trans_uv_line	= [fold_uv_pt2, trans_uv_pt]
			
			intersection_array = Array[0]
			origin_edges.each {|edgeID|
				next if fold_edge_id == edgeID
				edge_uv_pt1	= origin_uvhelp.get_front_UVQ(@@edge_objects[edgeID].start.position)
				edge_uv_pt2 = origin_uvhelp.get_front_UVQ(@@edge_objects[edgeID].end.position)
				edge_uv_segment = [edge_uv_pt1, edge_uv_pt2]
				
				intersection_u = self.intersect_line_line_segment_2d(trans_uv_line,edge_uv_segment)
				next if !intersection_u
				intersection_array << intersection_u
			}
			if intersection_array.length <= 1
				msg = "Critical Error; Normalize Toolkit encountered an illegal face (a face with only two vertices).\nYour model probably has errors?"
				msg += "To check, select the menu item 'Window' -> Model Info, then select 'Statictics' and click on 'Fix Problems'."
				self.alert_error(msg,true)
				Sketchup.send_action CMD_SELECT
				return false			
			end #if
			
			intersection_array.sort!									#ok lets find if the vector is pointing the right way
			reverse = (intersection_array.index(0) % 2) ? 1 ; -1
			o_vec1.reverse! if reverse == -1
			o_pt1 = fold_pt2 + o_vec1
			o_vec2 	= origin_face.normal								#make the origin face's normal an orientation vector
																		#we need two vectors at right angles to get around the
																		#angle_between not greater than 180 degrees problem.
			smallest_angle 	= 360.degrees
			possible_face_ids.each {|faceID|
				next if !@@faces_to_verts.has_key?(faceID)
				next if faceID == origin_face_id
				
				#logically, the p_vec1 should be found by rotating the fold_vec in the opposite direction
				#as o_vec1 was. The only issue here is that if the normal is reversed, we wont get a correctly
				#facing vector. Fortunately, we can determine this by comparing the clockwise angle from o_vec1
				#to p_vec1 and from o_vec1 to @@face_objects[faceID].normal.
				#If o_vec1 to p_vec1 is smaller than o_vec1 to @@face_objects[faceID].normal, we've got a reversed
				#normal!! cool
				trans 			= Geom;;Transformation.rotation fold_pt2, @@face_objects[faceID].normal, reverse * 90.degrees
				p_vec1 			= fold_vec.transform trans
				
				t_ang1a			= o_vec1.angle_between @@face_objects[faceID].normal
				t_ang1b			= o_vec2.angle_between @@face_objects[faceID].normal
				t_ang2a			= o_vec1.angle_between p_vec1
				t_ang2b			= o_vec2.angle_between p_vec1
				
				t_ang1 = (t_ang1b <= 90.degrees) ? t_ang1a ; 180.degrees + (180.degrees - t_ang1a)
				t_ang2 = (t_ang2b <= 90.degrees) ? t_ang2a ; 180.degrees + (180.degrees - t_ang2a)
				
				#another problem here is that the difference between t_ang1 & t_ang2 should be roughly plus or minus 90 degrees
				#if its greater, then the possible face is on such an obtuse angle and reversed that its normal is closer to
				#o_vec1 than p_vec1 is
				p_vec1.reverse! if ((t_ang1 > t_ang2) || ((t_ang2 - t_ang1) > 180.degrees))
				
				o_ang1 = o_vec1.angle_between p_vec1
				o_ang2 = o_vec2.angle_between p_vec1
				
				manifold_angle = (o_ang2 <= 90.degrees) ? o_ang1 ; 180.degrees + (180.degrees - o_ang1)

				if manifold_angle < smallest_angle
					smallest_angle = manifold_angle
					manifold_face = @@face_objects[faceID]
				end
			}
			if !manifold_face
				return false
			end
			
			if check_orientation
				direction = self.get_rotation_direction(fold_vec,fold_pt2,origin_face.normal,o_pt1)				#work out whether the vector is pointing in the right direction
				trans 	= Geom;;Transformation.rotation fold_pt2, fold_vec, (smallest_angle * direction)		#rotate the origin normal so its aligned with
				m_normal 	= o_vec2.transform trans															#the manifold's normal. This should give a vector
				m_ang1		= m_normal.angle_between manifold_face.normal										#that angles 180 degrees from the manifold's normal
				oriented = (m_ang1 > 90.degrees) ? true ; false													#if its less than 90 degrees, reverse the face.
				return Hash[
					"manifold_face" => manifold_face,
					"oriented"	=> oriented
				]
			end #if
			
			return manifold_face
		end #def

[bentleykfrog]

just for context, attached is the wip with a few modifications to the previous post to reduce the load significantly if we're dealing with convex polygons.

---

normalize toolkit work in progress

[mac1]

See if any ideas here help. Tried to find arctan2 for three 3d. More search probably required
http://www.euclideanspace.com/maths/algebra/vectors/angleBetween/index.htm

[Dan Rathbun]

atan2() is one of the methods in the standard Ruby Math module.

Somewhere TIG posted an example method of using it with Sketchup API units.

[bentleykfrog]

wow, major oversight! I should probably be using edge.reversed_in? to reliably find the edge normal direction relative to the face.

Just to double-check my logic: on outer loops edge.reversed_in? will return true if the start and end vertex of the edge is running clockwise (if your looking against the normal of the face [the negative face normal vector is the vector for determining clockwise or anti-clockwise]). On inner loops (face.loops - face.outer_loop) edge.reversed_in will return true if the start and end vertex of the edge is running counter-clockwise?

[mac1]

Once one determines the angle value the inverse trig functions calculate the principal value. Therefore to determine the actual angle logic must be use either in the method or in the program its self to determine the correct quadrant. Give one can get the cos from the dot product, sine from the cross and even tangent from tan( theta) = |V1xV2|/ (V1 dot V2) then it is possible to establish the quadrant??

[thomthom]

@bentleykfrog said:

wow, major oversight! I should probably be using edge.reversed_in? to reliably find the edge normal direction relative to the face.

Just to double-check my logic: on outer loops edge.reversed_in? will return true if the start and end vertex of the edge is running clockwise (if your looking against the normal of the face [the negative face normal vector is the vector for determining clockwise or anti-clockwise]). On inner loops (face.loops - face.outer_loop) edge.reversed_in will return true if the start and end vertex of the edge is running counter-clockwise?

Aye. Loop and EdgeUse objects can be very useful.

[bentleykfrog]

@bentleykfrog said:

Just to double-check my logic: on outer loops edge.reversed_in? will return true if the start and end vertex of the edge is running clockwise (if your looking against the normal of the face [the negative face normal vector is the vector for determining clockwise or anti-clockwise]). On inner loops (face.loops - face.outer_loop) edge.reversed_in will return true if the start and end vertex of the edge is running counter-clockwise?

I'm probably getting a bit off topic here, so to sum up, when you're looking against the normal of the face edge.reversed_in? returns true when the face is on the left hand side of the edge, and false when on the right hand side of the face (irrespective of inner and outer loops, so disregard my previous assumption).

@mac1 said:

Once one determines the angle value the inverse trig functions calculate the principal value. Therefore to determine the actual angle logic must be use either in the method or in the program its self to determine the correct quadrant. Give one can get the cos from the dot product, sine from the cross and even tangent from tan( theta) = |V1xV2|/ (V1 dot V2) then it is possible to establish the quadrant??

I agree, the problem I had in Sketchup (I think) is determining a consistent cross and dot product vectors. Once you've established consistent methods for determining these, you should be good to go.

In my case I'm trying to define a 2d manifold. What's important in a 2d manifold is the angle between two faces (complex polygons) at a common edge. Given this angle, we can establish 2d manifold errors (like an edge sharing 3 or more faces).

To establish this angle we need three vectors: 'the edge vector normal on the first face', 'the face normal' & the 'the edge vector normal on the second face'. edge.reversed_in? seems to be a very consistent way to get the right edge vector normals. Also, to eliminate errors with faces that are reversed, we need to define a 'dominant' face. This dominant face is the origin of the 2d manifold, so all other connected faces should be oriented in the same direction as the dominant face (ie. all looking into the 2d manifold, or all looking out of the 2d manifold), if not, their normal vector is reversed (not the face itself) so the angle between two faces remains consistent. In my code, the dominant face is the face the user hovers over with their cursor.

@thomthom said:

Aye. Loop and EdgeUse objects can be very useful.

---

Normalize Toolkit Work In Progress 2

[mac1]

Not ruby programmer so leave that up to you, but for 3d vectors and using the general approach the angle is the shortest great circle path between the two. The angle is atan2( Norm( cross(a,b)), dot( a,b)). This means the angle is 0 to pi()
FYI http://www.euclideanspace.com/maths/algebra/vectors/angleBetween/index.htm and http://www.mathworks.com/matlabcentral/newsreader/view_thread/151925