Commit 1b56f6b9 authored by Célestin Marot's avatar Célestin Marot
Browse files

need to revert back to old predicates

parent 6e889870
......@@ -23,8 +23,6 @@ Author: Célestin Marot (celestin.marot@uclouvain.be) */
#include <math.h>
#include <assert.h>
// #define SPR_RECURSIVE_VERSION
#ifdef _MSC_VER
#define SPRNOINLINE __declspec(noinline)
......@@ -159,8 +157,7 @@ static SPRNOINLINE uint32_t add_orient3d(SPRCavity* SPR,
const unsigned v[4] = {v0, v1, v2, v3};
// compute the orientation
double orientation = orient3d(SPR->values,
SPR->points.array[v0].coord,
double orientation = orient3d(SPR->points.array[v0].coord,
SPR->points.array[v1].coord,
SPR->points.array[v2].coord,
SPR->points.array[v3].coord);
......@@ -226,8 +223,7 @@ static int get_orient3d(SPRCavity* SPR,
if(v0==v1 || v0==v2 || v0==v3 || v1==v2 || v1==v3 || v2==v3)
return 0;
return orient3d_inline(SPR->values,
SPR->points.array[v0].coord,
return orient3d_inline(SPR->points.array[v0].coord,
SPR->points.array[v1].coord,
SPR->points.array[v2].coord,
SPR->points.array[v3].coord);
......@@ -608,7 +604,6 @@ static int tet_contains_pt(SPRCavity* SPR,
return 1;
}
#ifdef SPR_SAVE_QUALITIES
static SPRNOINLINE double add_quality_map(SPRCavity* SPR,
unsigned v0,
unsigned v1,
......@@ -694,38 +689,6 @@ static double get_quality_map(SPRCavity* SPR,
}
}
#else // <-- SPR_SAVE_QUALITIES
static double get_quality_map(SPRCavity* SPR,
unsigned v0,
unsigned v1,
unsigned v2,
unsigned v3)
{
int orientation = get_orient3d(SPR, v0, v1, v2, v3);
if(orientation > 0)
return -1.0;
else if(orientation==0)
return 0.0;
else if(SPR->quality.function==NULL)
return 1.0;
double* c0 = SPR->points.array[v0].coord;
double* c1 = SPR->points.array[v1].coord;
double* c2 = SPR->points.array[v2].coord;
double* c3 = SPR->points.array[v3].coord;
double qual = SPR->quality.function(c0, c1, c2, c3,
SPR->quality.userData);
if(qual<0.0) {
fprintf(stderr, "WARNING: negative quality with correct orientation\n");
qual = DBL_MIN;
}
return qual;
}
#endif // <-- SPR_SAVE_QUALITIES
static void compute_bbox(SPRCavity* SPR, uint16_t triangleID)
{
......@@ -745,8 +708,6 @@ static inline uint16_t best_face_heuristic(SPRCavity* SPR, double* qualMax)
{
const uint16_t nfaces = SPR->triangles.num;
#ifdef SPR_USE_FACE_HEURISTIC
const uint8_t npts = SPR->points.num;
int minValid = 255;
uint16_t worstFaceID = UINT16_MAX;
......@@ -796,10 +757,6 @@ static inline uint16_t best_face_heuristic(SPRCavity* SPR, double* qualMax)
*qualMax = minimax;
return worstFaceID;
#else // <-- SPR_USE_FACE_HEURISTIC
*qualMax = 1.0;
return nfaces-1;
#endif // <-- SPR_USE_FACE_HEURISTIC
}
......@@ -816,10 +773,8 @@ static inline void compute_candidates(SPRCavity* SPR,
SPRBbox* triangleBbox = &SPR->triangles.bbox[triangleID];
uint8_t num = 0;
#ifdef SPR_SORT_CANDIDATES
uint8_t score[SPR_MAX_PTS];
uint8_t candidate[SPR_MAX_PTS];
#endif
for (tet->node[3]=0; tet->node[3]<nPts; tet->node[3]++)
{
......@@ -840,7 +795,6 @@ static inline void compute_candidates(SPRCavity* SPR,
SPRBbox tetBbox = *triangleBbox;
SPRBbox_add(&tetBbox, SPR->points.array[tet->node[3]].coord);
#ifndef SPR_SAVE_QUALITIES
for (int ptID=0; ptID<nPts; ptID++) {
if(tet_contains_pt(SPR, &tetBbox, tet->node, ptID)) {
stop = 1;
......@@ -855,7 +809,6 @@ static inline void compute_candidates(SPRCavity* SPR,
if(stop)
continue;
#endif
// test intersection with boundary triangles
for (int otherFaceID=0; otherFaceID<nfaces; otherFaceID++)
......@@ -869,7 +822,6 @@ static inline void compute_candidates(SPRCavity* SPR,
if(stop)
continue;
#ifdef SPR_SORT_CANDIDATES
// those heuristics can be changed in many way
// if you have some theory validating some choice,
// you can give a higher/lower score depending on
......@@ -903,18 +855,13 @@ static inline void compute_candidates(SPRCavity* SPR,
curScore+=2;
}
// the maximum score is 4+2+2+3*2 = 15
// the maximum score is 4+2+2+3*2 = 14
candidate[num] = tet->node[3];
score[num] = curScore;
num++;
#else
step->candidate[num] = tet->node[3];
num++;
#endif
}
#ifdef SPR_SORT_CANDIDATES
// simple radix sort pass
uint8_t count[16] = {0};
for (uint8_t i=0; i<num; i++) {
......@@ -933,7 +880,6 @@ static inline void compute_candidates(SPRCavity* SPR,
for (uint8_t i=0; i<num; i++) {
step->candidate[count[score[i]]++] = candidate[i];
}
#endif // <-- SPR_SORT_CANDIDATES
step->numCandidates=num;
}
......@@ -1009,18 +955,8 @@ static inline int SPR_init(SPRCavity* SPR)
assert(SPR->triangles.num <= SPR_MAX_FACES);
const size_t n = SPR->points.num;
{ // allocatinc steps
size_t max_tets = ((n*n - 3*n - 2)/2);
SPR->steps.array = malloc(sizeof(SPRStep) * max_tets);
if(SPR->steps.array==NULL)
goto allocation_error;
}
{ // allocating & initializing map.face
size_t face_len = SPR_MAX_PTS*SPR_MAX_PTS*SPR_MAX_PTS;
SPR->map.faces = malloc(sizeof(uint16_t) * face_len);
if(SPR->map.faces==NULL)
goto allocation_error;
{ // initializing map.face
size_t face_len = sizeof(SPR->map.faces)/sizeof(uint16_t);
// compiler should replace this by a memset (UINT16_MAX is all ones)
for (size_t i=0; i<face_len; i++) {
......@@ -1028,44 +964,31 @@ static inline int SPR_init(SPRCavity* SPR)
}
}
#ifdef SPR_SAVE_QUALITIES
{ // allocating & initializing map.qualities
size_t qual_len = n*(n-1)*(n-2)*(n-3)/24;
SPR->map.qualities = malloc(sizeof(double)*qual_len);
if(SPR->map.qualities==NULL)
goto allocation_error;
{ // initializing map.qualities
size_t qual_len = sizeof(SPR->map.qualities)/sizeof(double);
// compiler should also replace this by a memset (NAN is all ones)
for (size_t i=0; i<qual_len; i++) {
SPR->map.qualities[i] = NAN;
}
}
#endif
#ifdef SPR_SAVE_ORIENT3D
{ // allocating & setting to zero map.orient3d
size_t chir_len = n*n*n*n/16+1;
SPR->map.orient3d = calloc(chir_len, sizeof(uint32_t));
if(SPR->map.orient3d==NULL)
goto allocation_error;
size_t ori_len = sizeof(SPR->map.orient3d)/sizeof(uint32_t);
// compiler should also replace this by a memset
for (size_t i=0; i<ori_len; i++) {
SPR->map.orient3d[i] = 0;
}
}
#endif
SPRBbox bbox = {{ DBL_MAX, DBL_MAX, DBL_MAX},
{-DBL_MAX,-DBL_MAX,-DBL_MAX}};
for (int i=0; i<SPR->points.num; i++) {
SPRPoint* p = &SPR->points.array[i];
p->valence = 0;
SPRBbox_add(&bbox, p->coord);
}
exactinit(SPR->values,
bbox.max[0] - bbox.min[0],
bbox.max[1] - bbox.min[1],
bbox.max[2] - bbox.min[2]);
for (uint16_t i=0; i<SPR->triangles.num; i++) {
uint8_t* v = SPR->triangles.array[i].node;
v[3] = v[0];
......@@ -1085,27 +1008,6 @@ static inline int SPR_init(SPRCavity* SPR)
}
}
return 0;
allocation_error:
fprintf(stderr, "memory allocation failure\n");
return 1;
}
static inline void SPR_terminate(SPRCavity* SPR)
{
free(SPR->steps.array);
SPR->steps.array = NULL;
free(SPR->map.faces);
SPR->map.faces = NULL;
#ifdef SPR_SAVE_QUALITIES
free(SPR->map.qualities);
SPR->map.qualities = NULL;
#endif
#ifdef SPR_SAVE_ORIENT3D
free(SPR->map.orient3d);
SPR->map.orient3d = NULL;
#endif
}
......@@ -1129,12 +1031,8 @@ static inline void SPR_terminate(SPRCavity* SPR)
|*| | | \ \ / -_) _/ _ \ ' \| ' \/ -_) _| _| / _ \ ' \ |*|
|*| |_| \_\\___\__\___/_||_|_||_\___\__|\__|_\___/_||_| |*|
\*|=======================================================|*/
#ifndef SPR_RECURSIVE_VERSION
int hxtSPR(SPRCavity* SPR)
{
SPRExact valuesOnStack; // structure needed for the exact predicates to work
SPR->values = &valuesOnStack;
if(SPR_init(SPR))
return 2;
......@@ -1145,12 +1043,9 @@ int hxtSPR(SPRCavity* SPR)
while(1) {
if(SPR->num_search_nodes >= SPR->max_search_nodes){
SPR_terminate(SPR);
return 1;
}
SPR->num_search_nodes++;
int rewind = 0;
if(SPR->triangles.num==0) { // we found a new (better) tetrahedralization
new_tetrahedralization(SPR, step);
......@@ -1190,6 +1085,8 @@ int hxtSPR(SPRCavity* SPR)
if(qual <= SPR->tetrahedra.quality)
continue;
SPR->num_search_nodes++;
for (int f=0; f<3; f++) {
uint8_t p0 = step->tet.node[nodeFromFacet[0][f]];
uint8_t p1 = step->tet.node[nodeFromFacet[1][f]];
......@@ -1222,7 +1119,6 @@ int hxtSPR(SPRCavity* SPR)
if(rewind) {
if(SPR->steps.num==0) {
SPR_terminate(SPR);
return 0;
}
......@@ -1249,129 +1145,4 @@ int hxtSPR(SPRCavity* SPR)
}
}
}
}
#else
/*****************************************************************
* Here below is a recursive version. It works and is more easy *
* to understand (IMO) than the iterative one *
*****************************************************************/
void hxtSPR_recur(SPRCavity* SPR) {
if(SPR->num_search_nodes >= SPR->max_search_nodes){
return;
}
SPR->num_search_nodes++;
SPRStep* step = &SPR->steps.array[SPR->steps.num];
if(SPR->triangles.num==0) { // we found a new (better) tetrahedralization
new_tetrahedralization(SPR, step);
return;
}
{
double qualMax;
const uint16_t triangleID = best_face_heuristic(SPR, &qualMax);
if(triangleID == UINT16_MAX)
return;
step->tet = SPR->triangles.array[triangleID];
// create the candidate array
compute_candidates(SPR, step, qualMax, triangleID);
// remove the face from the map
remove_face(SPR, triangleID);
}
for (; step->numCandidates; step->numCandidates--)
{
step->tet.node[3] = step->candidate[step->numCandidates-1];
double qual = get_quality_map(SPR,
step->tet.node[0],
step->tet.node[1],
step->tet.node[2],
step->tet.node[3]);
if(qual <= SPR->tetrahedra.quality)
continue;
for (int f=0; f<3; f++) {
uint8_t p0 = step->tet.node[nodeFromFacet[0][f]];
uint8_t p1 = step->tet.node[nodeFromFacet[1][f]];
uint8_t p2 = step->tet.node[nodeFromFacet[2][f]];
uint16_t index = get_face_map(SPR, p0, p1, p2);
if(index==UINT16_MAX) {
// add the face to the map
add_face(SPR, p0, p2, p1);
}
else {
// remove the face from the map
remove_face(SPR, index);
}
}
SPRStep* nextStep = &SPR->steps.array[SPR->steps.num+1];
nextStep->quality = qual < step->quality ? qual : step->quality;
nextStep->numCandidates = -1;
SPR->steps.num++;
hxtSPR_recur(SPR);
SPR->steps.num--;
for (int f=0; f<3; f++) {
uint8_t p0 = step->tet.node[nodeFromFacet[0][f]];
uint8_t p1 = step->tet.node[nodeFromFacet[1][f]];
uint8_t p2 = step->tet.node[nodeFromFacet[2][f]];
uint16_t index = get_face_map(SPR, p0, p2, p1);
if(index==UINT16_MAX) {
// add the face to the map
add_face(SPR, p0, p1, p2);
}
else {
// remove the face from the map
remove_face(SPR, index);
}
}
if(step->quality <= SPR->tetrahedra.quality ||
SPR->num_search_nodes >= SPR->max_search_nodes)
break;
}
add_face(SPR,
step->tet.node[0],
step->tet.node[1],
step->tet.node[2]);
}
int hxtSPR(SPRCavity* SPR) {
SPRExact valuesOnStack; // structure needed for the exact predicates to work
SPR->values = &valuesOnStack;
if(SPR_init(SPR))
return 2;
SPR->steps.num = 0;
SPR->steps.array[0].quality = DBL_MAX; // we can hope for the biggest quality at the moment
hxtSPR_recur(SPR);
SPR_terminate(SPR);
if(SPR->num_search_nodes>=SPR->max_search_nodes)
return 1;
return 0;
}
#endif
\ No newline at end of file
}
\ No newline at end of file
......@@ -29,28 +29,22 @@ extern "C" {
#include <stdio.h>
#include <stdint.h>
/* /!\ Warning /!\ This version is a modified version of HXTSPR
*
* One must call exactinit() correctly before
* being able to use any function defined here
*/
// defines to set/unset the main optimizations
#define SPR_SORT_CANDIDATES
#define SPR_SAVE_QUALITIES
#define SPR_USE_FACE_HEURISTIC
// #define SPR_SAVE_ORIENT3D
// #define SPR_INPUT_VERIFICATION // TODO: not implemented yet
#define SPR_SAVE_ORIENT3D
#define SPR_MAX_PTS 48 // you can change only this value
#define SPR_MAX_PTS 32
#define SPR_MAX_EDGES ((SPR_MAX_PTS*SPR_MAX_PTS - SPR_MAX_PTS)/2)
#define SPR_MAX_TETS ((SPR_MAX_PTS*SPR_MAX_PTS - 3*SPR_MAX_PTS - 2)/2)
#define SPR_MAX_FACES (3*SPR_MAX_TETS)
#if SPR_MAX_PTS > 99
#warning "Using SPR with cavities containing more" \
"than a hundred points might be prohibitively slow"
#if SPR_MAX_PTS > 255
#error "Points are indexed with uint8_t, their number is thus limited to 255."
#endif
#endif
/* structure for the points (32 bytes)*/
typedef struct {
double coord[3]; // the coordinates of the vertex (unchanged by hxtSPR)
......@@ -162,21 +156,21 @@ typedef struct {
/* up to this point, you should not fill the structure yourself
* these are all internal stuffs */
struct {
SPRStep* array;
SPRStep array[SPR_MAX_TETS];
int num;
} steps;
struct {
#ifdef SPR_SAVE_QUALITIES
double* qualities;
#endif
double qualities[SPR_MAX_PTS*
(SPR_MAX_PTS-1)*
(SPR_MAX_PTS-2)*
(SPR_MAX_PTS-3)
/24];
uint16_t faces[SPR_MAX_PTS*SPR_MAX_PTS*SPR_MAX_PTS];
#ifdef SPR_SAVE_ORIENT3D
uint32_t* orient3d;
uint32_t orient3d[SPR_MAX_PTS*SPR_MAX_PTS*SPR_MAX_PTS*SPR_MAX_PTS/16+1];
#endif
uint16_t* faces;
} map;
SPRExact* values;
} SPRCavity;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment