Code:
/* main module */
#include "vector.h"
#include "mesh.h"
#include "kdtree.h"
#include <stdio.h>
#include <stdlib.h>
int main(void)
{
mesh *m;
kdtree *tree;
FILE *fp = fopen("sphere.zeus", "r");
m = malloc(sizeof(*m));
tree = malloc(sizeof(*tree));
tree->maxdepth = 1;
tree->maxtriangles = 5;
parse_zeus_file(fp, &m);
compute_bounding_box(&m->box, m->vert, m->nvert);
if(init_kdtree(&tree, m))
return 1;
return 0;
}
Code:
/* vector.h */
#ifndef VECTOR_H
#define VECTOR_H
#include <stdio.h>
#include <math.h>
typedef struct
{
double coord[3];
}vector;
double vector_squared_length(vector *);
double vector_length(vector *);
vector *vector_negate(vector *);
double vector_dot(vector *, vector *);
double vector_distance(vector *, vector *);
vector *vector_normalize(vector *);
vector *vector_add(vector *, vector *, vector *);
vector *vector_scale(vector *, double);
vector *vector_sub(vector *, vector *, vector*);
void vector_print(vector *);
vector *vector_cross(vector *, vector *, vector *);
vector *vector_mul_scalar(vector *, double);
#endif
Code:
#include "vector.h"
double vector_squared_length(vector *a)
{
return (a->coord[0] * a->coord[0] + a->coord[1] * a->coord[1] + a->coord[2] * a->coord[2]);
}
double vector_length(vector *a)
{
return (sqrt(vector_squared_length(a)));
}
vector *vector_negate(vector *a)
{
a->coord[0] = -a->coord[0];
a->coord[1] = -a->coord[1];
a->coord[2] = -a->coord[2];
return (a);
}
double vector_dot(vector *a, vector *b)
{
return (a->coord[0] * b->coord[0] + a->coord[1] * b->coord[1] + a->coord[2] * b->coord[2]);
}
double vector_distance(vector *a, vector *b)
{
vector ds;
vector_sub(a, b, &ds);
return (vector_length(&ds));
}
vector *vector_normalize(vector *a)
{
double len;
len = vector_length(a);
if(len != 0)
{
a->coord[0] /= len;
a->coord[1] /= len;
a->coord[2] /= len;
}
return (a);
}
vector *vector_add(vector *a, vector *b, vector *c)
{
c->coord[0] = a->coord[0] + b->coord[0];
c->coord[1] = a->coord[1] + b->coord[1];
c->coord[2] = a->coord[2] + b->coord[2];
return (c);
}
vector *vector_scale(vector *a, double newlength)
{
double length;
length = vector_length(a);
if(length != 0)
{
a->coord[0] *= newlength/length;
a->coord[1] *= newlength/length;
a->coord[2] *= newlength/length;
}
return (a);
}
vector *vector_sub(vector *a, vector *b, vector *c)
{
c->coord[0] = a->coord[0] - b->coord[0];
c->coord[1] = a->coord[1] - b->coord[1];
c->coord[2] = a->coord[2] - b->coord[2];
return (c);
}
void vector_print(vector *a)
{
printf("%f %f %f\n", a->coord[0], a->coord[1] , a->coord[2]);
}
vector *vector_cross(vector *a, vector *b, vector *c)
{
c->coord[0] = (a->coord[1] * b->coord[2]) - (a->coord[2] * b->coord[1]);
c->coord[1] = (a->coord[2] * b->coord[0]) - (a->coord[0] * b->coord[2]);
c->coord[2] = (a->coord[0] * b->coord[1]) - (a->coord[1] * b->coord[0]);
return (c);
}
vector *vector_mul_scalar(vector *a, double scalar)
{
a->coord[0] = a->coord[0] * scalar;
a->coord[1] = a->coord[1] * scalar;
a->coord[2] = a->coord[2] * scalar;
return (a);
}
Code:
/* mesh.h */
#ifndef MESH_H
#define MESH_H
#include "vector.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <float.h>
typedef vector vertex;
typedef struct
{
int v[3];
vector normal;
}triangle;
typedef struct
{
vector maxB, minB;
}bbox;
typedef struct
{
int nvert, ntri;
triangle *tri;
vertex *vert;
bbox box;
}mesh;
int skip_line(FILE *);
int parse_sizes(FILE *, mesh *);
int read_vertex(FILE *, vertex *);
int read_triangle(FILE *, triangle *, int);
int read_triangles(FILE *, mesh *);
int read_vertices(FILE *, mesh *);
int check_end_of_file(FILE *);
int parse_zeus_file(FILE *, mesh **);
void compute_bounding_box(bbox *, vertex *, int);
#endif
Code:
/* mesh.c */
#include "mesh.h"
int skip_line(FILE *fp)
{
int c;
int flag = 0;
while((c = fgetc(fp)) != EOF && c != '\n')
continue;
if(c == EOF)
{
flag = 1;
fprintf(stderr, "%s, %s(), line %d: Error while parsing the file\n", __FILE__, __func__, __LINE__);
return flag;
}
return flag;
}
int parse_sizes(FILE *fp, mesh *m)
{
int nvert, ntri;
int flag = 0;
if(fscanf(fp, "%d%*d%d", &nvert, &ntri) != 2)
{
flag = 1;
fprintf(stderr, "%s, %s(), line %d: Error while reading the number of vertices and/or triangles\n",__FILE__,__func__,__LINE__);
return flag;
}
m->ntri = ntri;
m->nvert = nvert;
if((m->vert = malloc(sizeof(vertex) * m->nvert)) == NULL)
{
flag = 1;
fprintf(stderr, "%s, %s(), line %d: Error while allocating memory for the vertices\n",__FILE__, __func__, __LINE__);
return flag;
}
if((m->tri = malloc(sizeof(triangle) * m->ntri)) == NULL)
{
flag = 1;
fprintf(stderr, "%s, %s(), line %d: Error while allocating memory for the triangles\n",__FILE__, __func__, __LINE__);
return flag;
}
if(skip_line(fp) && skip_line(fp))
{
flag = 1;
fprintf(stderr, "%s, %s(), line %d: Error while reading the file\n",__FILE__, __func__, __LINE__);
return flag;
}
return flag;
}
int read_vertex(FILE *fp, vertex *vp)
{
int flag = 0;
if(fscanf(fp, "%lf %lf %lf", &vp->coord[0], &vp->coord[1], &vp->coord[2]) != 3)
{
flag = 1;
fprintf(stderr, "%s, %s(), line %d: Error while reading the vertex coordinates\n", __FILE__, __func__, __LINE__);
return flag;
}
return flag;
}
int read_triangle(FILE *fp, triangle *tp, int nvert)
{
int flag = 0;
int temp;
if(fscanf(fp, "%d %d %d", &tp->v[0], &tp->v[1], &tp->v[2]) != 3)
{
flag = 1;
fprintf(stderr, "%s, %s(), line %d: Error while reading the triangle vertex indices\n", __FILE__, __func__, __LINE__);
return flag;
}
temp = tp->v[0];
tp->v[0] = tp->v[2];
tp->v[2] = temp;
return flag;
}
int read_vertices(FILE *fp, mesh *m)
{
int nv = 0;
int flag = 0;
while(nv < m->nvert && !read_vertex(fp, &m->vert[nv]))
nv += 1;
if(nv != m->nvert)
{
flag = 1;
fprintf(stderr, "%s, %s(), line %d: Number of vertices read is not consistent with the expected value\n", __FILE__, __func__, __LINE__);
return flag;
}
return flag;
}
int read_triangles(FILE *fp, mesh *m)
{
int nt = 0;
int flag = 0;
vector e1, e2;
while(nt < m->ntri && !read_triangle(fp, &m->tri[nt], m->nvert))
{
vector_sub(&m->vert[m->tri[nt].v[1]], &m->vert[m->tri[nt].v[0]], &e1);
vector_sub(&m->vert[m->tri[nt].v[2]], &m->vert[m->tri[nt].v[0]], &e2);
vector_cross(&e1,&e2, &m->tri[nt].normal);
vector_normalize(&m->tri[nt].normal);
nt += 1;
}
if(nt != m->ntri)
{
flag = 1;
fprintf(stderr, "%s, %s(), line %d: Number of triangles read is not consistent with the expected value\n", __FILE__, __func__, __LINE__);
return flag;
}
return flag;
}
int check_end_of_file(FILE *fp)
{
char eofstring[] = "\t0\tEND\n";
char line[sizeof eofstring];
int flag = 0;
if(skip_line(fp))
{
flag = 1;
fprintf(stderr, "%s, %s(), line %d: Error while reading the file\n", __FILE__, __func__, __LINE__);
return flag;
}
if(fgets(line, sizeof line, fp) == NULL)
{
flag = 1;
fprintf(stderr, "%s, %s(), line %d: Error while trying to read the end of file\n",__FILE__, __func__, __LINE__);
return flag;
}
if(strcmp(line, eofstring) != 0)
{
flag = 1;
fprintf(stderr, "%s, %s(), line %d: Error while trying to read the end of file\n",__FILE__, __func__, __LINE__);
return flag;
}
return flag;
}
int parse_zeus_file(FILE *fp, mesh **m)
{
return (skip_line(fp) || skip_line(fp) || parse_sizes(fp, *m)
|| read_vertices(fp, *m) ||read_triangles(fp, *m) || check_end_of_file(fp));
}
void compute_bounding_box(bbox *b, vertex *vert, int nvert)
{
int i;
vertex p;
b->minB.coord[0] = b->minB.coord[1] = b->minB.coord[2] = DBL_MAX;
b->maxB.coord[0] = b->maxB.coord[1] = b->maxB.coord[2] = -DBL_MAX;
for(i = 0; i < nvert; i++)
{
p = vert[i];
if(b->minB.coord[0] > p.coord[0])
b->minB.coord[0] = p.coord[0];
if(b->maxB.coord[0] < p.coord[0])
b->maxB.coord[0] = p.coord[0];
if(b->minB.coord[1] > p.coord[1])
b->minB.coord[1] = p.coord[1];
if(b->maxB.coord[1] < p.coord[1])
b->maxB.coord[1] = p.coord[1];
if(b->minB.coord[2] > p.coord[2])
b->minB.coord[2] = p.coord[2];
if(b->maxB.coord[2] < p.coord[2])
b->maxB.coord[2] = p.coord[2];
}
}
Code:
/* kdtree.h */
#ifndef KD_TREE_H
#define KD_TREE_H
#include "vector.h"
#include "mesh.h"
#include <stdbool.h>
#define STACK_MAX_SIZE 50
typedef struct trinode_s
{
int tindex;
struct trinode_s *next;
}trinode;
typedef struct kdnode_s
{
bbox box;
int ntriangles;
int nvertices;
int splitplane;
struct kdnode_s *child[2];
trinode *trifirst;
vector **vertexlist;
double split;
}kdnode;
typedef struct
{
kdnode *root;
int maxdepth;
int maxtriangles;
}kdtree;
int init_kdtree(kdtree **, mesh *);
double find_split_point(vertex **, int , int);
int subdivide_kdtree(kdnode *, mesh *, int , int , int , int );
void quicksort(vertex **, int, int , int );
bool vertex_inside_box(vector *, double, int , int);
bool triangle_inside_box(mesh *, double , int , int , int );
int add_to_triangle_list(trinode **, trinode **, int );
#endif
Code:
/* kdtree.c - Error is in this module line no 162 */
#include "kdtree.h"
int add_to_triangle_list(trinode **first, trinode **last, int i)
{
trinode *p;
int flag = 0;
p = malloc(sizeof(*p));
if(p == NULL)
{
fprintf(stderr, "%s, %s(), line %d: Couldn't allocate memory for the kd tree\n", __FILE__, __func__, __LINE__);
flag = 1;
return flag;
}
p->next = NULL;
p->tindex = i;
if(*first == NULL)
{
*first = p;
*last = p;
}
else
{
(*last)->next = p;
*last = p;
}
return flag;
}
bool triangle_inside_box(mesh *m, double split, int index, int axis, int i)
{
bool rc = true;
if(i == 0)
{
if(m->vert[m->tri[index].v[0]].coord[axis] <= split || m->vert[m->tri[index].v[1]].coord[axis] <= split || m->vert[m->tri[index].v[2]].coord[axis] <= split)
;
else
rc = false;
}
else
{
if(m->vert[m->tri[index].v[0]].coord[axis] > split || m->vert[m->tri[index].v[1]].coord[axis] > split || m->vert[m->tri[index].v[2]].coord[axis] > split)
;
else
rc = false;
}
return rc;
}
bool vertex_inside_box(vector *v, double split, int axis, int i)
{
bool rc = true;
if(i == 0)
{
if(v->coord[axis] <= split)
;
else
rc = false;
}
else
{
if(v->coord[axis] > split)
;
else
rc = false;
}
return rc;
}
void quicksort(vertex **a, int left, int right, int axis)
{
int i, j;
vector *pivotpoint;
vector *tempstore;
i = left;
j = right;
pivotpoint = a[(left+right)/2];
while(i<=j)
{
while(a[i]->coord[axis] < pivotpoint->coord[axis])
i++;
while (a[j]->coord[axis] > pivotpoint->coord[axis])
j--;
if(i<=j)
{
tempstore = a[i];
a[i] = a[j];
a[j] = tempstore;
i++;
j--;
}
}
if(left<j)
quicksort(a, left, j, axis);
if(i<right)
quicksort(a, i, right, axis);
return;
}
double find_split_point(vertex **a, int last, int axis)
{
double split;
if(last % 2 == 0)
split = a[last / 2]->coord[axis];
else
split = (a[last/2]->coord[axis] + a[last/2 + 1]->coord[axis])/2.0;
return split;
}
int subdivide_kdtree(kdnode *kd, mesh *m, int depth, int maxdepth, int maxtriangles, int axis)
{
int i, nextaxis;
trinode *ptr, *last;
int j,k;
int flag = 0;
kd->child[0] = kd->child[1] = NULL;
kd->splitplane = -1;
if(kd->ntriangles > maxtriangles && (depth < maxdepth))
{
kd->splitplane = axis;
quicksort(kd->vertexlist, 0, kd->nvertices - 1, axis);
kd->split = find_split_point(kd->vertexlist, kd->nvertices - 1, axis);
for(i = 0; i < 2; i++)
{
kd->child[i] = malloc(sizeof(kdnode));
if(kd->child[i] == NULL)
{
fprintf(stderr, "%s , %s(), line %d: Memory allocation failure\n",__FILE__, __func__, __LINE__);
flag = 1;
return flag;
}
kd->child[i]->box = kd->box;
kd->child[i]->trifirst = NULL;
kd->child[i]->ntriangles = 0;
kd->child[i]->nvertices = 0;
kd->child[i]->box.minB.coord[axis] = kd->box.minB.coord[axis] + i * (kd->split - kd->box.minB.coord[axis]);
kd->child[i]->box.maxB.coord[axis] = kd->split + i * (kd->box.maxB.coord[axis] - kd->split);
ptr = kd->trifirst;
while(ptr != NULL)
{
if(triangle_inside_box(m, kd->split, ptr->tindex, axis, i))
{
if(add_to_triangle_list(&kd->child[i]->trifirst, &last, ptr->tindex))
{
flag = 1;
return flag;
}
else
kd->child[i]->ntriangles += 1;
}
ptr = ptr->next;
}
for(j = 0; j < kd->nvertices; j++)
{
if(vertex_inside_box(kd->vertexlist[j], kd->split, axis, i))
kd->child[i]->nvertices += 1;
}
kd->child[i]->vertexlist = malloc(sizeof(vector *) * (kd->child[i]->nvertices));
k = 0;
for(j = 0; j < kd->nvertices; j++)
{
if(vertex_inside_box(kd->vertexlist[j], kd->split, axis, i))
kd->child[i]->vertexlist[k] = kd->vertexlist[j];
k += 1;
}
nextaxis = (axis + 1) % 3;
flag = subdivide_kdtree(kd->child[i], m, depth + 1, maxdepth, maxtriangles, nextaxis);
}
}
return flag;
}
int init_kdtree(kdtree **tree, mesh *m)
{
int i;
trinode *last;
int flag = 0;
(*tree)->root = malloc(sizeof(kdnode));
if((*tree)->root == NULL)
{
fprintf(stderr, "%s, %s(), line %d: Memory allocation failure\n", __FILE__, __func__, __LINE__);
flag = 1;
return flag;
}
(*tree)->root->box = m->box;
(*tree)->root->ntriangles = m->ntri;
(*tree)->root->trifirst = NULL;
for(i = 0; i < m->ntri; i++)
{
if(add_to_triangle_list(&(*tree)->root->trifirst, &last, i))
{
flag = 1;
return flag;
}
}
(*tree)->root->nvertices = m->nvert;
(*tree)->root->vertexlist = malloc(sizeof(vector *) * m->nvert);
if((*tree)->root->vertexlist == NULL)
{
fprintf(stderr, "%s, %s(), line %d: Memory allocation failure\n", __FILE__, __func__, __LINE__);
flag = 1;
return flag;
}
for(i = 0; i < m->nvert; i++)
(*tree)->root->vertexlist[i] = &m->vert[i];
if(subdivide_kdtree((*tree)->root, m, 0 , (*tree)->maxdepth, (*tree)->maxtriangles, 0))
{
flag = 1;
return flag;
}
return flag;
}