<commit_before>//
//    Paraboloid Fitting of a point cloud
// 
//  first we center and line up the point cloud along its principle axes (higher two eigen values determine x,y plane)
//  then we fit parabaloid in local space using z as the up/height axis
//
//  Note:  when you move a point, it end up uniformly shifting all the points since we constantly re-center on 0,0,0, align as mentioned, and shift z to be tangent to xy plane.  this can seem a bit confusing. 
// 

#include <exception>
#include <iostream>
#include <random>

#include <geometric.h>
#include <glwin.h>
#include <mesh.h>


// this routine taken from testcov:  
std::pair<Pose, float3> PrincipalAxes(const std::vector<float3> &points)  // returns principal axes as a pose and population's variance along pose's local x,y,z
{
	float3   com(0, 0, 0);
	float3x3 cov;
	for (auto p : points)
		com += p;
	com /= (float)points.size();
	for (auto p : points)
		cov += outerprod(p - com, p - com);
	cov /= (float)points.size();
	auto q = Diagonalizer(cov);
	return std::make_pair<Pose, float3>({ com, q }, Diagonal(mul(transpose(qgetmatrix(q)), cov, qgetmatrix(q))));
}

float4 ParabloidFit(const std::vector<float3> &pts)  // least squares fitting quadratic patch centered at x,y==0,0
{
	float4x4 m;
	float4 b;
	for (auto &p : pts)
	{
		float4 v(p.x*p.x, p.y*p.y, p.x*p.y, 1.0f);
		m += outerprod(v, v);
		b += v*p.z;
	}
	return mul(inverse(m), b);   // returns 'h' best fits   z = h.x * x*x + h.y * y*y + h.z * x*y + h.w  hessian and zoffset 
}


std::default_random_engine random_number_generator;
float3 vrand() // output range is from -1 to 1
{
	auto r = [](){return std::uniform_real<float>(-1.0f, 1.0f)(random_number_generator); };
	return{ r(), r(), r() };
}
float3 vrand(const float3 &scale) { return cmul(scale, vrand()); }
float3 vsqrt(const float3 &v)     { return{ sqrtf(v.x), sqrtf(v.y), sqrtf(v.z) }; }

std::vector<float3> RandomParabolicCloud(float3 range = float3(1.5f, 1.0f, 0.25f), int n = 15)
{
	float2 k = vrand().xy()*4.0f;   // we randomize the curvatures used to generate the point cloud.
	Pose randpose({ 1, 1, 1 }, normalize(float4(vrand(), 1.0f)));  // randomly pose the initial data to show that we are indeed normalizing a potentially generic point cloud 
	std::vector<float3> points(n);
	std::transform(points.begin(), points.end(), points.begin(), [&randpose, &range,&k](const float3)->float3 {auto v = vrand().xy() ; return randpose*cmul(range, float3(v, v.x*v.x*k.x + v.y*v.y*k.y)); });
	return points;
}


// some typical opengl drawing support routines
//

void glWirePatch(std::function<float3(float2)> m)
{
	auto f = [](std::function<float3(float2)> m)  
	{
		for (float x = 0; x <= 1.0f; x += 1.0f / 16)
		{
			glBegin(GL_LINE_STRIP);
			for (float y = 0; y <= 1.0f; y += 1.0f / 16)
				glVertex3fv(m({ x,y }));
			glEnd();
		}
	};
	f(m); 
	f([m](float2 c) {return m({ c.y,c.x });});
}

void glHeightField(float2 bmin, float2 bmax, std::function<float(float2)> h, float3 c = { 1,1,1 })
{
	glColor3fv(c);
	glWirePatch( [&h,&bmin,&bmax](float2 c)->float3 {auto p = bmin + cmul(c, bmax - bmin); return{ p.x,p.y,h(p) };} );
}


void glAxis()
{
	glPushAttrib(GL_ALL_ATTRIB_BITS);
	glLineWidth(1.5f);
	glBegin(GL_LINES);
	for (int i : {0, 1, 2})
	{
		float3 v(0, 0, 0);
		v[i] = 0.6f;
		glColor3fv(v);
		glVertex3fv({ 0, 0, 0 });
		glVertex3fv(v);
	}
	glEnd();
	glPopAttrib();
}

void glcolorbox(const float3 &r, const Pose &p)   // p = { { 0, 0, 0 }, { 0, 0, 0, 1 } })
{
	glPushMatrix();
	glMultMatrixf(p.Matrix());
	glBegin(GL_QUADS);
	for (int m : {0, 1}) for (int i : {0, 1, 2})
	{
		int i1 = (i + 1 + m) % 3;
		int i2 = (i + 2 - m) % 3;
		float3 u, v, w;
		u[i1] = r[i1];
		v[i2] = r[i2];
		w[i] = (m) ? -1.0f : 1.0f;
		float3 a((float)m, (float)m, (float)m);
		a[i] = 1 - a[i];
		glColor3fv(a);
		glNormal3fv(w);
		float2 corners[] = { { -1.0f, -1.0f }, { 1.0f, -1.0f }, { 1.0f, 1.0f }, { -1.0f, 1.0f } };  // ccw order
		for (float2 t : corners)
			glTexCoord2fv(t), glVertex3fv(w*r[i] + u*t.x + v*t.y);
	}
	glEnd();
	glPopMatrix();
}
void glcolorbox(const float &r, const Pose &p) { glcolorbox({ r, r, r }, p); }



float DegToRad(float a) { return a*3.14159f / 180.0f; }

int APIENTRY WinMain(HINSTANCE hCurrentInst, HINSTANCE hPreviousInst, LPSTR lpszCmdLine, int nCmdShow) try
{ 
	Pose camera = Pose({0,0,0},{sin(DegToRad(60.0f/2.0f)),0,0,cos(DegToRad(60.0f/2.0f))}) * Pose({0,0,6},{0,0,0,1}) ;
	float3 focuspoint(0, 0, 0);
	float3 mousevec_prev;
	float4 model_orientation(0, 0, 0, 1);
	float dt = 0.01f, t = 0;
	float3 *selected = NULL;
	bool moving_enabled = 0;  // to move vertices with Left mouse drag
	float boxr = 0.025f;
	std::vector<float4> planes = { { 1, 0, 0, 0 }, { 0, 1, 0, 0 }, { 0, 0, 1, 0 }, { -1, 0, 0, 0 }, { 0, -1, 0, 0 }, { 0, 0, -1, 0 } };
	for (auto &p : planes)
		p.w = -boxr;
	std::vector<float3> points = RandomParabolicCloud();

	GLWin glwin("ParaboloidFit  least Squares quadratic surface fit",800,600);
	glwin.keyboardfunc = [&](int key, int, int)
	{
		moving_enabled = key == 'm' != moving_enabled;  // user vertex dragging
		if (key == ' ')
			points = RandomParabolicCloud();
	};
	while (glwin.WindowUp())
	{
		t = t + dt;  // advance our global time    t is in 0..1
		if (t > 1.0f)
			t = 0.0f;

		// user interaction: 
		float3 ray = qrot(camera.orientation, normalize(glwin.MouseVector));   // for mouse selection
		float3 v1 = camera.position + ray*100.0f;
		if (!glwin.MouseState)  // note that we figure out what is being selected only when the mouse is up
		{
			selected = NULL;
			for (float3 &p : points) if (moving_enabled)
			{
				if (auto h = ConvexHitCheck(planes, Pose(p, { 0, 0, 0, 1 }), camera.position, v1))<commit_msg>more uniform randomish paraboloidic point cloud<commit_after>//
//    Paraboloid Fitting of a point cloud
// 
//  first we center and line up the point cloud along its principle axes (higher two eigen values determine x,y plane)
//  then we fit parabaloid in local space using z as the up/height axis
//
//  Note:  when you move a point, it end up uniformly shifting all the points since we constantly re-center on 0,0,0, align as mentioned, and shift z to be tangent to xy plane.  this can seem a bit confusing. 
// 

#include <exception>
#include <iostream>
#include <random>

#include <geometric.h>
#include <glwin.h>
#include <mesh.h>


// this routine taken from testcov:  
std::pair<Pose, float3> PrincipalAxes(const std::vector<float3> &points)  // returns principal axes as a pose and population's variance along pose's local x,y,z
{
	float3   com(0, 0, 0);
	float3x3 cov;
	for (auto p : points)
		com += p;
	com /= (float)points.size();
	for (auto p : points)
		cov += outerprod(p - com, p - com);
	cov /= (float)points.size();
	auto q = Diagonalizer(cov);
	return std::make_pair<Pose, float3>({ com, q }, Diagonal(mul(transpose(qgetmatrix(q)), cov, qgetmatrix(q))));
}

float4 ParabloidFit(const std::vector<float3> &pts)  // least squares fitting quadratic patch centered at x,y==0,0
{
	float4x4 m;
	float4 b;
	for (auto &p : pts)
	{
		float4 v(p.x*p.x, p.y*p.y, p.x*p.y, 1.0f);
		m += outerprod(v, v);
		b += v*p.z;
	}
	return mul(inverse(m), b);   // returns 'h' best fits   z = h.x * x*x + h.y * y*y + h.z * x*y + h.w  hessian and zoffset 
}


std::default_random_engine random_number_generator;
float3 vrand() // output range is from -1 to 1
{
	auto r = [](){return std::uniform_real<float>(-1.0f, 1.0f)(random_number_generator); };
	return{ r(), r(), r() };
}
float3 vrand(const float3 &scale) { return cmul(scale, vrand()); }
float3 vsqrt(const float3 &v)     { return{ sqrtf(v.x), sqrtf(v.y), sqrtf(v.z) }; }

std::vector<float3> RandomParabolicCloud(float3 range = float3(1.5f, 1.0f, 0.25f))
{
	float3 k = vrand()*4.0f;   // we randomize the curvatures used to generate the point cloud.
	Pose randpose({ 1, 1, 1 }, normalize(float4(vrand(), 1.0f)));  // randomly pose the initial data to show that we are indeed normalizing a potentially generic point cloud 
	std::vector<float3> points;
	for (float y = -0.5f; y <= 0.5f; y += 0.25f) for (float x = -0.5f; x <= 0.5f; x += 0.25f)
		points.push_back(vrand()*0.1f+ float3(x,y,0.0f));
	std::transform(points.begin(), points.end(), points.begin(), [&randpose, &range,&k](const float3 v)->float3 { return randpose*cmul(range, float3(v.xy(), v.x*v.x*k.x + v.y*v.y*k.y +k.z*v.x*v.y)); });
	return points;
}


// some typical opengl drawing support routines
//

void glWirePatch(std::function<float3(float2)> m)
{
	auto f = [](std::function<float3(float2)> m)  
	{
		for (float x = 0; x <= 1.0f; x += 1.0f / 16)
		{
			glBegin(GL_LINE_STRIP);
			for (float y = 0; y <= 1.0f; y += 1.0f / 16)
				glVertex3fv(m({ x,y }));
			glEnd();
		}
	};
	f(m); 
	f([m](float2 c) {return m({ c.y,c.x });});
}

void glHeightField(float2 bmin, float2 bmax, std::function<float(float2)> h, float3 c = { 1,1,1 })
{
	glColor3fv(c);
	glWirePatch( [&h,&bmin,&bmax](float2 c)->float3 {auto p = bmin + cmul(c, bmax - bmin); return{ p.x,p.y,h(p) };} );
}


void glAxis()
{
	glPushAttrib(GL_ALL_ATTRIB_BITS);
	glLineWidth(1.5f);
	glBegin(GL_LINES);
	for (int i : {0, 1, 2})
	{
		float3 v(0, 0, 0);
		v[i] = 0.6f;
		glColor3fv(v);
		glVertex3fv({ 0, 0, 0 });
		glVertex3fv(v);
	}
	glEnd();
	glPopAttrib();
}

void glcolorbox(const float3 &r, const Pose &p)   // p = { { 0, 0, 0 }, { 0, 0, 0, 1 } })
{
	glPushMatrix();
	glMultMatrixf(p.Matrix());
	glBegin(GL_QUADS);
	for (int m : {0, 1}) for (int i : {0, 1, 2})
	{
		int i1 = (i + 1 + m) % 3;
		int i2 = (i + 2 - m) % 3;
		float3 u, v, w;
		u[i1] = r[i1];
		v[i2] = r[i2];
		w[i] = (m) ? -1.0f : 1.0f;
		float3 a((float)m, (float)m, (float)m);
		a[i] = 1 - a[i];
		glColor3fv(a);
		glNormal3fv(w);
		float2 corners[] = { { -1.0f, -1.0f }, { 1.0f, -1.0f }, { 1.0f, 1.0f }, { -1.0f, 1.0f } };  // ccw order
		for (float2 t : corners)
			glTexCoord2fv(t), glVertex3fv(w*r[i] + u*t.x + v*t.y);
	}
	glEnd();
	glPopMatrix();
}
void glcolorbox(const float &r, const Pose &p) { glcolorbox({ r, r, r }, p); }



float DegToRad(float a) { return a*3.14159f / 180.0f; }

int APIENTRY WinMain(HINSTANCE hCurrentInst, HINSTANCE hPreviousInst, LPSTR lpszCmdLine, int nCmdShow) try
{ 
	Pose camera = Pose({0,0,0},{sin(DegToRad(60.0f/2.0f)),0,0,cos(DegToRad(60.0f/2.0f))}) * Pose({0,0,6},{0,0,0,1}) ;
	float3 focuspoint(0, 0, 0);
	float3 mousevec_prev;
	float4 model_orientation(0, 0, 0, 1);
	float3 *selected = NULL;
	bool moving_enabled = 0;  // to move vertices with Left mouse drag
	float boxr = 0.025f;
	std::vector<float4> planes = { { 1, 0, 0, 0 }, { 0, 1, 0, 0 }, { 0, 0, 1, 0 }, { -1, 0, 0, 0 }, { 0, -1, 0, 0 }, { 0, 0, -1, 0 } };
	for (auto &p : planes)
		p.w = -boxr;
	std::vector<float3> points = RandomParabolicCloud();

	GLWin glwin("ParaboloidFit  least Squares quadratic surface fit",800,600);
	glwin.keyboardfunc = [&](int key, int, int)
	{
		moving_enabled = key == 'm' != moving_enabled;  // user vertex dragging
		if (key == ' ')
			points = RandomParabolicCloud();
	};
	while (glwin.WindowUp())
	{

		// user interaction: 
		float3 ray = qrot(camera.orientation, normalize(glwin.MouseVector));   // for mouse selection
		float3 v1 = camera.position + ray*100.0f;
		if (!glwin.MouseState)  // note that we figure out what is being selected only when the mouse is up
		{
			selected = NULL;
			for (float3 &p : points) if (moving_enabled)
			{
				if (auto h = ConvexHitCheck(planes, Pose(p, { 0, 0, 0, 1 }), camera.position, v1))
{'fraction_non_alphanumeric': 0.18379706837531817, 'fraction_numerical': 0.046256473273062405, 'mean_word_length': 3.5376344086021505, 'pattern_counts': {'":': 0, '<': 45, '<?xml version=': 0, '>': 44, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1164883', 'n_tokens_mistral': 5174, 'n_tokens_neox': 4738, 'n_words': 1313}