# Выпуклая оболочка точек на плоскости

• теория?

```#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define EPS 1E-10

int debug = 0;

typedef struct point
{
double x, y;
} point_t;

typedef struct line
{
point_t a, b;
} line_t;

typedef struct ppoly
{
int n;         // number of points in polyline
point_t **ppts;      //array of pointers to points
} ppoly_t;

typedef struct poly
{
int n;         // number of points in polyline
point_t *pts;         //array of points
} poly_t;

/*
Current axes origin for cmp function
*/
point_t *orig;

/*
Calculates sign of VectorProduct [ orig->a, orig->b ].
*/
int
cmp (const void *aa, const void *bb)
{
point_t *a = *(point_t **) aa;
point_t *b = *(point_t **) bb;
double s = (a->x - orig->x) * (b->y - orig->y) - (b->x - orig->x) * (a->y - orig->y);
if (fabs (s) < EPS)
{
s = (a->x - orig->x) * (a->x - orig->x) + (a->y - orig->y) * (a->y - orig->y)
- (b->x - orig->x) * (b->x - orig->x) + (b->y - orig->y) * (b->y - orig->y);
if (fabs (s) < EPS)
return 0;
if (s < 0)
return 1;
return -1;
};
if (s < 0)
return 1;
return -1;
}

/*
Calculates convex hull (polyline) of given array of points
*/
ppoly_t
convex_hull (point_t * pts, int n)
{
int sp;         /* Stack pointer */
int i, j, i_min = 0;
point_t **stack = (point_t **) malloc ((n + 1) * sizeof (point_t *));
ppoly_t result;

for (i = 0; i < n; i++)
{
if (pts[i].y < pts[i_min].y)
i_min = i;
stack[i] = pts + i;
}

orig = stack[i_min];
stack[i_min] = stack[0];
stack[0] = orig;
qsort (stack + 1, n - 1, sizeof (point_t *), cmp);
sp = 2;
i = 2;
while (i < n)
{
orig = stack[sp - 1];
if (cmp (&stack[sp], &stack[(i + 1) % n]) <= 0)
stack[++sp] = stack[(++i) % n];
else
sp--;
}
result.ppts = stack;
result.n = sp;
return result;
}

void
free_poly (poly_t p)
{
if (p.pts)
free (p.pts);
}

void
free_ppoly (ppoly_t p)
{
if (p.ppts)
free (p.ppts);
}

/*
Interesection of line ..--a0--a1--..  and  segment [b0, b1].
Return TRUE if there is the intersection and put intersection point in variable *res;
Otherwise return FALSE;
*/

int
line_segment_intersect (point_t a0, point_t a1, point_t b0, point_t b1,
point_t * res)
{
point_t da, db;
double det = 0;
double t_a = 0.0, t_b = 0.0;
da.x = a1.x - a0.x;
da.y = a1.y - a0.y;
db.x = b1.x - b0.x;
db.y = b1.y - b0.y;

// line 1 = a0 + da*t_a
// line 2 = b0 + db*t_b

// Solve for t_a and t_b.

// Check for divide by zero
det = da.x * db.y - da.y * db.x;
if (fabs (det) < EPS)
return 0;

t_a = db.x * (a0.y - b0.y) - db.y * (a0.x - b0.x);
t_a /= det;

t_b = da.x * (a0.y - b0.y) - da.y * (a0.x - b0.x);
t_b /= det;

if (t_b >= -EPS && t_b <= 1.0 + EPS)
{
// There is intersection
res->x = db.x * t_b + b0.x;
res->y = db.y * t_b + b0.y;
return 1;
}
else
{
return 0;
}
}

void
no_memory ()
{
fprintf (stderr, "No memory\n");
exit (2);
}

inline double
distance_rec (point_t a, point_t b)
{
double t1 = fabs (a.x - b.x);
double t2 = fabs (a.y - b.y);
if (t1 > t2)
return t1;
return t2;
}

poly_t
line_poly_intersect (line_t line, ppoly_t ppoly)
{
poly_t res;
int i;
int size = 2;
int side;
int prev_side;
point_t *pb = &line.b;

res.pts = (point_t *) malloc (size * sizeof (point_t));
res.n = 0;
if (res.pts == NULL)
{
no_memory ();
}
orig = &line.a;

prev_side = cmp (&pb, &ppoly.ppts[ppoly.n - 1]);
for (i = 0; i < ppoly.n; i++)
{
side = cmp (&pb, &ppoly.ppts[i]);
//   fprintf (stderr, "%11g %11g : Side =%d prev = %d\n", ppoly.ppts[i]->x, ppoly.ppts[i]->y, side, prev_side);
if (prev_side != side)
{
if (res.n >= size)
{
int size2 = (res.n * 3 + 1) / 2;
if (res.pts = realloc (res.pts, size2 * sizeof (point_t)))
{
size = size2;
}
else
{
no_memory ();
}
}
if (line_segment_intersect (line.a, line.b, *ppoly.ppts[i - 1], *ppoly.ppts[i], &res.pts[res.n]))
{
if (res.n == 0 || distance_rec (res.pts[res.n], res.pts[res.n - 1]) > EPS)
res.n++;
}
}
prev_side = side;
}

return res;
}

/*
Find intersection points of convex_hull of given points array and a line
Result is represented as poly_t.
inters = convex_intersect(pts, n, line);
inters.n — number of intersection points; could be 0, 1 or 2
inters.pts — array of inters.n intersection points
Call free_poly(inters) at the end;
*/
poly_t
convex_intersect (point_t * pts, int n, line_t line)
{
ppoly_t hull = convex_hull (pts, n);
poly_t inters = line_poly_intersect (line, hull);
free_ppoly (hull);
return inters;
}

int
main (int argc, char *argv[])
{
int n, i;
point_t *pts;
ppoly_t hull;
poly_t inters;
line_t line;
fprintf (stderr,
"Usage: ./convex\n   Read from STDIN points and finds convex hull.\n   Then reads coordinates of two points (point A and point B) of a line\n   and finds intersection points of the line with found convex hull.\n   INPUT FORMAT:\n      n\n      X1 Y1\n      X2 Y2\n      ..\n      Xn Yn\n   \n      A.x A.y B.x B.y\n");

if (scanf ("%d", &n) != 1 || n < 0 || n > 100000000)
{
fprintf (stderr, "Error:  First line should contain number of points\n Now n=%d\n", n);
exit (1);
}
i = 0;
pts = (point_t *) malloc (n * sizeof (point_t));
if (pts == NULL)
no_memory ();

while (i < n && scanf ("%lf%lf", &pts[i].x, &pts[i].y) == 2) i++;

if (i < n)
{
fprintf (stderr, "Warning: wrong number of point in the input.\nDeclared %d points, read only %d points", n, i);
n = i;
}

hull = convex_hull (pts, n);

printf ("Found convex hull consisting of %d points\n", hull.n);
for (i = 0; i < hull.n; i++)
printf ("%11g %11g\n", hull.ppts[i]->x, hull.ppts[i]->y);

printf ("Enter two points (4 numbers: A.x A.y B.x B.y  ) of the line\n");
if( scanf ("%lf%lf%lf%lf", &line.a.x, &line.a.y, &line.b.x, &line.b.y) != 4 )
{
fprintf(stderr, "Cant read four real numbers. Stop\n");
exit(3);
}

inters = line_poly_intersect (line, hull);
free_ppoly (hull);

printf ("Found %d intersections\n", inters.n);
for (i = 0; i < inters.n; i++)
printf ("%11g %11g\n", inters.pts[i].x, inters.pts[i].y);
free_poly (inters);
free (pts);
return 0;
}

Старая версия :
%CODE{"cpp"}%
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define EPS 1E-10
#define N 1000

int n;          //   число точек всего
double x[N][2]; // массив с координатами точек

/*
текущий центр координат (для вычисления вектороного произведения в функции cmp)
*/
double *xt;

double xtstatic[2];

/*
Находит знак векторного произведения векторов, проведеденных из
точки xt в точки a и b
*/
int
cmp (const void *a, const void *b)
{
double *p1, *p2;
double s;
p1 = (double *) a;
p2 = (double *) b;

s = (p1[0] - xt[0]) * (p2[1] - xt[1]) - (p2[0] - xt[0]) * (p1[1] - xt[1]);
if (fabs (s) < EPS)
{
s =
(p1[0] - xt[0]) * (p1[0] - xt[0]) + (p1[1] - xt[1]) * (p1[1] - xt[1]);
s -=
(p2[0] - xt[0]) * (p2[0] - xt[0]) + (p2[1] - xt[1]) * (p2[1] - xt[1]);
if (abs (s) < EPS)
return 0;
if (s < 0)
return 1;
return -1;
};
if (s < 0)
return 1;
return -1;
}

/*
Для точек из глобального массива x находит выпуклую оболочку и кладёт
индексы точек выпуклой оболочки в массив stack
*/
int
convexhull (int *stack)
{
int i_ymin; //  индек точки с минимальной координатой y
int sp;     //  индекс последнего элемента в стеке; в конце
// — число точек в выпуклой оболочке

for (i = 0; i < n; i++){
if (x[i][1] < x[i_ymin][1]) i_ymin = i;
}

xt = xtstatic;
xt[0] = x[i_ymin][0];
xt[1] = x[i_ymin][1];
x[i_ymin][0] = x[0][0];
x[i_ymin][1] = x[0][1];
x[0][0] = xt[0];
x[0][1] = xt[1];

qsort (&x[1], n, 2 * sizeof (double), cmp);

stack[0] = 0;
stack[1] = 1;
stack[2] = 2;

sp = 2;
i = 2;
while (i < n + 1)
{
xt = x[stack[sp - 1]];
if (cmp (&x[stack[sp]][0], &x[(i + 1) % n][0]) < 0)
{
stack[++sp] = (++i) % n;
} else {
sp--;
}
}
return sp;
}

int
main ()
{
int stack[N];
// Scan points from std input
scanf ("%d", &n);
for (i = 0; i < n; i++) scanf ("%lf %lf", &x[i][0], &x[i][1]);

// Place id's of convexhull points into stack
int length = convexhull (stack);
return 0;
}```

-- ArtemVoroztsov - 17 Mar 2004

AlgorithmClasifyForm
Type: Код
Scope: Геометрия
Strategy:
Complexity: Medium