/*
Szymon Rusinkiewicz
Princeton University

pointsample.cc
Return n randomly-distributed samples on the surface of a triangle mesh
*/


#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "trimesh.h"

#include <vector>
#include <numeric>
#include <algorithm>
using std::vector;
using std::accumulate;
using std::lower_bound;


// Compute the area of each face
void find_face_areas(const TriMesh *m, vector<float> &faceareas)
{
	int nf = m->numfaces;
	faceareas.resize(nf);

	for (int i = 0; i < nf; i++) {
		const face &f = m->faces[i];
		vec facenormal;
		FindNormal(m->vertices[f[0]], m->vertices[f[1]],
			   m->vertices[f[2]], facenormal); // Area-weighted
		faceareas[i] = 0.5f * Len(facenormal);
	}
}


// Compute the cumulative probability distribution over face areas
void compute_cumul_prob(const vector<float> &faceareas,
			vector<float> &cumul_prob)
{
	float totarea = accumulate(faceareas.begin(), faceareas.end(), 0.0f);
	float scale = 1.0f / totarea;

	int nf = faceareas.size();
	cumul_prob.resize(nf);
	float area_so_far = 0;
	for (int i = 0; i < nf; i++) {
		area_so_far += scale * faceareas[i];
		cumul_prob[i] = area_so_far;
	}
}


// Print out a random point on a random triangle
void print_random_point(const TriMesh *m, const vector<float> &cumul_prob,
			FILE *f)
{
	// Pick a random triangle (area-weighted)
	float key = drand48();
	int tri = lower_bound(cumul_prob.begin(), cumul_prob.end(), key) -
		  cumul_prob.begin();

	// Pick a random point
	float w0 = drand48(), w1 = drand48();
	if (w0 + w1 > 1.0f) {
		w0 = 1.0f - w0;
		w1 = 1.0f - w1;
	}
	float w2 = 1.0f - w0 - w1;

	const point &v0 = m->vertices[m->faces[tri][0]];
	const point &v1 = m->vertices[m->faces[tri][1]];
	const point &v2 = m->vertices[m->faces[tri][2]];
	point p = { w0*v0[0] + w1*v1[0] + w2*v2[0],
		    w0*v0[1] + w1*v1[1] + w2*v2[1],
		    w0*v0[2] + w1*v1[2] + w2*v2[2] };

	const vec &n0 = m->normals[m->faces[tri][0]];
	const vec &n1 = m->normals[m->faces[tri][1]];
	const vec &n2 = m->normals[m->faces[tri][2]];
	vec n = { w0*n0[0] + w1*n1[0] + w2*n2[0],
		  w0*n0[1] + w1*n1[1] + w2*n2[1],
		  w0*n0[2] + w1*n1[2] + w2*n2[2] };
	Normalize(n);

	fprintf(f, "%f %f %f %f %f %f\n", p[0], p[1], p[2], n[0], n[1], n[2]);
}


// Usage
void usage(const char *myname)
{
	fprintf(stderr, "Usage: %s model.ply nsamples outfile\n", myname);
	exit(1);
}


// Main: read mesh, compute CPDF of areas, sample
int main(int argc, char *argv[])
{
	if (argc != 4)
		usage(argv[0]);

	const char *meshfile = argv[1];
	int nsamples = atoi(argv[2]);
	const char *outfile = argv[3];

	TriMesh *m = TriMesh::ReadPly(meshfile);
	if (!m)
		usage(argv[0]);
	m->need_normals();

	vector<float> faceareas, cumul_prob;
	find_face_areas(m, faceareas);
	compute_cumul_prob(faceareas, cumul_prob);

	FILE *f = fopen(outfile, "w");
	for (int i = 0; i < nsamples; i++)
		print_random_point(m, cumul_prob, f);
	fclose(f);
}