31 static int comp_double(
double *,
double *);
32 static int V__within(
double,
double,
double);
35 static void destroy_links(
struct link_head *,
struct Slink *);
36 static int Vect__divide_and_conquer(
struct Slink *,
const struct line_pnts *,
60 static int first_time = 1;
61 static int isl_allocated = 0;
64 G_debug(3,
"Vect_get_point_in_area()");
72 if (n_isles > isl_allocated) {
75 for (i = isl_allocated; i < n_isles; i++)
77 isl_allocated = n_isles;
83 for (i = 0; i < n_isles; i++) {
91 (
const struct line_pnts**) IPoints, n_isles, X, Y));
96 static int comp_double(
double *i,
double *j)
104 static int V__within(
double a,
double x,
double b)
107 return (
x >= a &&
x < b);
109 return (
x > b &&
x <= a);
132 double a,
b, c, d,
x;
135 for (i = 1; i < Points->
n_points; i++) {
136 a = Points->
y[i - 1];
139 c = Points->
x[i - 1];
142 if (V__within(a, y, b)) {
146 perc = (y - a) / (b - a);
147 x = perc * (d - c) + c;
176 double a,
b, c, d,
y;
179 for (i = 1; i < Points->
n_points; i++) {
180 a = Points->
x[i - 1];
183 c = Points->
y[i - 1];
186 if (V__within(a, x, b)) {
190 perc = (x - a) / (b - a);
191 y = perc * (d - c) + c;
213 double cent_x, cent_y;
217 static int first_time = 1;
232 G_debug(3,
"Vect_get_point_in_poly(): divide and conquer");
235 x_max = x_min = Points->
x[0];
236 for (i = 0; i < Points->
n_points; i++) {
237 if (x_min > Points->
x[i])
238 x_min = Points->
x[i];
239 if (x_max < Points->
x[i])
240 x_max = Points->
x[i];
252 Head = (
struct Slink *)
link_new(Token);
253 tmp = (
struct Slink *)
link_new(Token);
262 ret = Vect__divide_and_conquer(Head, Points, Token, X, Y, 10);
264 destroy_links(Token, Head);
267 G_warning(
"Vect_get_point_in_poly(): %s",
268 _(
"Unable to find point in polygon"));
272 G_debug(3,
"Found point in %d iterations", 10 - ret);
300 Vect__divide_and_conquer(
struct Slink *Head,
303 double *
X,
double *
Y,
int levels)
305 struct Slink *A, *B, *C;
307 G_debug(3,
"Vect__divide_and_conquer(): LEVEL %d", levels);
312 C = (
struct Slink *)
link_new(Token);
316 C->x = (A->x + B->x) / 2.;
336 return Vect__divide_and_conquer(Head, Points, Token,
X,
Y, --levels);
339 static void destroy_links(
struct link_head *Token,
struct Slink *Head)
341 struct Slink *p, *tmp;
364 double *cent_x,
double *cent_y)
367 double *xptr1, *yptr1;
368 double *xptr2, *yptr2;
369 double cent_weight_x, cent_weight_y;
378 xptr2 = points->
x + 1;
379 yptr2 = points->
y + 1;
382 for (i = 1; i < points->
n_points; i++) {
383 len = hypot(*xptr1 - *xptr2, *yptr1 - *yptr2);
384 cent_weight_x += len * ((*xptr1 + *xptr2) / 2.);
385 cent_weight_y += len * ((*yptr1 + *yptr2) / 2.);
396 *cent_x = cent_weight_x / tot_len;
397 *cent_y = cent_weight_y / tot_len;
458 const struct line_pnts **IPoints,
int n_isles,
459 double *att_x,
double *att_y)
462 static int first_time = 1;
463 double cent_x, cent_y;
465 double max, hi_x, lo_x, hi_y, lo_y;
469 int point_in_sles = 0;
473 G_debug(3,
"Vect_get_point_in_poly_isl(): n_isles = %d", n_isles);
482 *att_x = Points->
x[0];
483 *att_y = Points->
y[0];
494 for (i = 0; i < n_isles; i++) {
500 if (!point_in_sles) {
520 for (i = 0; i < Points->
n_points; i++) {
521 if ((lo_y < cent_y) && (hi_y >= cent_y) &&
522 (lo_x < cent_x) && (hi_x >= cent_x))
524 if (Points->
y[i] < cent_y)
526 if (Points->
y[i] >= cent_y)
529 if (Points->
x[i] < cent_x)
531 if (Points->
x[i] >= cent_x)
535 for (i = 0; i < Points->
n_points; i++) {
536 if ((Points->
y[i] < cent_y) &&
537 ((cent_y - Points->
y[i]) < (cent_y - lo_y)))
539 if ((Points->
y[i] >= cent_y) &&
540 ((Points->
y[i] - cent_y) < (hi_y - cent_y)))
543 if ((Points->
x[i] < cent_x) &&
544 ((cent_x - Points->
x[i]) < (cent_x - lo_x)))
546 if ((Points->
x[i] >= cent_x) &&
547 ((Points->
x[i] - cent_x) < (hi_x - cent_x)))
550 for (i = 0; i < n_isles; i++) {
551 for (j = 0; j < IPoints[i]->
n_points; j++) {
552 if ((IPoints[i]->
y[j] < cent_y) &&
553 ((cent_y - IPoints[i]->
y[j]) < (cent_y - lo_y)))
554 lo_y = IPoints[i]->
y[j];
555 if ((IPoints[i]->
y[j] >= cent_y) &&
556 ((IPoints[i]->
y[j] - cent_y) < (hi_y - cent_y)))
557 hi_y = IPoints[i]->
y[j];
559 if ((IPoints[i]->
x[j] < cent_x) &&
560 ((cent_x - IPoints[i]->
x[j]) < (cent_x - lo_x)))
561 lo_x = IPoints[i]->
x[j];
562 if ((IPoints[i]->
x[j] >= cent_x) &&
563 ((IPoints[i]->
x[j] - cent_x) < (hi_x - cent_x)))
564 hi_x = IPoints[i]->
x[j];
571 *att_y = (hi_y + lo_y) / 2.0;
577 for (i = 0; i < n_isles; i++) {
586 qsort(Intersects->
x, (
size_t) Intersects->
n_points,
sizeof(
double),
587 (
void *)comp_double);
593 for (i = 0; i < Intersects->
n_points; i += 2) {
594 diff = Intersects->
x[i + 1] - Intersects->
x[i];
604 fa = fabs(Intersects->
x[maxpos]);
605 fb = fabs(Intersects->
x[maxpos + 1]);
607 dmax = frexp(fa, &exp);
609 dmax = frexp(fb, &exp);
611 dmax = ldexp(dmax, exp);
614 *att_x = (Intersects->
x[maxpos] + Intersects->
x[maxpos + 1]) / 2.;
618 G_debug(3,
"Vect_get_point_in_poly_isl(): trying x intersect");
623 *att_x = (hi_x + lo_x) / 2.0;
629 for (i = 0; i < n_isles; i++) {
638 qsort(Intersects->
y, (
size_t) Intersects->
n_points,
sizeof(
double),
639 (
void *)comp_double);
645 for (i = 0; i < Intersects->
n_points; i += 2) {
646 diff = Intersects->
y[i + 1] - Intersects->
y[i];
654 fa = fabs(Intersects->
y[maxpos]);
655 fb = fabs(Intersects->
y[maxpos + 1]);
657 dmax = frexp(fa, &exp);
659 dmax = frexp(fb, &exp);
661 dmax = ldexp(dmax, exp);
663 *att_y = (Intersects->
y[maxpos] + Intersects->
y[maxpos + 1]) / 2.;
667 G_warning(
"Vect_get_point_in_poly_isl(): collapsed area");
680 G_warning(
"Vect_get_point_in_poly_isl(), the hard way: centroid is on outer ring, max dist is %g", max);
685 for (i = 0; i < n_isles; i++) {
688 G_warning(
"Vect_get_point_in_poly_isl(), the hard way: centroid is in isle, max dist is %g", max);
692 if (!point_in_sles) {
705 static int segments_x_ray(
double X,
double Y,
const struct line_pnts *Points)
707 double x1, x2, y1, y2;
712 G_debug(3,
"segments_x_ray(): x = %f y = %f n_points = %d",
X,
Y,
719 for (n = 1; n < Points->
n_points; n++) {
720 x1 = Points->
x[n - 1];
721 y1 = Points->
y[n - 1];
738 if (y1 >
Y && y2 >
Y)
742 if (y1 <
Y && y2 <
Y)
746 if (x1 <
X && x2 <
X)
750 if ((x1 ==
X && y1 ==
Y) || (x2 ==
X && y2 ==
Y))
754 if (x1 == x2 && x1 ==
X) {
755 if ((y1 <= Y && y2 >=
Y) || (y1 >=
Y && y2 <=
Y))
760 if (y1 == y2 && y1 ==
Y) {
761 if ((x1 <= X && x2 >=
X) || (x1 >=
X && x2 <=
X))
772 if ((y1 ==
Y && y2 >
Y) || (y2 ==
Y && y1 >
Y))
778 if (y1 ==
Y && y2 <
Y) {
783 if (y2 ==
Y && y1 <
Y) {
790 if ((y1 < Y && y2 >
Y) || (y1 > Y && y2 < Y)) {
791 if (x1 >=
X && x2 >=
X) {
798 G_debug(3,
"x_inter = %f", x_inter);
801 else if (x_inter >
X)
809 (
"segments_x_ray() %s: X = %f Y = %f x1 = %f y1 = %f x2 = %f y2 = %f",
810 _(
"conditions failed"),
X, Y, x1, y1, x2, y2);
830 G_debug(3,
"Vect_point_in_poly(): x = %f y = %f n_points = %d", X, Y,
833 n_intersects = segments_x_ray(X, Y, Points);
835 if (n_intersects == -1)
840 return (n_intersects & 1);
859 static int first = 1;
860 int n_intersects, inter;
868 G_debug(3,
"Vect_point_in_area_outer_ring(): x = %f y = %f area = %d",
877 Area = Plus->
Area[area];
880 if (X < box->
W || X > box->
E || Y > box->
N || Y < box->S)
884 for (i = 0; i < Area->
n_lines; i++) {
885 line = abs(Area->
lines[i]);
902 inter = segments_x_ray(X, Y, Points);
905 n_intersects += inter;
910 return (n_intersects & 1);
928 static int first = 1;
929 int n_intersects, inter;
937 G_debug(3,
"Vect_point_in_island(): x = %f y = %f isle = %d",
946 Isle = Plus->
Isle[isle];
949 if (X < box->
W || X > box->
E || Y > box->
N || Y < box->S)
953 for (i = 0; i < Isle->
n_lines; i++) {
954 line = abs(Isle->
lines[i]);
971 inter = segments_x_ray(X, Y, Points);
974 n_intersects += inter;
979 return (n_intersects & 1);
int Vect_point_in_poly(double X, double Y, const struct line_pnts *Points)
Determines if a point (X,Y) is inside a polygon.
int Vect_get_area_isle(const struct Map_info *, int, int)
Returns isle id for area.
int alloc_points
Allocated space for points.
int Vect__intersect_x_line_with_poly()
struct P_area ** Area
Array of areas.
int n_points
Number of points.
void link_dispose(struct link_head *, VOID_T *)
int Vect_get_area_num_isles(const struct Map_info *, int)
Returns number of isles for given area.
struct P_isle ** Isle
Array of isles.
plus_t n_lines
Number of boundary lines.
int Vect_point_in_area_outer_ring(double X, double Y, const struct Map_info *Map, int area, struct bound_box *box)
Determines if a point (X,Y) is inside an area outer ring. Islands are not considered.
double dig_x_intersect(double, double, double, double, double)
double * x
Array of X coordinates.
int Vect_get_point_in_area(const struct Map_info *Map, int area, double *X, double *Y)
Get point inside area and outside all islands.
Feature geometry info - coordinates.
int Vect_point_in_island(double X, double Y, const struct Map_info *Map, int isle, struct bound_box *box)
Determines if a point (X,Y) is inside an island.
VOID_T * link_new(struct link_head *)
plus_t * lines
List of boundary lines.
void link_exit_on_error(int)
Basic topology-related info.
struct line_pnts * Vect_new_line_struct(void)
Creates and initializes a line_pnts structure.
int Vect_get_isle_points(const struct Map_info *, int, struct line_pnts *)
Returns polygon array of points for given isle.
int Vect_append_point(struct line_pnts *, double, double, double)
Appends one point to the end of a line.
plus_t * lines
List of boundary lines.
struct Plus_head plus
Plus info (topology, version, ...)
int Vect_get_point_in_poly(const struct line_pnts *Points, double *X, double *Y)
Get point inside polygon.
int Vect_find_poly_centroid(const struct line_pnts *points, double *cent_x, double *cent_y)
Get centroid of polygon.
struct link_head * link_init(int)
int Vect_get_point_in_poly_isl(const struct line_pnts *Points, const struct line_pnts **IPoints, int n_isles, double *att_x, double *att_y)
Get point inside polygon but outside the islands specifiled in IPoints.
double * y
Array of Y coordinates.
void G_warning(const char *,...) __attribute__((format(printf
int Vect__intersect_y_line_with_poly()
int Vect_get_area_points(const struct Map_info *, int, struct line_pnts *)
Returns polygon array of points (outer ring) of given area.
plus_t n_lines
Number of boundary lines.
int Vect_read_line(const struct Map_info *, struct line_pnts *, struct line_cats *, int)
Read vector feature (topological level required)
int G_debug(int, const char *,...) __attribute__((format(printf