Skip to content

Commit

Permalink
fixed a bug in convex hull calculation, closes #805
Browse files Browse the repository at this point in the history
  • Loading branch information
@ntamas
ntamas committed Jan 9, 2015
1 parent 5a9f51c commit c4b53ab
Show file tree
Hide file tree
Showing 3 changed files with 152 additions and 65 deletions.
120 changes: 104 additions & 16 deletions examples/simple/igraph_convex_hull.c
View file
Expand Up @@ -23,22 +23,14 @@

#include <igraph.h>

int main() {
igraph_real_t coords_array[][2] =
{{3, 2}, {5, 1}, {4, 4}, {6, 4}, {4, 3},
{2, 5}, {1, 3}, {2, 4}, {6, 3}, {9, 2}
};

igraph_matrix_t coords, resmat;
int check_convex_hull(igraph_matrix_t* coords) {
igraph_vector_t result;
long i;

igraph_matrix_init(&coords, 10, 2);
for (i=0; i<20; i++) MATRIX(coords, i/2, i%2) = coords_array[i/2][i%2];

igraph_matrix_t resmat;
long int i;

/* Testing with index output mode */
igraph_vector_init(&result, 1);
if (igraph_convex_hull(&coords, &result, 0))
if (igraph_convex_hull(coords, &result, 0))
return 1;

for (i=0; i<igraph_vector_size(&result); i++)
Expand All @@ -48,15 +40,111 @@ int main() {

/* Testing with coordinate output mode */
igraph_matrix_init(&resmat, 0, 0);
if (igraph_convex_hull(&coords, 0, &resmat))
if (igraph_convex_hull(coords, 0, &resmat))
return 1;

for (i=0; i<igraph_matrix_nrow(&resmat); i++)
printf("%ld %ld ", (long)MATRIX(resmat, i, 0), (long)MATRIX(resmat, i, 1));
printf("%.3f %.3f ", MATRIX(resmat, i, 0), MATRIX(resmat, i, 1));
printf("\n");

igraph_matrix_destroy(&resmat);

return 0;
}

int test_simple() {
igraph_real_t coords_array[][2] =
{{3, 2}, {5, 1}, {4, 4}, {6, 4}, {4, 3},
{2, 5}, {1, 3}, {2, 4}, {6, 3}, {9, 2}
};
igraph_matrix_t coords;
int i, result;

printf("test_simple\n");

igraph_matrix_init(&coords, 10, 2);
for (i=0; i<20; i++) MATRIX(coords, i/2, i%2) = coords_array[i/2][i%2];
result = check_convex_hull(&coords);
igraph_matrix_destroy(&coords);

return result;
}

int test_collinear() {
igraph_real_t coords_array[][2] =
{{3, 2}, {5, 1}, {7, 0}, {9, -1}, {11, -2}};
igraph_matrix_t coords;
int i, result;

printf("test_collinear\n");

igraph_matrix_init(&coords, 5, 2);
for (i=0; i<10; i++) MATRIX(coords, i/2, i%2) = coords_array[i/2][i%2];
result = check_convex_hull(&coords);
igraph_matrix_destroy(&coords);
return result;
}

int test_degenerate() {
igraph_matrix_t coords;
int i, result;

printf("test_degenerate\n");

igraph_matrix_init(&coords, 2, 2);
MATRIX(coords, 0, 0) = 3; MATRIX(coords, 0, 1) = 2;
MATRIX(coords, 1, 0) = 5; MATRIX(coords, 1, 1) = 1;
result = check_convex_hull(&coords);

igraph_matrix_resize(&coords, 1, 2);
MATRIX(coords, 0, 0) = 3; MATRIX(coords, 0, 1) = 2;
result = check_convex_hull(&coords);

igraph_matrix_resize(&coords, 0, 2);
result = check_convex_hull(&coords);

igraph_matrix_destroy(&coords);
return result;
}

int test_bug_805() {
igraph_real_t coords_array[][2] = {
{0, 0}, {1, 0}, {0.707, 0.707}, {0, 1}, {-0.707, 0.707}, {-1, 0},
{-0.707, -0.707}, {0, -1}, {0.707, -0.707}, {2, 0}, {1.414, 1.414}, {0, 2},
{-1.414, 1.414}, {-2, 0}, {-1.414, -1.414}, {0, -2}, {1.414, -1.414}, {3, 0},
{2.121, 2.121}, {0, 3}, {-2.121, 2.121}, {-3, 0}, {-2.121, -2.121}, {0, -3},
{2.121, -2.121}, {4, 0}, {2.828, 2.828}, {0, 4}, {-2.828, 2.828}, {-4, 0},
{-2.828, -2.828}, {0, -4}, {2.828, -2.828}
};
igraph_matrix_t coords;
int i, result;

printf("test_bug_805\n");

igraph_matrix_init(&coords, 33, 2);
for (i=0; i<66; i++) MATRIX(coords, i/2, i%2) = coords_array[i/2][i%2];
result = check_convex_hull(&coords);
igraph_matrix_destroy(&coords);
return 0;
}

int main() {
int result;

result = test_simple();
if (result != 0)
return result;

result = test_collinear();
if (result != 0)
return result;

result = test_degenerate();
if (result != 0)
return result;

result = test_bug_805();
if (result != 0)
return result;

return 0;
}
16 changes: 15 additions & 1 deletion examples/simple/igraph_convex_hull.out
View file
@@ -1,2 +1,16 @@
test_simple
1 6 5 3 9
5 1 1 3 2 5 6 4 9 2
5.000 1.000 1.000 3.000 2.000 5.000 6.000 4.000 9.000 2.000
test_collinear
4 0
11.000 -2.000 3.000 2.000
test_degenerate
1 0
5.000 1.000 3.000 2.000
0
3.000 2.000


test_bug_805
31 30 29 28 27 26 25 32
0.000 -4.000 -2.828 -2.828 -4.000 0.000 -2.828 2.828 0.000 4.000 2.828 2.828 4.000 0.000 2.828 -2.828
81 changes: 33 additions & 48 deletions src/other.c
View file
Expand Up @@ -179,55 +179,40 @@ int igraph_convex_hull(const igraph_matrix_t *data, igraph_vector_t *resverts,
igraph_Free(angles);
IGRAPH_FINALLY_CLEAN(1);

if (no_of_nodes == 1) {
IGRAPH_CHECK(igraph_vector_push_back(&stack, 0));
igraph_indheap_delete_max(&order);
} else {
/* Do the trick */
IGRAPH_CHECK(igraph_vector_push_back(&stack, igraph_indheap_max_index(&order)-1));
igraph_indheap_delete_max(&order);
IGRAPH_CHECK(igraph_vector_push_back(&stack, igraph_indheap_max_index(&order)-1));
igraph_indheap_delete_max(&order);

j=2;
while (!igraph_indheap_empty(&order)) {
/* Determine whether we are at a left or right turn */
last_idx=(long int) VECTOR(stack)[j-1];
before_last_idx=(long int) VECTOR(stack)[j-2];
next_idx=(long)igraph_indheap_max_index(&order)-1;
igraph_indheap_delete_max(&order);
j=0;
last_idx=-1;
before_last_idx=-1;
while (!igraph_indheap_empty(&order)) {
next_idx=(long)igraph_indheap_max_index(&order)-1;
/* Determine whether we are at a left or right turn */
if (j < 2) {
/* Pretend that we are turning into the right direction if we have less
* than two items in the stack */
cp=-1;
} else {
cp=(MATRIX(*data, last_idx, 0)-MATRIX(*data, before_last_idx, 0))*
(MATRIX(*data, next_idx, 1)-MATRIX(*data, before_last_idx, 1))-
(MATRIX(*data, next_idx, 0)-MATRIX(*data, before_last_idx, 0))*
(MATRIX(*data, last_idx, 1)-MATRIX(*data, before_last_idx, 1));
/*
printf("B L N cp: %d, %d, %d, %f [", before_last_idx, last_idx, next_idx, (float)cp);
for (k=0; k<j; k++) printf("%ld ", (long)VECTOR(stack)[k]);
printf("]\n");
*/
if (cp == 0) {
/* The last three points are collinear. Replace the last one in
* the stack to the newest one */
VECTOR(stack)[j-1]=next_idx;
} else if (cp < 0) {
/* We are turning into the right direction */
IGRAPH_CHECK(igraph_vector_push_back(&stack, next_idx));
j++;
} else {
/* No, skip back until we're okay */
while (cp >= 0 && j > 2) {
igraph_vector_pop_back(&stack);
j--;
last_idx=(long int) VECTOR(stack)[j-1];
before_last_idx=(long int) VECTOR(stack)[j-2];
cp=(MATRIX(*data, last_idx, 0)-MATRIX(*data, before_last_idx, 0))*
(MATRIX(*data, next_idx, 1)-MATRIX(*data, before_last_idx, 1))-
(MATRIX(*data, next_idx, 0)-MATRIX(*data, before_last_idx, 0))*
(MATRIX(*data, last_idx, 1)-MATRIX(*data, before_last_idx, 1));
}
IGRAPH_CHECK(igraph_vector_push_back(&stack, next_idx));
j++;
}
(MATRIX(*data, next_idx, 1)-MATRIX(*data, before_last_idx, 1))-
(MATRIX(*data, next_idx, 0)-MATRIX(*data, before_last_idx, 0))*
(MATRIX(*data, last_idx, 1)-MATRIX(*data, before_last_idx, 1));
}
/*
printf("B L N cp: %ld, %ld, %ld, %f [", before_last_idx, last_idx, next_idx, (float)cp);
for (int k=0; k<j; k++) printf("%ld ", (long)VECTOR(stack)[k]);
printf("]\n");
*/
if (cp < 0) {
/* We are turning into the right direction */
igraph_indheap_delete_max(&order);
IGRAPH_CHECK(igraph_vector_push_back(&stack, next_idx));
before_last_idx = last_idx;
last_idx = next_idx;
j++;
} else {
/* No, skip back and try again in the next iteration */
igraph_vector_pop_back(&stack);
j--;
last_idx = before_last_idx;
before_last_idx = (j >= 2) ? (long int) VECTOR(stack)[j-2] : -1;
}
}

Expand Down

0 comments on commit c4b53ab

Please sign in to comment.