evenly distributed random points on a mesh

Mesh2Points.h

//
//  Mesh2Points.h
//
//  Created by lars berg on 8/1/14.
//

#pragma once

#include "ofMain.h"

namespace Mesh2Points
{
	
	static float areaOfTriangle(ofVec3f p0, ofVec3f p1, ofVec3f p2)
	{
		return (p2 - p1).cross(p0 - p1).length() * .5;
	}
	
	template <class T>
	static T pointInTriangle(T t0, T t1, T t2 , float u, float v)
	{
		if( ( u + v ) > 1 )
		{
			u = 1. - u;
			v = 1. - v;
		}
		
		return t0 * u + t1 * v + t2 * (1 - u - v);
	}
	
	class Face3
	{
	public:
		Face3(ofVec3f a, ofVec3f b, ofVec3f c) :
		a(a), b(b), c(c), area(areaOfTriangle(a,b,c))
		{};
		
		static bool compareAreas(Face3 a, Face3 b )
		{
			return a.area < b.area;
		}
		
		ofVec3f a, b, c;
		float area, normalizedArea, u;
	};
	
	static vector<Face3> getSortedAndWeightedMeshFaces(ofMesh& mesh, bool normalizeFaceAreas = true)
	{
		vector<Face3> faces;
		float totalArea = 0;
		
		ofVec3f minBound, maxBound;
		
		//center the mesh
		auto& v = mesh.getVertices();
		minBound = maxBound = v[0];
		for(auto& vert: v)
		{
			for(int i=0; i<3; i++)
			{
				minBound[i] = fmin(minBound[i], vert[i]);
				maxBound[i] = fmax(maxBound[i], vert[i]);
			}
		}
		
		ofVec3f centroid = (minBound + maxBound) * .5;
		for(auto &vert: v)	vert -= centroid;
		
		//create our faces
		if(mesh.getIndices().size() == 0)
		{
			for(int i=0; i<v.size(); i+=3)
			{
				Face3 f(v[i], v[i+1], v[i+2]);
				totalArea += f.area;
				faces.push_back(f);
			}
		}
		else
		{
			auto& indices = mesh.getIndices();
			for(int i=0; i<indices.size(); i+=3)
			{
				Face3 f(v[indices[i]], v[indices[i+1]], v[indices[i+2]]);
				totalArea += f.area;
				faces.push_back(f);
			}
		}
		
		//sort the faces by are
		sort(faces.begin(), faces.end(), Face3::compareAreas);
		
		//get the normalized face areas and a 'u' sampling value
		float u = 0;
		for(auto &f: faces)
		{
			f.normalizedArea = f.area / totalArea;
			u += f.normalizedArea;
			f.u = u;
		}
		
		return faces;
	}
	
	static vector<ofVec3f> getRandomSurfacePoints(ofMesh& mesh, int numPoints = 100000, float scale = 1.)
	{
		vector<ofVec3f> points(numPoints);
		
		//get our faces
		auto faces = getSortedAndWeightedMeshFaces(mesh, true);
		
		//find a bunch of random points on the surface
		for(auto &p: points)
		{
			float lowThreshold = ofRandom(0., 1.);
			for(auto& f: faces)
			{
				if(lowThreshold < f.u)
				{
					p = pointInTriangle(f.a, f.b, f.c, ofRandom(0.,1.), ofRandom(0.,1.)) * scale;
					break;
				}
			}
		}
		
		return points;
	}
}