MeshAlgo::distributePoints : Avoiding duplicate points

danieldresser-ie · danieldresser-ie · commit 3bb898d26d20 · 2023-05-04T17:03:49.000-07:00
diff --git a/src/IECoreScene/MeshAlgoDistributePoints.cpp b/src/IECoreScene/MeshAlgoDistributePoints.cpp
@@ -136,6 +136,101 @@ T triangleInterpolatedPrimVarValue(
 	}
 }
 
+class TriangleTester
+{
+public:
+	TriangleTester( const V2f (&points)[3] )
+	{
+		for( int i = 0; i < 3; i++ )
+		{
+			V2f a = points[ ( i + 1 ) % 3];
+			V2f b = points[ ( i + 2 ) % 3 ];
+
+			// Swap the vertices so that the edge is always handled in the same order, regardless of
+			// which triangle it belongs to.  This ensures consistency that if two triangles share an
+			// edge, any point will always be one side of the edge or the other
+			bool swap;
+			if( a.y != b.y )
+			{
+				swap = a.y > b.y;
+			}
+			else
+			{
+				swap = a.x < b.x;
+			}
+
+			if( swap )
+			{
+				V2f temp = a;
+				a = b;
+				b = temp;
+			}
+
+			// Store values that will make it easier to compute the signed area of a triangle formed
+			// by connected a point to this edge
+			m_edgeOrigins[i] = a;
+			m_edgeNormals[i] = V2f( b.y - a.y, a.x - b.x );
+
+			// Compute the signed edge of the triangle relative to this edge.
+			// Note that all signed areas in this class are actually stored as 2 times the area
+			// of the triangle, since the area of the extended parallelogram is what naturally
+			// falls out of the cross product, and all the areas are relative so it doesn't matter
+			V2f c = points[ i ];
+			float doubleSignedArea = ( c - m_edgeOrigins[i] ).dot( m_edgeNormals[i] );
+
+			// Store which side of the edge is inside the triangle
+			m_insideDir[i] = doubleSignedArea >= 0;
+
+			if( i == 0 )
+			{
+				// Having computed the signed area of the whole triangle to get the sign, we can store it
+				// so we have our divisor for when we divide the areas to get barycentric coordinates.
+				m_areaNormalize = 1.0f / fabsf( doubleSignedArea );
+			}
+		}
+	}
+
+	inline bool contains( const V2f &p, V3f &bary )
+	{
+		// Compute the signed areas formed by connecting this point to the 3 edges
+		float doubleSignedAreas[3];
+		for( int i = 0; i < 3; i++ )
+		{
+			doubleSignedAreas[i] = ( p - m_edgeOrigins[i] ).dot( m_edgeNormals[i] );
+		}
+
+		// To be inside the triangle, all the comparisons must match.  Note that
+		// m_insideDir stores which side we need to be on, but the comparison is
+		// always a >= comparison that is the same on either side of the edge.
+		// This ensures that a point near the edge will appear in exactly one
+		// of two triangles sharing the edge.
+		if(
+			( ( doubleSignedAreas[0] >= 0 ) != m_insideDir[0] ) ||
+			( ( doubleSignedAreas[1] >= 0 ) != m_insideDir[1] ) ||
+			( ( doubleSignedAreas[2] >= 0 ) != m_insideDir[2] )
+		)
+		{
+			return false;
+		}
+
+		// Compute 2 barycentric coordinates by using the ratios of the areas to the total area.
+		bary[0] = fabsf( doubleSignedAreas[0] ) * m_areaNormalize;
+		bary[1] = fabsf( doubleSignedAreas[1] ) * m_areaNormalize;
+
+		// The 3rd barycentric coordinate is determined by the first two.
+		bary[2] = 1.0f - bary[0] - bary[1];
+
+		return true;
+	}
+
+private:
+
+	float m_areaNormalize;
+	V2f m_edgeOrigins[3];
+	V2f m_edgeNormals[3];
+	bool m_insideDir[3];
+};
+
 void processInputs(
 	const MeshPrimitive *mesh,
 	const std::string &refPosition, const std::string &uvSet, const std::string &densityMask,
@@ -255,16 +350,16 @@ void distributePointsInTriangle(
 	std::vector< BaryAndFaceIdx >& results, const Canceller *canceller
 )
 {
-	Imath::V2f uv0, uv1, uv2;
-	triangleCornerPrimVarValues< Imath::V2f >( uvInterpolation, uvView, vertexIds, faceIdx, uv0, uv1, uv2 );
-	uv0 += offset;
-	uv1 += offset;
-	uv2 += offset;
+	Imath::V2f uvs[3];
+	triangleCornerPrimVarValues< Imath::V2f >( uvInterpolation, uvView, vertexIds, faceIdx, uvs[0], uvs[1], uvs[2] );
+	uvs[0] += offset;
+	uvs[1] += offset;
+	uvs[2] += offset;
 
 	Imath::Box2f uvBounds;
-	uvBounds.extendBy( uv0 );
-	uvBounds.extendBy( uv1 );
-	uvBounds.extendBy( uv2 );
+	uvBounds.extendBy( uvs[0] );
+	uvBounds.extendBy( uvs[1] );
+	uvBounds.extendBy( uvs[2] );
 
 	const float maxCandidatePoints = 1e9;
 	const float approxCandidatePoints = uvBounds.size().x * uvBounds.size().y * textureDensity;
@@ -307,6 +402,7 @@ void distributePointsInTriangle(
 		cornerDensities[2] *= invMaxDensity;
 	}
 
+	TriangleTester triTester( uvs );
 	PointDistribution::defaultInstance()(
 		uvBounds, finalDensity,
 		[&]( const Imath::V2f pos, float densityThreshold )
@@ -318,7 +414,7 @@ void distributePointsInTriangle(
 			}
 
 			Imath::V3f bary;
-			if( triangleContainsPoint( uv0, uv1, uv2, pos, bary ) )
+			if( triTester.contains( pos, bary ) )
 			{
 				if( hasDensityVar )
 				{
diff --git a/test/IECoreScene/MeshAlgoDistributePointsTest.py b/test/IECoreScene/MeshAlgoDistributePointsTest.py
@@ -36,6 +36,7 @@
 import re
 import sys
 import unittest
+import random
 import imath
 import time
 import threading
@@ -45,7 +46,46 @@
 
 class MeshAlgoDistributePointsTest( unittest.TestCase ) :
 
-	def pointTest( self, mesh, points, density, error=0.05 ) :
+	def checkForDuplicates( self, positions ) :
+		duplicatePositions = []
+		uniquePositions = set()
+		for p in positions:
+			prevNumUnique = len( uniquePositions )
+			uniquePositions.add( tuple( p ) )
+			if not len( uniquePositions ) > prevNumUnique:
+				duplicatePositions.append( p )
+
+		self.assertEqual( duplicatePositions, [] )
+
+	def checkVectorListsAlmostEqualUnordered( self, aIn, bIn, tolerance ) :
+		l = len( aIn )
+		self.assertEqual( len( bIn ), l )
+		a = IECore.V3fVectorData( sorted( aIn, key = lambda v : v[0] ) )
+		b = IECore.V3fVectorData( sorted( bIn, key = lambda v : v[0] ) )
+		aMatched = IECore.BoolVectorData( l )
+
+		aStartIndex = 0
+
+		for i in b:
+			while aStartIndex < l and a[aStartIndex][0] < i[0] - tolerance:
+				aStartIndex += 1
+
+			matched = False
+			aIndex = aStartIndex
+			while aIndex < l and a[aIndex][0] < i[0] + tolerance:
+				if not aMatched[aIndex]:
+					if i.equalWithAbsError( a[aIndex], tolerance ):
+						aMatched[aIndex] = True
+						matched = True
+						continue
+				aIndex += 1
+
+			if not matched:
+				raise AssertionError( "No match for point ", i )
+
+		self.assertTrue( all( aMatched ) )
+
+	def pointTest( self, mesh, points, density, error=0.05, checkDuplicates = True ) :
 
 		self.assertTrue( "P" in points )
 		self.assertEqual( points.numPoints, points['P'].data.size() )
@@ -58,6 +98,9 @@ def pointTest( self, mesh, points, density, error=0.05 ) :
 		pointsPerFace = [ 0 ] * mesh.verticesPerFace.size()
 		positions = points["P"].data
 
+		if checkDuplicates:
+			self.checkForDuplicates( positions )
+
 		## test that the points are on the mesh
 		for p in positions :
 			self.assertAlmostEqual( meshEvaluator.signedDistance( p, result ), 0.0, 6 )
@@ -82,6 +125,30 @@ def testSimple( self ) :
 		p = IECoreScene.MeshAlgo.distributePoints( mesh = m, density = 100 )
 		self.pointTest( m, p, 100 )
 
+	def testActualDistribution( self ) :
+		# Most of these tests target the properties of the generated points, rather than the specific points
+		# chosen. This is good for the generality of the tests, and will come in handy if we ever switch
+		# the implementation.  But it's probably worth having one test that actually asserts the exact
+		# distribution we are currently using.
+
+		# This also checks the handling of points that lie on edges - we accept <0,0> as lying within the
+		# rect from 0 -> 1, but <1,0>, <0,1> and <1,1> are considered outside ( this ensures that if adjacent
+		# polygons are added, the point will be included exactly once )
+
+		m = IECoreScene.MeshPrimitive.createPlane( imath.Box2f( imath.V2f( 0 ), imath.V2f( 1 ) ), imath.V2i( 1 ) )
+		points = IECoreScene.MeshAlgo.distributePoints( mesh = m, density = 5 )["P"].data
+		refPoints = [
+			imath.V3f( 0 ),
+			imath.V3f( 0.612342954, 0.365361691, 0 ),
+			imath.V3f( 0.309612036, 0.182079479, 0 ),
+			imath.V3f( 0.0924436674, 0.497548342, 0 ),
+			imath.V3f( 0.512147605, 0.871189654, 0 )
+		]
+		self.assertEqual( len( points ), len( refPoints ) )
+		for a, b in zip( points, refPoints ):
+			with self.subTest( a=a, b=b ):
+				self.assertTrue( a.equalWithRelError( b, 0.0000001 ) )
+
 	def testRaisesExceptionIfInvalidUVs( self ) :
 
 		m = IECore.Reader.create( os.path.join( "test", "IECore", "data", "cobFiles", "pCubeShape1.cob" ) ).read()
@@ -129,7 +196,10 @@ def testDensityMaskPrimVar( self ) :
 
 		m['density'] = IECoreScene.PrimitiveVariable( IECoreScene.PrimitiveVariable.Interpolation.Constant, IECore.FloatData( 2 ) )
 		p = IECoreScene.MeshAlgo.distributePoints( mesh = m, density = 1000 )
-		self.pointTest( m, p, 1000, error=0.1 )
+		# We happen to hit an edge here which has two different UV space locations, so we don't avoid creating
+		# duplicates - that's probably OK, we just want to make sure that we don't produce the same UV point
+		# multiple times
+		self.pointTest( m, p, 1000, error=0.1, checkDuplicates = False )
 
 		m['density'] = IECoreScene.PrimitiveVariable( IECoreScene.PrimitiveVariable.Interpolation.Constant, IECore.FloatData( -0.001 ) )
 		p = IECoreScene.MeshAlgo.distributePoints( mesh = m, density = 1000 )
@@ -167,8 +237,8 @@ def testRefPosition( self ):
 		m = IECoreScene.MeshPrimitive.createPlane( imath.Box2f( imath.V2f( -1 ), imath.V2f( 1 ) ), imath.V2i( 1 ) )
 		m["altP"] = IECoreScene.PrimitiveVariable( IECoreScene.PrimitiveVariable.Interpolation.Vertex, IECore.V3fVectorData( [ i * 2 + imath.V3f( 1, 7, 10 ) for i in m["P"].data ] ) )
 
-		p1 = IECoreScene.MeshAlgo.distributePoints( mesh = m, density = 40 )
-		p2 = IECoreScene.MeshAlgo.distributePoints( mesh = m, density = 10, refPosition = "altP" )
+		p1 = IECoreScene.MeshAlgo.distributePoints( mesh = m, density = 400 )
+		p2 = IECoreScene.MeshAlgo.distributePoints( mesh = m, density = 100, refPosition = "altP" )
 
 		# Using a reference position with a scale affects the density ( which is compensated by changing the
 		# density argument, but otherwise does not affect anything )
@@ -179,7 +249,7 @@ def testRefPosition( self ):
 		# across a much larger area, but the point counts are kept constant
 		m["P"] = IECoreScene.PrimitiveVariable( IECoreScene.PrimitiveVariable.Interpolation.Vertex, IECore.V3fVectorData( [ i * 100 for i in m["P"].data ] ) )
 
-		p3 = IECoreScene.MeshAlgo.distributePoints( mesh = m, density = 10, refPosition = "altP" )
+		p3 = IECoreScene.MeshAlgo.distributePoints( mesh = m, density = 100, refPosition = "altP" )
 		self.assertLess( p3.bound().min()[0], -100 )
 		self.assertGreater( p3.bound().max()[0], 100 )
 		self.assertEqual( p3.numPoints, p2.numPoints )
@@ -227,6 +297,50 @@ def testDensityRange( self ) :
 		m = IECore.Reader.create( os.path.join( "test", "IECore", "data", "cobFiles", "pCubeShape1.cob" ) ).read()
 		self.assertRaises( RuntimeError, IECoreScene.MeshAlgo.distributePoints, m, -1.0 )
 
+	def testEdgeOverlaps( self ):
+		# As long as we don't change the edges or stretch the UVs, we should be able to add more subdivisions,
+		# or more vertices around within the mesh, and still get identical points out ( down to floating point
+		# error ). This test is run with enough points and polygons to have a good chance of triggering
+		# special case interactions of points lying on or near edges and vertices. We want to make sure that
+		# every point gets output exactly once, not being doubled up when it's on a shared edge between two
+		# polygons, or being missed.
+
+		dens = 10000
+		m = IECoreScene.MeshPrimitive.createPlane( imath.Box2f( imath.V2f( 0 ), imath.V2f( 1 ) ), imath.V2i( 1 ) )
+		refPoints = IECoreScene.MeshAlgo.distributePoints( mesh = m, density = dens )["P"].data
+
+		div = 256
+		m = IECoreScene.MeshPrimitive.createPlane( imath.Box2f( imath.V2f( 0 ), imath.V2f( 1 ) ), imath.V2i( div ) )
+		pointsFromHighDensity = IECoreScene.MeshAlgo.distributePoints( mesh = m, density = dens )["P"].data
+
+		m = IECoreScene.MeshPrimitive.createPlane( imath.Box2f( imath.V2f( 0 ), imath.V2f( 1 ) ), imath.V2i( div ) )
+		jitterUVs = m["uv"].data.copy()
+		jitterPs = m["P"].data.copy()
+
+		jitter = 0.25 / div
+		random.seed( 42 )
+		for y in range( 1, div ):
+			for x in range( 1, div ):
+				i = y * ( div + 1 ) + x
+				jitterUVs[ i ] = jitterUVs[ i ] + imath.V2f(
+					random.uniform( -jitter, jitter ), random.uniform( -jitter, jitter )
+				)
+				jitterPs[ i ] = imath.V3f( jitterUVs[ i ][0], jitterUVs[ i ][1], 0 )
+		m["uv"] = IECoreScene.PrimitiveVariable(
+			IECoreScene.PrimitiveVariable.Interpolation.FaceVarying, jitterUVs, m["uv"].indices
+		)
+		m["P"] = IECoreScene.PrimitiveVariable(
+			IECoreScene.PrimitiveVariable.Interpolation.Vertex, jitterPs
+		)
+		pointsFromJittered = IECoreScene.MeshAlgo.distributePoints( mesh = m, density = dens )["P"].data
+
+		self.checkForDuplicates( refPoints )
+		self.checkForDuplicates( pointsFromHighDensity )
+		self.checkForDuplicates( pointsFromJittered )
+
+		self.checkVectorListsAlmostEqualUnordered( pointsFromHighDensity, refPoints, 0.0000002 )
+		self.checkVectorListsAlmostEqualUnordered( pointsFromJittered, refPoints, 0.0000002 )
+
 	def testVertexUVs( self ) :
 
 		m = IECoreScene.MeshPrimitive.createPlane( imath.Box2f( imath.V2f( -1 ), imath.V2f( 1 ) ), imath.V2i( 4 ) )
@@ -259,7 +373,7 @@ def testPrimitiveVariables( self ):
 		self.assertEqual( p.keys(), ['N', 'P', 'cA', 'cB', 'fvA', 'fvB', 'uA', 'uB', 'uv', 'vA', 'vB', 'vC'] )
 		self.assertEqual( p['cA'], IECoreScene.PrimitiveVariable( IECoreScene.PrimitiveVariable.Interpolation.Constant, IECore.FloatData( 42 ) ) )
 		self.assertEqual( p['cB'], IECoreScene.PrimitiveVariable( IECoreScene.PrimitiveVariable.Interpolation.Constant, IECore.M44fData( imath.M44f( 7 ) ) ) )
-		self.assertEqual( p.numPoints, 15 )
+		self.assertEqual( p.numPoints, 10 )
 		for i in range( p.numPoints ):
 			v = { k : p[k].data[i] for k in p.keys() if k[0] != 'c' }
 			self.assertAlmostEqual( v['uv'][0], v['P'][0] + 0.5, places = 6 )
@@ -288,7 +402,7 @@ def testPrimitiveVariables( self ):
 		self.assertEqual( p.keys(), ['N', 'P', 'cA', 'cB', 'fvA', 'fvB', 'uA', 'uB', 'uv', 'vA', 'vB', 'vC'] )
 
 		# Check that the overrides applied correctly
-		self.assertEqual( p.numPoints, 410 )
+		self.assertEqual( p.numPoints, 406 )
 
 		# Test variable types that can't be interpolated
 		m['invalid1'] = IECoreScene.PrimitiveVariable( IECoreScene.PrimitiveVariable.Interpolation.FaceVarying, IECore.StringVectorData( [ "foo" ] * 16 ) )
