SaraP
cbec90699f
- creato il progetto in Visual Studio per compilare come libreria statica - modifiche al codice originale per integrarlo nelle nostre librerie.
362 lines
13 KiB
C++
362 lines
13 KiB
C++
/*****************************************************************************/
|
|
/* */
|
|
/* F I S T : Fast, Industrial-Strength Triangulation */
|
|
/* */
|
|
/*****************************************************************************/
|
|
/* */
|
|
/* (C) Martin Held */
|
|
/* (C) Universitaet Salzburg, Salzburg, Austria */
|
|
/* */
|
|
/* This code is not in the public domain. All rights reserved! Please make */
|
|
/* sure to read the full copyright statement contained in api_functions.cpp. */
|
|
/* */
|
|
/*****************************************************************************/
|
|
|
|
/* */
|
|
/* get standard libraries */
|
|
/* */
|
|
#include <stdlib.h>
|
|
#include <stdio.h>
|
|
#include <math.h>
|
|
#include <assert.h>
|
|
|
|
/* */
|
|
/* get my header files */
|
|
/* */
|
|
#include "fpkernel.h"
|
|
#include "martin.h"
|
|
#include "defs.h"
|
|
#include "list.h"
|
|
#include "header.h"
|
|
#include "basic.h"
|
|
#include "vertex.h"
|
|
|
|
/* */
|
|
/* prototypes of functions provided in this file */
|
|
/* */
|
|
void ProjectFace(global_struct *all, int loop_min, int loop_max);
|
|
|
|
|
|
#define BLOCK_SIZE 32768
|
|
|
|
typedef struct {
|
|
double element[4][4];
|
|
} Matrix4;
|
|
|
|
|
|
/* note: pnt and new_pnt must not be the same memory cell! */
|
|
#define MatrixTransformation(matrix, pnt, new_pnt) \
|
|
{ \
|
|
(new_pnt).x = (pnt).x * (matrix).element[0][0] \
|
|
+ (pnt).y * (matrix).element[1][0] \
|
|
+ (pnt).z * (matrix).element[2][0] \
|
|
+ (matrix).element[3][0]; \
|
|
(new_pnt).y = (pnt).x * (matrix).element[0][1] \
|
|
+ (pnt).y * (matrix).element[1][1] \
|
|
+ (pnt).z * (matrix).element[2][1] \
|
|
+ (matrix).element[3][1]; \
|
|
(new_pnt).z = (pnt).x * (matrix).element[0][2] \
|
|
+ (pnt).y * (matrix).element[1][2] \
|
|
+ (pnt).z * (matrix).element[2][2] \
|
|
+ (matrix).element[3][2]; \
|
|
}
|
|
|
|
/* */
|
|
/* this function projects the vertices of the polygons referenced by */
|
|
/* loops[loop_min,..,loop_max-1] to an approximating plane. */
|
|
/* */
|
|
void ProjectFace(global_struct *all, int loop_min, int loop_max)
|
|
{
|
|
boolean DetermineNormal(global_struct *all, list_ind ind, vertex *normal);
|
|
void ProjectPoints(global_struct *all, int i1, int i2, vertex normal);
|
|
vertex normal, nr;
|
|
int i;
|
|
double d;
|
|
|
|
normal.x = normal.y = normal.z = C_0_0;
|
|
|
|
/* */
|
|
/* determine the normal of the plane onto which the points get projected */
|
|
/* */
|
|
for (i = loop_min; i < loop_max; ++i) {
|
|
if (DetermineNormal(all, all->c_list.loops[i], &nr)) {
|
|
if (DotProduct(normal, nr) < C_0_0) {
|
|
InvertVector(nr);
|
|
}
|
|
VectorAdd(normal, nr, normal);
|
|
}
|
|
}
|
|
|
|
if (loop_min <= (loop_max - 1)) {
|
|
d = Length_l2(normal);
|
|
if (gt(d, ZERO)) {
|
|
DivScalar(d, normal);
|
|
}
|
|
else {
|
|
if (all->rt_opt.verbose) {
|
|
FIST_Warning("ProjectFace() - zero-length normal vector!\n");
|
|
}
|
|
normal.x = normal.y = C_0_0;
|
|
normal.z = C_1_0;
|
|
}
|
|
}
|
|
|
|
/* COREBACKEND: Simplify the expression tree of the normal vector, */
|
|
/* otherwise the final point coordinates' expression tree is sometimes */
|
|
/* too big for at least CORE.version <= 2.0.8 */
|
|
#ifdef WITH_COREBACKEND
|
|
normal.x = TO_MDOUBLE(normal.x);
|
|
normal.y = TO_MDOUBLE(normal.y);
|
|
normal.z = TO_MDOUBLE(normal.z);
|
|
#endif
|
|
|
|
/* */
|
|
/* project the points onto this plane. the projected points are stored in */
|
|
/* the array `points[0,..,num_pnts]' */
|
|
/* */
|
|
ProjectPoints(all, loop_min, loop_max, normal);
|
|
|
|
return;
|
|
}
|
|
|
|
|
|
/* */
|
|
/* this function computes the average of all normals defined by triples of */
|
|
/* successive vertices of the polygon. we'll see whether this is a good */
|
|
/* heuristic for finding a suitable plane normal... */
|
|
/* */
|
|
boolean DetermineNormal(global_struct *all, list_ind ind, vertex *normal)
|
|
{
|
|
listdef *list = &all->c_list;
|
|
vertexdef *vert = &all->c_vertex;
|
|
|
|
#ifdef SUM_OF_NORMALS
|
|
vertex nr, pq, pr;
|
|
list_ind ind0, ind1, ind2;
|
|
int i0, i1, i2;
|
|
double d;
|
|
|
|
ind1 = ind;
|
|
i1 = GetIndex(list, ind1);
|
|
ind0 = GetPrevNode(list, ind1);
|
|
i0 = GetIndex(list, ind0);
|
|
ind2 = GetNextNode(list, ind1);
|
|
i2 = GetIndex(list, ind2);
|
|
VectorSub(vert->vertices[i0], vert->vertices[i1], pq);
|
|
VectorSub(vert->vertices[i2], vert->vertices[i1], pr);
|
|
VectorProduct(pq, pr, nr);
|
|
d = Length_l2(nr);
|
|
if (gt(d, ZERO)) {
|
|
DivScalar(d, nr);
|
|
*normal = nr;
|
|
}
|
|
else {
|
|
normal->x = normal->y = normal->z = C_0_0;
|
|
}
|
|
pq = pr;
|
|
ind1 = ind2;
|
|
ind2 = GetNextNode(list, ind1);
|
|
i2 = GetIndex(list, ind2);
|
|
while (ind1 != ind) {
|
|
VectorSub(vert->vertices[i2], vert->vertices[i1], pr);
|
|
VectorProduct(pq, pr, nr);
|
|
d = Length_l2(nr);
|
|
if (gt(d, ZERO)) {
|
|
DivScalar(d, nr);
|
|
if (DotProduct(*normal, nr) < C_0_0) {
|
|
InvertVector(nr);
|
|
}
|
|
VectorAdd(*normal, nr , *normal);
|
|
}
|
|
pq = pr;
|
|
ind1 = ind2;
|
|
ind2 = GetNextNode(list, ind1);
|
|
ind2 = GetIndex(list, ind2);
|
|
}
|
|
|
|
d = Length_l2(*normal);
|
|
if (gt(d, ZERO)) {
|
|
DivScalar(d, *normal);
|
|
return true;
|
|
}
|
|
else {
|
|
if (all->rt_opt.verbose) {
|
|
FIST_Warning("DetermineNormal() - zero-length normal vector!\n");
|
|
}
|
|
normal->x = normal->y = C_0_0;
|
|
normal->z = C_1_0;
|
|
return false;
|
|
}
|
|
#else
|
|
vertex nr, pq, pr;
|
|
vertex x_normal, y_normal, z_normal;
|
|
list_ind ind0, ind1, ind2;
|
|
int i0, i1, i2;
|
|
double d;
|
|
|
|
x_normal.x = x_normal.y = x_normal.z = C_0_0;
|
|
y_normal.x = y_normal.y = y_normal.z = C_0_0;
|
|
z_normal.x = z_normal.y = z_normal.z = C_0_0;
|
|
|
|
ind1 = ind;
|
|
i1 = GetIndex(list, ind1);
|
|
ind0 = GetPrevNode(list, ind1);
|
|
i0 = GetIndex(list, ind0);
|
|
VectorSub(vert->vertices[i0], vert->vertices[i1], pq);
|
|
ind2 = GetNextNode(list, ind1);
|
|
i2 = GetIndex(list, ind2);
|
|
do {
|
|
VectorSub(vert->vertices[i2], vert->vertices[i1], pr);
|
|
VectorProduct(pq, pr, nr);
|
|
d = Length_l2(nr);
|
|
if (gt(d, ZERO)) {
|
|
DivScalar(d, nr);
|
|
if (Abs(nr.x) <= Abs(nr.y)) {
|
|
if (Abs(nr.y) <= Abs(nr.z)) {
|
|
if (nr.z >= C_0_0) {VectorAdd(z_normal, nr, z_normal);}
|
|
else {VectorSub(z_normal, nr, z_normal);}
|
|
}
|
|
else {
|
|
if (nr.y >= C_0_0) {VectorAdd(y_normal, nr, y_normal);}
|
|
else {VectorSub(y_normal, nr, y_normal);}
|
|
}
|
|
}
|
|
else {
|
|
if (Abs(nr.x) <= Abs(nr.z)) {
|
|
if (nr.z >= C_0_0) {VectorAdd(z_normal, nr, z_normal);}
|
|
else {VectorSub(z_normal, nr, z_normal);}
|
|
}
|
|
else {
|
|
if (nr.x >= C_0_0) {VectorAdd(x_normal, nr, x_normal);}
|
|
else {VectorSub(x_normal, nr, x_normal);}
|
|
}
|
|
}
|
|
}
|
|
pq = pr;
|
|
ind1 = ind2;
|
|
ind2 = GetNextNode(list, ind1);
|
|
i2 = GetIndex(list, ind2);
|
|
} while (ind1 != ind);
|
|
|
|
if (DotProduct(z_normal, x_normal) < C_0_0) {
|
|
InvertVector(x_normal);
|
|
}
|
|
VectorAdd(z_normal, x_normal, z_normal);
|
|
if (DotProduct(z_normal, y_normal) < C_0_0) {
|
|
InvertVector(y_normal);
|
|
}
|
|
VectorAdd(z_normal, y_normal, *normal);
|
|
d = Length_l2(*normal);
|
|
if (gt(d, ZERO)) {
|
|
DivScalar(d, *normal);
|
|
return true;
|
|
}
|
|
else {
|
|
if (all->rt_opt.verbose) {
|
|
FIST_Warning("DetermineNormal() - zero-length normal vector!\n");
|
|
}
|
|
normal->x = normal->y = C_0_0;
|
|
normal->z = C_1_0;
|
|
return false;
|
|
}
|
|
#endif
|
|
}
|
|
|
|
|
|
/* */
|
|
/* this function maps the vertices of the polygons referenced by */
|
|
/* loops[i1,..,i2-1] to the plane n3.x * x + n3.y * y + n3.z * z = 0. */
|
|
/* every mapped vertex (x,y,z) is then expressed in terms of (x',y',z'), */
|
|
/* where z'=0. this is achieved by transforming the original vertices into */
|
|
/* a coordinate system whose z-axis coincides with n3, and whose two other */
|
|
/* coordinate axes n1 and n2 are orthonormal on n3. note that n3 is */
|
|
/* supposed to be of unit length! */
|
|
/* */
|
|
void ProjectPoints(global_struct *all, int i1, int i2, vertex n3)
|
|
{
|
|
listdef *list = &all->c_list;
|
|
vertexdef *vert = &all->c_vertex;
|
|
datadef *data = &all->c_data;
|
|
|
|
Matrix4 matrix;
|
|
vertex vtx, n1, n2;
|
|
double d;
|
|
list_ind ind, ind1;
|
|
int i, j1;
|
|
|
|
/* */
|
|
/* choose n1 and n2 appropriately */
|
|
/* */
|
|
if ((Abs(n3.x) > C_0_1) || (Abs(n3.y) > C_0_1)) {
|
|
n1.x = -n3.y;
|
|
n1.y = n3.x;
|
|
n1.z = C_0_0;
|
|
}
|
|
else {
|
|
n1.x = n3.z;
|
|
n1.z = -n3.x;
|
|
n1.y = C_0_0;
|
|
}
|
|
d = Length_l2(n1);
|
|
DivScalar(d, n1);
|
|
VectorProduct(n1, n3, n2);
|
|
d = Length_l2(n2);
|
|
DivScalar(d, n2);
|
|
|
|
/* */
|
|
/* initialize the transformation matrix */
|
|
/* */
|
|
matrix.element[0][0] = n1.x;
|
|
matrix.element[1][0] = n1.y;
|
|
matrix.element[2][0] = n1.z;
|
|
matrix.element[3][0] = C_0_0; /* translation of the coordinate system */
|
|
matrix.element[0][1] = n2.x;
|
|
matrix.element[1][1] = n2.y;
|
|
matrix.element[2][1] = n2.z;
|
|
matrix.element[3][1] = C_0_0; /* translation of the coordinate system */
|
|
matrix.element[0][2] = n3.x;
|
|
matrix.element[1][2] = n3.y;
|
|
matrix.element[2][2] = n3.z;
|
|
matrix.element[3][2] = C_0_0; /* translation of the coordinate system */
|
|
matrix.element[0][3] = C_0_0;
|
|
matrix.element[1][3] = C_0_0;
|
|
matrix.element[2][3] = C_0_0;
|
|
matrix.element[3][3] = C_1_0;
|
|
|
|
/* */
|
|
/* transform the vertices and store the transformed vertices in the array */
|
|
/* `points' */
|
|
/* */
|
|
|
|
// MODIF : alloco solo lo spazio necessario per l'array dei punti ( verrà eventualmente esteso da StorePnt)
|
|
// InitPnts(data, BLOCK_SIZE);
|
|
InitPnts( data, data->max_num_pnts) ;
|
|
for (i = i1; i < i2; ++i) {
|
|
ind = list->loops[i];
|
|
ind1 = ind;
|
|
j1 = GetIndex(list, ind1);
|
|
MatrixTransformation(matrix, vert->vertices[j1], vtx);
|
|
#ifdef EXT_APPL_SITES
|
|
j1 = StorePnt(data, vtx.x, vtx.y, GetExtApplVtx(j1));
|
|
#else
|
|
j1 = StorePnt(data, vtx.x, vtx.y);
|
|
#endif
|
|
UpdateIndex(list, ind1, j1);
|
|
ind1 = GetNextNode(list, ind1);
|
|
j1 = GetIndex(list, ind1);
|
|
while (ind1 != ind) {
|
|
MatrixTransformation(matrix, vert->vertices[j1], vtx);
|
|
#ifdef EXT_APPL_SITES
|
|
j1 = StorePnt(data, vtx.x, vtx.y, GetExtApplVtx(j1));
|
|
#else
|
|
j1 = StorePnt(data, vtx.x, vtx.y);
|
|
#endif
|
|
UpdateIndex(list, ind1, j1);
|
|
ind1 = GetNextNode(list, ind1);
|
|
j1 = GetIndex(list, ind1);
|
|
}
|
|
}
|
|
|
|
return;
|
|
}
|