Files
manigraph/include/klein/norm.h
2024-12-03 22:51:50 -06:00

133 lines
2.5 KiB
C

#ifdef KLEIN_NORM_H
#error file included twice
#endif
#define KLEIN_NORM_H
#ifndef KLEIN_H
#warning Please include klein/klein.h
#endif
void klein_normalize(struct klein * klein);
#ifdef KLEIN_IMPLEMENT
static inline
void __calculate_normal( float *p1, float *p2, float *p3, float *normal, unsigned char n)
{ unsigned char i;
float alpha;
float *v1, *v2, *v3;
float *u1, *u2, *u3;
v1 = malloc(n * sizeof(float));
v2 = malloc(n * sizeof(float));
v3 = malloc(n * sizeof(float));
u1 = malloc(n * sizeof(float));
u2 = malloc(n * sizeof(float));
u3 = malloc(n * sizeof(float));
/*
Calculate a normal vector of a plain using Gram-Schmidt process
*/
{
for (i = 0; i < n; ++i)
{
v1[i] = p2[i] - p1[i];
v2[i] = p3[i] - p1[i];
v3[i] = p1[i];
}
for (i = 0; i < n; ++i)
{
u1[i] = v1[i];
}
{
float proj[n];
float dot_v2_u1 = 0.0f, dot_u1_u1 = 0.0f;
for (i = 0; i < n; ++i)
{
dot_v2_u1 += v2[i] * u1[i];
dot_u1_u1 += u1[i] * u1[i];
}
alpha = dot_v2_u1 / dot_u1_u1;
for (i = 0; i < n; ++i)
{
proj[i] = u1[i] * alpha;
u2[i] = v2[i] - proj[i];
}
}
{
float proj1[n], proj2[n];
float dot_v3_u1 = 0.0f, dot_u1_u1 = 0.0f;
float dot_v3_u2 = 0.0f, dot_u2_u2 = 0.0f;
for (i = 0; i < n; ++i)
{
dot_v3_u1 += v3[i] * u1[i];
dot_u1_u1 += u1[i] * u1[i];
}
for (i = 0; i < n; ++i)
{
proj1[i] = u1[i] * (dot_v3_u1 / dot_u1_u1);
}
for (i = 0; i < n; ++i)
{
dot_v3_u2 += v3[i] * u2[i];
dot_u2_u2 += u2[i] * u2[i];
}
for (i = 0; i < n; ++i)
{
proj2[i] = u2[i] * (dot_v3_u2 / dot_u2_u2);
u3[i] = v3[i] - proj1[i] - proj2[i];
}
}
float magnitude = 0.0f;
for (i = 0; i < n; ++i)
{
magnitude += u3[i] * u3[i];
}
magnitude = sqrtf(magnitude);
for (i = 0; i < n; ++i)
{
normal[i] = u3[i] / magnitude;
}
free(v1);
free(v2);
free(v3);
free(u1);
free(u2);
free(u3);
return;
}
}
void klein_normalize(struct klein * klein)
{
unsigned long i;
unsigned char j;
float *norm_vec;
klein->normals = malloc((klein->dim * klein->vertex_size) * sizeof(float));
norm_vec = malloc(klein->dim * sizeof(float));
for (i = 0; i < klein->vertex_size; i += 3 * klein->dim)
{
__calculate_normal(klein->vertex + i, klein->vertex + i + klein->dim, klein->vertex + i + 2 * klein->dim, norm_vec, klein->dim);
for (j = 0; j < klein->dim; ++j )
{
(klein->normals + i)[j]=norm_vec[j];
(klein->normals + i + klein->dim)[j]=norm_vec[j];
(klein->normals + i + 2*klein->dim)[j]=norm_vec[j];
}
}
free(norm_vec);
}
#endif