subdiv.c

subdiv.c

#include <stdio.h>
#include <stdlib.h>
#include <stdarg.h>
#include <string.h>
#include <assert.h>

#include <signal.h>
#include <execinfo.h>

void panic(const char* format, ...) {
	fprintf(stderr, "Panic: ");

	va_list args;
	va_start(args, format);
	vfprintf(stderr, format, args);
	va_end(args);

	#define BACKTRACE_MAX 256
	void* data[BACKTRACE_MAX];
	int n = backtrace(data, BACKTRACE_MAX);

	char** syms = backtrace_symbols(data, n);

	for (int i = 0; i < n; i++) {
		printf("%s\n", syms[i]);
	}

	free(syms);

	exit(EXIT_FAILURE);
}

void* xmalloc(size_t size) {
	void* ptr = malloc(size);
	if (ptr == NULL) {
		panic("Failed to allocate %zu bytes of memory\n", size);
	}
	return ptr;
}

void* xrealloc(void* ptr, size_t size) {
	void* p2 = realloc(ptr, size);
	if (p2 == NULL) {
		panic("Failed to reallocate %ld bytes of memory\n", size);
	}
	return p2;
}

void xfree(void* ptr) {
	free(ptr);
}

typedef struct Vector Vector;

struct Vector {
	void** data;
	size_t capacity;
	size_t count;
};

Vector* Vector_new() {
	Vector* vec = xmalloc(sizeof(Vector));
	vec->data = NULL;
	vec->capacity = 0;
	vec->count = 0;
	return vec;
}

void Vector_delete(Vector* vec) {
	if (vec->data != NULL) {
		xfree(vec->data);
	}
	xfree(vec);
}

void Vector_add(Vector* vec, void* item) {
	if (vec->data == NULL) {
		vec->capacity = 4;
		vec->data = xmalloc(vec->capacity * sizeof(void*));
	}

	if (vec->count == vec->capacity) {
		vec->capacity *= 2;
		vec->data = xrealloc(vec->data, vec->capacity * sizeof(void*));
	}

	vec->data[vec->count++] = item;
}

typedef struct vec3 vec3;

struct vec3 {
	float x;
	float y;
	float z;
};

vec3 vec3_add(vec3 a, vec3 b) {
	return (vec3){a.x + b.x, a.y + b.y, a.z + b.z};
}

vec3 vec3_mul(vec3 a, float s) {
	return (vec3){a.x * s, a.y * s, a.z * s};
}

vec3 vec3_div(vec3 a, float s) {
	return vec3_mul(a, 1 / s);
}

vec3 vec3_mid(vec3 a, vec3 b) {
	return vec3_div(vec3_add(a, b), 2);
}

typedef struct Mesh Mesh;
typedef struct Vertex Vertex;
typedef struct Edge Edge;
typedef struct Face Face;

struct Mesh {
	Vector* vertices;
	Vector* edges;
	Vector* faces;
};

struct Vertex {
	vec3 pos;
	Vector* edges;
	Vector* faces;
	int index;

	Vertex* vertex_point;
};

struct Edge {
	Face* faces[2];
	Vertex* vertices[2];

	Vertex* edge_point;
};

struct Face {
	Vector* vertices;
	Vector* edges;

	Vertex* face_point;
};

Mesh* Mesh_new();
void Mesh_delete(Mesh* mesh);

Vertex* Vertex_new(vec3 pos);
void Vertex_delete(Vertex* vert);

Edge* Edge_new(Face* f0, Face* f1, Vertex* v0, Vertex* v1);
void Edge_delete(Edge* edge);

Face* Face_new();
void Face_delete(Face* face);

Mesh* Mesh_new() {
	Mesh* mesh = xmalloc(sizeof(Mesh));
	mesh->vertices = Vector_new();
	mesh->edges = Vector_new();
	mesh->faces = Vector_new();
	return mesh;
}

void Mesh_delete(Mesh* mesh) {
	for (int i = 0; i < mesh->vertices->count; i++) {
		Vertex_delete(mesh->vertices->data[i]);
	}
	Vector_delete(mesh->vertices);

	for (int i = 0; i < mesh->edges->count; i++) {
		Edge_delete(mesh->edges->data[i]);
	}
	Vector_delete(mesh->edges);

	for (int i = 0; i < mesh->faces->count; i++) {
		Face_delete(mesh->faces->data[i]);
	}
	Vector_delete(mesh->faces);
	xfree(mesh);
}

Vertex* Vertex_new(vec3 pos) {
	// printf("%f %f %f\n", pos.x, pos.y, pos.z);
	
	Vertex* vert = xmalloc(sizeof(Vertex));
	vert->pos = pos;
	vert->edges = Vector_new();
	vert->faces = Vector_new();

	// vert->index = -1;

	return vert;
}

void Vertex_delete(Vertex* vert) {
	Vector_delete(vert->edges);
	Vector_delete(vert->faces);
	xfree(vert);
}

Edge* Edge_new(Face* f0, Face* f1, Vertex* v0, Vertex* v1) {
	// printf("Edge_new %p %p\n", f0, f1);

	Edge* edge = xmalloc(sizeof(Edge));
	edge->faces[0] = f0;
	edge->faces[1] = f1;
	edge->vertices[0] = v0;
	edge->vertices[1] = v1;
	return edge;
}

void Edge_delete(Edge* edge) {
	xfree(edge);
}

vec3 Edge_mid(Edge* edge) {
	return vec3_mid(edge->vertices[0]->pos, edge->vertices[1]->pos);
}

Face* Face_new() {
	Face* face = xmalloc(sizeof(Face));
	face->vertices = Vector_new();
	face->edges = Vector_new();
	return face;
}

void Face_delete(Face* face) {
	Vector_delete(face->vertices);
	Vector_delete(face->edges);
	xfree(face);
}

void Face_compute_face_point(Face* face) {
	vec3 face_point = {0};
	for (int i = 0; i < face->vertices->count; i++) {
		Vertex* vert = face->vertices->data[i];
		face_point = vec3_add(face_point, vert->pos);
	}
	face_point = vec3_div(face_point, face->vertices->count);
	face->face_point = Vertex_new(face_point);
}

void Edge_compute_edge_point(Edge* edge) {
	// printf("edge %p\n", edge->faces[1]);

	vec3 mid = Edge_mid(edge);
	vec3 edge_point;

	if (edge->faces[1] == NULL) {
		edge_point = mid;
	}
	else {
		vec3 fp_mid = vec3_mid(edge->faces[0]->face_point->pos, edge->faces[1]->face_point->pos);
		edge_point = vec3_mid(mid, fp_mid);
	}

	edge->edge_point = Vertex_new(edge_point);
}

void Vertex_compute_vertex_point(Vertex* vert) {
	if (vert->faces->count == vert->edges->count) {
		vec3 face_avg = {0};
		for (int i = 0; i < vert->faces->count; i++) {
			Face* face = vert->faces->data[i];
			face_avg = vec3_add(face_avg, face->face_point->pos);
		}
		face_avg = vec3_div(face_avg, vert->faces->count);

		vec3 edge_avg = {0};
		for (int i = 0; i < vert->edges->count; i++) {
			Edge* edge = vert->edges->data[i];
			edge_avg = vec3_add(edge_avg, Edge_mid(edge));
		}
		edge_avg = vec3_div(edge_avg, vert->edges->count);

		vec3 vertex_point = vec3_div(vec3_add(face_avg, vec3_add(vec3_mul(edge_avg, 2), vec3_mul(vert->pos, (float)vert->faces->count - 3))), vert->faces->count);
		vert->vertex_point = Vertex_new(vertex_point);
	}
	else {
		vec3 vertex_point = vert->pos;
		int n = 1;
		for (int i = 0; i < vert->edges->count; i++) {
			Edge* edge = vert->edges->data[i];
			if (edge->faces[1] == NULL) {
				vertex_point = vec3_add(vertex_point, Edge_mid(edge));
				n++;
			}
		}
		vertex_point = vec3_div(vertex_point, n);
		vert->vertex_point = Vertex_new(vertex_point);
	}
}

Mesh* Mesh_catmull_clark_subdivide(Mesh* mesh) {
	// printf("face point\n");
	for (int i = 0; i < mesh->faces->count; i++) {
		Face_compute_face_point(mesh->faces->data[i]);
	}

	// printf("edge point\n");
	for (int i = 0; i < mesh->edges->count; i++) {
		Edge_compute_edge_point(mesh->edges->data[i]);
	}

	// printf("vertex point\n");
	int numVerts = mesh->vertices->count;
	for (int i = 0; i < numVerts; i++) {
		Vertex_compute_vertex_point(mesh->vertices->data[i]);
	}

	// printf("done\n");

	Mesh* subdiv = Mesh_new();

	// Add vertices

	int vertexi = 0;

	for (int i = 0; i < mesh->vertices->count; i++) {
		Vertex* vert = mesh->vertices->data[i];
		vert->vertex_point->index = vertexi++;
		Vector_add(subdiv->vertices, vert->vertex_point);
	}

	for (int i = 0; i < mesh->edges->count; i++) {
		Edge* edge = mesh->edges->data[i];
		edge->edge_point->index = vertexi++;
		Vector_add(subdiv->vertices, edge->edge_point);
	}

	for (int i = 0; i < mesh->faces->count; i++) {
		Face* face = mesh->faces->data[i];
		face->face_point->index = vertexi++;
		Vector_add(subdiv->vertices, face->face_point);
	}

	// Add faces

	for (int i = 0; i < mesh->faces->count; i++) {
		Face* face = mesh->faces->data[i];
		for (int j = 0; j < face->edges->count; j++) {
			Edge* edge0 = face->edges->data[j];
			Edge* edge1;
			if (j == face->edges->count - 1) {
				edge1 = face->edges->data[0];
			}
			else {
				edge1 = face->edges->data[j + 1];
			}

			// printf("%p, %p, %p, %p\n", edge0->vertices[0], edge0->vertices[1], edge1->vertices[0], edge1->vertices[1]);

			Vertex* vert = NULL;
			if (edge0->vertices[0] == edge1->vertices[0]) {
				vert = edge0->vertices[0];
			}
			else if (edge0->vertices[0] == edge1->vertices[1]) {
				vert = edge0->vertices[0];
			}
			else if (edge0->vertices[1] == edge1->vertices[0]) {
				vert = edge0->vertices[1];
			}
			else if (edge0->vertices[1] == edge1->vertices[1]) {
				vert = edge0->vertices[1];
			}
			assert(vert != NULL);
			vert = vert->vertex_point;

			Face* newFace = Face_new();

			Vector_add(newFace->vertices, face->face_point);
			Vector_add(newFace->vertices, edge0->edge_point);
			Vector_add(newFace->vertices, vert);
			Vector_add(newFace->vertices, edge1->edge_point);

			Vector_add(subdiv->faces, newFace);
		}
	}

	return subdiv;
}

Mesh* Mesh_load(const char* path) {
	FILE* file = fopen(path, "r");
	if (file == NULL) {
		panic("Unable to open file %s\n", path);
	}

	Mesh* mesh = Mesh_new();

	char line[256];
	while (fgets(line, sizeof(line), file) != NULL) {
		if (line[0] == 'v') {
			vec3 pos;
			sscanf(line, "v %f %f %f", &pos.x, &pos.y, &pos.z);
			Vector_add(mesh->vertices, Vertex_new(pos));
		}
		else if (line[0] == 'f') {
			Vector* indices = Vector_new();

			char* ptr = strchr(line, ' ');
			while (ptr != NULL) {
				int x = strtol(ptr, &ptr, 10) - 1;
				
				int* xptr = xmalloc(sizeof(int));
				*xptr = x;
				Vector_add(indices, xptr);

				ptr = strchr(ptr, ' ');
			}

			Face* face = Face_new();

			for (int i = 0; i < indices->count; i++) {
				int index = *((int*)indices->data[i]);
				Vertex* vert = mesh->vertices->data[index];
				
				Vector_add(face->vertices, vert);
				Vector_add(vert->faces, face);

				xfree(indices->data[i]);
			}
			Vector_delete(indices);

			Vector_add(mesh->faces, face);
		}
	}

	for (int i = 0; i < mesh->faces->count; i++) {
		Face* face = mesh->faces->data[i];

		// printf("face has %zu vertices\n", face->vertices->count);

		for (int j = 0; j < face->vertices->count; j++) {
			Vertex* vert0 = face->vertices->data[j];
			Vertex* vert1;
			if (j == face->vertices->count - 1) {
				vert1 = face->vertices->data[0];
			}
			else {
				vert1 = face->vertices->data[j + 1];
			}

			Edge* edge = NULL;
			for (int k = 0; k < mesh->edges->count; k++) {
				Edge* e = mesh->edges->data[k];
				if ((e->vertices[0] == vert0 && e->vertices[1] == vert1) || (e->vertices[0] == vert1 && e->vertices[1] == vert0)) {
					edge = e;
					break;
				}
			}
			if (edge == NULL) {
				// printf("find shared faces %zu %zu ", vert0->faces->count, vert1->faces->count);

				Face* face0 = NULL;
				Face* face1 = NULL;

				for (int k = 0; k < vert0->faces->count; k++) {
					Face* f0 = vert0->faces->data[k];

					for (int l = 0; l < vert1->faces->count; l++) {
						Face* f1 = vert1->faces->data[l];

						if (f0 == f1) {
							if (face0 == NULL) {
								face0 = f0;
								// printf("first ");
							}
							else {
								face1 = f0;
								// printf("second ");
							}
						}
					}
				}

				// printf("done ");

				edge = Edge_new(face0, face1, vert0, vert1);
				Vector_add(mesh->edges, edge);
			}

			Vector_add(face->edges, edge);
			Vector_add(vert0->edges, edge);
		}
	}

	fclose(file);
	
	return mesh;
}

void Mesh_save(Mesh* mesh, const char* path) {
	FILE* file = fopen(path, "w");
	if (file == NULL) {
		panic("Unable to open file %s\n", path);
	}

	// printf("save verts %d\n", mesh->vertices->count);
	for (int i = 0; i < mesh->vertices->count; i++) {
		Vertex* vert = mesh->vertices->data[i];
		vec3 pos = vert->pos;
		// printf("%p %f %f %f\n", vert, pos.x, pos.y, pos.z);
		fprintf(file, "v %f %f %f\n", pos.x, pos.y, pos.z);
	}

	// printf("save faces\n");
	for (int i = 0; i < mesh->faces->count; i++) {
		Face* face = mesh->faces->data[i];
		fprintf(file, "f");
		// printf("count is %zu\n", face->vertices->count);
		for (int j = 0; j < face->vertices->count; j++) {
			Vertex* vert = face->vertices->data[j];
			fprintf(file, " %d", vert->index + 1);
		}
		fprintf(file, "\n");
	}

	fclose(file);
}

void sigsegv_func(int signo) {
	panic("Segmentation fault\n");
}

int main(int argc, char** argv) {
	signal(SIGSEGV, sigsegv_func);

	if (argc != 3) {
		return 1;
	}

	// printf("loading\n");
	Mesh* mesh = Mesh_load(argv[1]);
	// printf("%zu vertices %zu edges %zu faces\n", mesh->vertices->count, mesh->edges->count, mesh->faces->count);
	
	// printf("subdiv\n");
	Mesh* mesh2 = Mesh_catmull_clark_subdivide(mesh);
	// printf("%zu vertices %zu edges %zu faces\n", mesh2->vertices->count, mesh2->edges->count, mesh2->faces->count);
	// printf("save\n");
	Mesh_save(mesh2, argv[2]);

	Mesh_delete(mesh);
	Mesh_delete(mesh2);

	return 0;
}