Координаты середины полилинии

[Prot_XT]

Добрый день!

Помогите решить задачу: имеется полилиния, допустим с тремя точками (координаты известны), как найти координаты середины полилинии?
Для линии с двумя точками это делаю просто - из векторов по алгебре - координаты середины вектора = (x1+x2)/2 и (y1+y2)/2
А как это сделать когда много точек? То есть линия состоит из множества векторов?

[Boris]

для случая, когда полилиния состоит множества векторов, понятие средней точки не может быть определено однозначно. За исключением тех случаев, когда реально это одна протяженная линия, но нарисованная участками. В случае ветвления - средняя точка зависит от предпочтений смотрящего. На форуме уже велось обсуждение того как разные ГИС вычисляют центральную точку полигона, для случая ветвящейся линии, решения точно такие же. Для случая же единой гладкой полилинии без ветвлений, решение такое же как и для отрезка - это точка отстоящая от начала - конца на расстояние равное половине длины линии. Т.е. сумма длин всех отрезков, деленная пополам. В зависимости от того как производится вычисление, с учетом кривизны Земли или нет, и в какой проекции, эта точка может меняться в некоторых пределах.

[Prot_XT]

Сумма длин всех отрезков пополам - это мне понятно. Причем это мое понимание длится до тех пор пока линия прямая. А если она на какой то точке уходит в сторону? Тогда, мое понимание пропадает. Э... проекция в моем случае ни причем. Нужно общее понимание. То есть в любой проекции простой вектор a1-a2 с двумя точками вот с таким вычислением: (x1+x2)/2 и (y1+y2)/2 отображает мне точку точно посередине линии. А вот имеется ли такая подобная формула для линии с тремя точками, причем допустим первый отрезок короче второго и посередине во второй точке линия имеет излом на некий угол?

[Boris]

вам для начала надо определиться с тем, что вы понимаете под серединой полилинии, или выбрать себе из имеющихся определений.
То о чем вы говорите, есть 1/2 суммы векторов, задающих сегменты полилинии. То о чем говорю я есть:
- точка лежащая на полилинии.
- точка расстояние до которой, с обеих концов полилинии одинаково.
Для таких операций существует специальный раздел ГИС, называемый "системы линейных координат", т.е. координат, отсчитываемых вдоль некоторой линии (1 мерного объекта), расположенного в N-мерном пространстве. В вашем случае 2-х или 3-х, если вы учитываете и высоту.Это полный аналог ситуации с автомобильной дорогой - где находится середина пути из Москвы в СПб, зависит от того как пролегает путь, а не от того, где находится середина прямой их соединяющей.
На сайте ESRI есть книжка посвященная этому, есть она и на сайте Дата+. Может какие добрые люди на форуме знают еще материалы, посвященные «Системы линейных координат», в но уже открытом доступе.
Если же все изложенное не есть ваша задача, то тогда вам надо ее как-то точнее сформулировать.

[Denis Rykov]

Видимо, автор темы хочет вычислить координаты точки аналогично данной функции PostGIS? Если так, то может стоит глянуть исходники как оно там реализовано?

[Prot_XT]

Да, возможно я не точно сформулировал свой вопрос. Мне нужно найти координаты точки находящейся посередине некой ломаной линии или полилинии. Пространство 2-х мерное, высоту учитывать не нужно. И действительно это аналог с серединой пути от Москвы до СПб, в зависимости от того как он пролегает.
Буду искать литературу на тему систем линейных координат. После нескольких попыток в интернете понял, что найти не так уж и легко будет...
А я почему то думал, что это буквально простая формула, которые все ГИС гуру знают наизусть - а я один не знаю...

[Olga_@@@]

Автокадовская команда _divide не подойдет?
Поставит точку посередине полилинии с учетом всех загогулин, останется только получить координаты точки.
Если нужен алгоритм, спросите здесь или на dwg.ru

[Prot_XT]

Нет, мне нужен алгоритм...

[Olga_@@@]

ссылка в сообщении выше, там ответят

[Prot_XT]

Видимо, автор темы хочет вычислить координаты точки аналогично данной функции PostGIS? Если так, то может стоит глянуть исходники как оно там реализовано?

Да, это то что мне нужно...
Как бы узнать алгоритм этого "чудесного" вычисления...

[Denis Rykov]

Иcходный код PostGIS доступен в svn.

Спойлер

/***********************************************************************
* --strk@keybit.net;
***********************************************************************/

/***********************************************************************
* Interpolate a point along a line, useful for Geocoding applications
* SELECT line_interpolate_point( 'LINESTRING( 0 0, 2 2'::geometry, .5 )
* returns POINT( 1 1 ).
* Works in 2d space only.
*
* Initially written by: jsunday@rochgrp.com
* Ported to LWGEOM by: strk@refractions.net
***********************************************************************/

Datum LWGEOM_line_interpolate_point(PG_FUNCTION_ARGS);

PG_FUNCTION_INFO_V1(LWGEOM_line_interpolate_point);
Datum LWGEOM_line_interpolate_point(PG_FUNCTION_ARGS)
{
PG_LWGEOM *geom = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
double distance = PG_GETARG_FLOAT8(1);
LWLINE *line;
LWPOINT *point;
POINTARRAY *ipa, *opa;
POINT4D pt;
int nsegs, i;
double length, slength, tlength;

if ( distance < 0 || distance > 1 )
{
elog(ERROR,"line_interpolate_point: 2nd arg isnt within [0,1]");
PG_RETURN_NULL();
}

if ( pglwgeom_get_type(geom) != LINETYPE )
{
elog(ERROR,"line_interpolate_point: 1st arg isnt a line");
PG_RETURN_NULL();
}

line = lwgeom_as_lwline(pglwgeom_deserialize(geom));
ipa = line->points;

/* If distance is one of the two extremes, return the point on that
* end rather than doing any expensive computations
*/
if ( distance == 0.0 || distance == 1.0 )
{
if ( distance == 0.0 )
getPoint4d_p(ipa, 0, &pt);
else
getPoint4d_p(ipa, ipa->npoints-1, &pt);

opa = ptarray_construct_reference_data(FLAGS_GET_Z(line->flags), FLAGS_GET_M(line->flags), 1, (uchar*)&pt);

point = lwpoint_construct(line->srid, NULL, opa);
PG_RETURN_POINTER(pglwgeom_serialize(lwpoint_as_lwgeom(point)));
}

/* Interpolate a point on the line */
nsegs = ipa->npoints - 1;
length = ptarray_length_2d(ipa);
tlength = 0;
for ( i = 0; i < nsegs; i++ )
{
POINT4D p1, p2;
POINT4D *p1ptr=&p1, *p2ptr=&p2; /* don't break
* strict-aliasing rules
*/

getPoint4d_p(ipa, i, &p1);
getPoint4d_p(ipa, i+1, &p2);

/* Find the relative length of this segment */
slength = distance2d_pt_pt((POINT2D*)p1ptr, (POINT2D*)p2ptr)/length;

/* If our target distance is before the total length we've seen
* so far. create a new point some distance down the current
* segment.
*/
if ( distance < tlength + slength )
{
double dseg = (distance - tlength) / slength;
interpolate_point4d(&p1, &p2, &pt, dseg);
opa = ptarray_construct_reference_data(FLAGS_GET_Z(line->flags), FLAGS_GET_M(line->flags), 1, (uchar*)&pt);
point = lwpoint_construct(line->srid, NULL, opa);
PG_RETURN_POINTER(pglwgeom_serialize(lwpoint_as_lwgeom(point)));
}
tlength += slength;
}

/* Return the last point on the line. This shouldn't happen, but
* could if there's some floating point rounding errors. */
getPoint4d_p(ipa, ipa->npoints-1, &pt);
opa = ptarray_construct_reference_data(FLAGS_GET_Z(line->flags), FLAGS_GET_M(line->flags), 1, (uchar*)&pt);
point = lwpoint_construct(line->srid, NULL, opa);
PG_RETURN_POINTER(pglwgeom_serialize(lwpoint_as_lwgeom(point)));
}
/***********************************************************************
* --jsunday@rochgrp.com;
***********************************************************************/

Думаю, знактоки С(си) смогут прокомментировать данный код.

[Boris]

Видимо задачи где-то выдают по 2-3 одинаковых. Только неделю назад вопрошавший искал на форуме как в MapInfo поставить точку на удалении A от начала линии.
В чем собственно состоит сложность алгоритма? Все формулы, за исключение формулы расчета расстояний на шаре, представляют собой курс средней школы:

• Надо посчитать число сегментов в линии, что бы было = 1, иначе решение не однозначно, или линию требуется "сшить".
• Надо посчитать N-1 расстояний между точками, где N - количество вершин задающих линию. Что бы узнать полную длину линии. Во многих ГИС эту величину возвращает некоторая функция от гео-объекта.
• Результат измерений поделить пополам и получить некоторое число равное длине полу-линии.
• Выбрать начало отсчета. Но для случая "пополам" - это все равно.
• Найти пару вершин, для которых сумма расстояний отрезков от начала отсчета будет для вершины M меньше искомой величины, а для вершины M+1 - больше. Если найдется вершина для которой равно, то и считать дальше не надо.
• На этом последнем отрезке [M; M+1] отмерить сколько осталось, если вычесть из искомой величины расстояние до вершины M в направлении вершины M+1.

Единственная сложность - как считать расстояния, спроецированными или на глобусе.
Проблема, если координаты вершин заданы в градусной СК.

[Prot_XT]

Спасибо, вроде бы понял.
Теперь нужно все это отработать на практике...