Реализация алгоритма Hoey Shamos с C#

Хорошо, теперь я получаю правильную информацию из моего текущего алгоритма! Однако, с проверкой 700000 полигонов, это слишком медленно! Предыдущая проблема исправлена ​​(My Line2D intersectsWith метод был неправильным)

Теперь нужно определить мое узкое место! Предполагается, что этот алгоритм будет O(nlog-n), поэтому он должен быть намного быстрее. Мой метод intersectsWith выглядит так, как будто он не может работать быстрее, однако я опубликую его код на случай, если я ошибаюсь

РЕДАКТИРОВАТЬ: Добавлен интерфейс IComparable

Мой метод для чтения пересечения отрезков. Некоторый код был опущен для удобства чтения.

    public class Line2D : IComparable
    {


    public Line2D(XYPoints p1, XYPoints p2)
    {

    }

    public bool intersectsLine(Line2D comparedLine)
    {

        if ((X2 == comparedLine.X1) && (Y2 == comparedLine.Y1)) return false;
        if ((X1 == comparedLine.X2) && (Y1 == comparedLine.Y2)) return false;

        if (X2 == comparedLine.X1 && Y2 == comparedLine.Y1)
        {
            return false;
        }

        if (X1 == comparedLine.X2 && Y1 == comparedLine.Y2)
        {
            return false;
        }
        double firstLineSlopeX, firstLineSlopeY, secondLineSlopeX, secondLineSlopeY;

        firstLineSlopeX = X2 - X1;
        firstLineSlopeY = Y2 - Y1;

        secondLineSlopeX = comparedLine.getX2() - comparedLine.getX1();
        secondLineSlopeY = comparedLine.getY2() - comparedLine.getY1();

        double s, t;
        s = (-firstLineSlopeY * (X1 - comparedLine.getX1()) + firstLineSlopeX * (getY1() - comparedLine.getY1())) / (-secondLineSlopeX * firstLineSlopeY + firstLineSlopeX * secondLineSlopeY);
        t = (secondLineSlopeX * (getY1() - comparedLine.getY1()) - secondLineSlopeY * (getX1() - comparedLine.getX1())) / (-secondLineSlopeX * firstLineSlopeY + firstLineSlopeX * secondLineSlopeY);

        if (s >= 0 && s <= 1 && t >= 0 && t <= 1)
        {
            return true;
        }

        return false; // No collision 
    }

    int IComparable.CompareTo(object obj)
    {

        //return Y1.GetHashCode();
        Line2D o1 = this;
        Line2D o2 = (Line2D)obj;
        if (o1.getY1() < o2.getY1())
        {
            return -1;
        }
        else if (o1.getY1() > o2.getY2())
        {
            return 1;
        }
        else
        {
            if (o1.getY2() < o2.getY2())
            {
                return -1;
            }
            else if (o1.getY2() > o2.getY2())
            {
                return 1;
            }
            else
            {
                return 0;
            }
        } 
    }


}

Большая часть моей реализации алгоритма, я понимаю, что List не самый быстрый для алгоритма, однако мне нужно индексирование!:

//Create a new list, sort by Y values.

 List<AlgEvent> SortedList = events.OrderBy(o => o.getY()).ToList();                
 List<Line2D> sweepline = new List<Line2D>();

 for (var g = 0; g < SortedList.Count; g++)
 {
 if (SortedList[g].isStart)
 {
    Line2D nl = SortedList[g].line;
    Line2D above;
    /* Start generating above */
    try
    {
        //grab index in sweepline
        int index = sweepline.IndexOf(nl);
        //add 1 to get above line
        if (index == -1)
        {
            above = null;
        }
        else
        {
            above = sweepline[index + 1];
        }


    }
    catch (ArgumentOutOfRangeException)
    {
        above = null;
    }
    /* End generating above */
    if (above != null)
    {
        if (above.intersectsLine(nl))
        {
            return true;
        }
    }
    Line2D below;
    /* Start generating below */
    try
    {
        //grab index in sweepline
        int index = sweepline.IndexOf(nl);
        //add 1 to get above line
        below = sweepline[index - 1];

    }
    catch (ArgumentOutOfRangeException)
    {

        below = null;

    }
    /* End generating below */
    if (below != null)
    {
        if (below.intersectsLine(nl))
        {
            return true;
        }
    }
    sweepline.Add(nl);
    sweepline = sweepline.OrderBy(o => o.getY1()).ToList();

}
else
{
    Line2D nl = SortedList[g].line;
    Line2D above;
    Line2D below;
    /* Start generating above */
    try
    {
        //grab index in sweepline
        int index = sweepline.IndexOf(nl);
        Console.Out.WriteLine("index:" + index);
        //add 1 to get above line
        above = sweepline[index + 1];

    }
    catch (ArgumentOutOfRangeException)
    {

        above = null;

    }
    /* End generating above */
    /* Start generating below */
    try
    {
        //grab index in sweepline
        int index = sweepline.IndexOf(nl);
        //add 1 to get above line
        below = sweepline[index - 1];

    }
    catch (ArgumentOutOfRangeException)
    {

        below = null;

    }
    /* End generating below */
    sweepline = sweepline.OrderBy(o => o.getY1()).ToList();
    sweepline.Remove(nl);
    if (above != null && below != null)
    {
        if (above.intersectsLine(below))
        {
            return true;
        }
    }
}
Console.WriteLine("");
  }



   } // end numofparts for-loop

   return false;

============================================

ОБНОВЛЕНИЕ: 12 сентября: Реализовал TreeSet из C5, реализовал IComparable для моих классов и еще больше замедлил его? Я все еще индексирую это, если это имеет значение?

http://www.itu.dk/research/c5/

Код с использованием TreeSet:

TreeSet<Line2D> sweepline = new TreeSet<Line2D>();
for (var g = 0; g < SortedList.Count; g++)
{
if (SortedList[g].isStart)
{
    Line2D nl = SortedList[g].line;
    Line2D above;
    /* Start generating above */
    try
    {
        //grab index in sweepline
        int index = sweepline.IndexOf(nl);
        //add 1 to get above line
        above = sweepline[index + 1];

    }
    catch (IndexOutOfRangeException)
    {

        above = null;

    }
    /* End generating above */
    if (above != null)
    {
        if (above.intersectsLine(nl))
        {
            return false;
        }
    }
    Line2D below;
    /* Start generating below */
    try
    {
        //grab index in sweepline
        int index = sweepline.IndexOf(nl);
        //add 1 to get above line
        below = sweepline[index - 1];

    }
    catch (IndexOutOfRangeException)
    {

        below = null;

    }
    /* End generating below */
    if (below != null)
    {
        if (below.intersectsLine(nl))
        {
            return false;
        }
    }
    sweepline.Add(nl);
    //sweepline = sweepline.OrderBy(o => o.getY1()).ToList();

}
else
{
    Line2D nl = SortedList[g].line;
    Line2D above;
    Line2D below;
    /* Start generating above */
    try
    {
        //grab index in sweepline
        int index = sweepline.IndexOf(nl);
        //Console.Out.WriteLine("index:" + index);
        //add 1 to get above line
        above = sweepline[index + 1];

    }
    catch (IndexOutOfRangeException)
    {

        above = null;

    }
    /* End generating above */
    /* Start generating below */
    try
    {
        //grab index in sweepline
        int index = sweepline.IndexOf(nl);
        //add 1 to get above line
        below = sweepline[index - 1];

    }
    catch (IndexOutOfRangeException)
    {

        below = null;

    }
    /* End generating below */
    //sweepline = sweepline.OrderBy(o => o.getY1()).ToList();
    sweepline.Remove(nl);
    if (above != null && below != null)
    {
        if (above.intersectsLine(below))
        {
            return false;
        }
    }
}
//Console.WriteLine("");

}

Источник

2 ответа

Решение

Во-первых, что касается пересечения линии: вам не нужна фактическая точка пересечения, только чтобы знать, пересекаются ли они. См. http://www.geeksforgeeks.org/check-if-two-given-line-segments-intersect/ для алгоритма, который делает именно это.

О List реализация:

В вашей реализации используется Listс, вы звоните indexOf на подметать, чтобы найти nl, Это ищет линию стрельбы от начала до конца. Увидеть List(T).IndexOf, Если бы вы использовали BinarySearch метод, который должен значительно ускорить поиск.

В документации List есть параграф, который называется соображениями производительности. Они призывают вас использовать тип значения, который реализует IEquatable<T> а также IComparable<T>, Так что ваши Line2D должно быть struct и реализовать эти интерфейсы.

Если вы последуете этому совету, извлечение конечной точки из линии развертки должно быть O(log n), что достаточно для вашей цели, и память должна использоваться более эффективно.

Вставка и удаление - O(n) для Списков, потому что основной массив должен быть перемещен в памяти. SortedSet имеет более быструю вставку и удаление, но я не совсем понимаю, как найти соседей элемента в O (log n) там. Кто-нибудь? (См. Также Почему SortedSet .GetViewBetween не O (log N)?)

В любом случае, С5 TreeSet должен решить это.

Я посмотрел производительность IndexOf и [i] в ​​руководстве пользователя, и они оба перечислены как O(log n). Так что это не должно быть проблемой. Вероятно, все еще несколько быстрее, но не более чем с фиксированным коэффициентом, вызывать специальные методы для поиска соседей на линии поворота, то есть, наследники и предшественники, которые также являются O(log n).

Так

[...]
try 
{
    Line2D above = sweepline.Successor(nl);
    if (above.intersectsLine(nl))
    {
        return false;
    }
}
catch (NoSuchItemException ignore) { }
[...]

Мне не нравится, что у них нет метода, который не генерирует исключение, так как создание исключений очень дорого. В общем, ваша стрела будет довольно полной, поэтому я думаю, что неудача в ее поиске будет редкой и вызывающей Successor это самый эффективный способ. В качестве альтернативы, вы можете продолжать звонить IndexOf как и сейчас, но проверьте, равно ли оно Count минус один перед извлечением [index + 1] и предотвратите генерацию исключения:

[...]
int index = sweepline.IndexOf(nl);
if( index < sweepline.Count-1 )
{
    Line2D above = sweepline[index + 1];
    if (above.intersectsLine(nl))
    {
        return false;
    }
}
[...]

Вторая глава руководства описывает равенство и сравнение для коллекций C5. Здесь тоже, по крайней мере, вы должны реализовать IEquatable<T> а также IComparable<T>!

Еще одна мысль: вы сообщаете алгоритму подачи 700000 строк. Не могли бы вы начать с синхронизации, например, 1000, 2500, 5000, 10000 строк и посмотреть, как алгоритм масштабируется для случаев, когда они не пересекаются?

О том, как сравнить линии на линии разметки:

Вам нужно найти какой-то естественный порядок для Line2D на Sweepline TreeSet, так как метод CompareTo просит вас сравнить один Line2D с другим. Один из Line2D уже находится в Sweepline TreeSet, другой только что обнаружен и добавляется.

Я думаю, ваша линия подметания проходит снизу вверх:

List<AlgEvent> SortedList = events.OrderBy(o => o.getY()).ToList();

Sweepline проходит через отрезки

Итак, скажем, сегмент S1 был добавлен в TreeSet при событии 1, и мы хотим сравнить его с S2, который добавляется при событии 2, прямо сейчас.

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

В любом случае... так что проще всего будет сравнить нижние конечные точки обоих отрезков. Слева меньше, справа больше. Однако нам нужно взглянуть на порядок в позиции линии поворота, и с тех пор они могли изменить положение, как на рисунке.

Таким образом, нам нужно сравнить нижнюю конечную точку S2 с красной точкой на S1 и посмотреть, находится ли она слева или справа от этой точки. Он расположен слева, поэтому S2 считается меньше, чем S1.

Обычно это проще, чем это: если весь S1 находится слева от нижней конечной точки S2, S1 меньше, чем S2. Если весь S1 находится справа от нижней конечной точки S2, S1 больше, чем S2.

Я думаю, что вы ищете более безопасную версию интерфейса:

public class Line2D : IComparable<Line2D>

Предполагая два свойства BottomY самое низкое из двух значений Y, и BottomX, значение X самой низкой конечной точки, несколько проверенная попытка:

int IComparable<Line2D>.CompareTo(Line2D other)
{
    if( BottomY < other.BottomY )
    {
        return -other.CompareTo(this);
    }

    // we're the segment being added to the sweepline
    if( BottomX >= other.X1 && BottomX >= other.X2 )
    {
        return 1;
    }
    if( BottomX <= other.X1 && BottomX <= other.X2 )
    {
        return -1;
    }

    if( other.Y2 == other.Y1 )
    {
        // Scary edge case: horizontal line that we intersect with. Return 0?
        return 0;
    }

    // calculate the X coordinate of the intersection of the other segment
    // with the sweepline
    // double redX = other.X1 + 
    //    (BottomY - other.Y1) * (other.X2 - other.X1) / (other.Y2 - other.Y1);
    //
    // return BottomX.CompareTo(redX);

    // But more efficient, and more along the lines of the orientation comparison:
    return Comparer<Double>.Default.Compare(
        (BottomX - other.X1) * (other.Y2 - other.Y1),
        (BottomY - other.Y1) * (other.X2 - other.X1) );

}
Источник

[оригинальный ответ]

Я не пользователь C#, но это должно немного ускорить процесс.

  • меньше кучи мусора
  • не вычисляйте одно и то же дважды
  • избегайте всех дополнительных вызовов, если можете (удалите функции)

код:

public bool intersectsLine(const Line2D &comparedLine)
    {
    if ((X2==comparedLine.X1)&&(Y2==comparedLine.Y1)) return false;
    if ((X1==comparedLine.X2)&&(Y1==comparedLine.Y2)) return false;
    double dx1,dy1,dx2,dy2;
    dx1 = X2 - X1;
    dy1 = Y2 - Y1;
    dx2 = comparedLine.X2 - comparedLine.X1;
    dy2 = comparedLine.Y2 - comparedLine.Y1;
    double s,t,ax,ay,b;
    ax=X1-comparedLine.X1;
    ay=Y1-comparedLine.Y1;
    b=1.0/(-(dx2*dy1)+(dx1*dy2));
    s = (-(dy1*ax)+(dx1*ay))*b;
    t = ( (dx2*ay)-(dy2*ax))*b;
    if ((s>=0)&&(s<=1)&&(t>=0)&&(t<=1)) return true;
    return false; // No collision
    }

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

[edit1]

После некоторых исследований случайных данных я пришел к выводу:

  1. если по всей области слишком много строк, тогда никакие оптимизации недопустимы
  2. если больше маленьких строк, чем ускорение для любых оптимизаций
  3. грубая сила T((N*N-N)/2) который до сих пор O(N*N)приблизительно 700 часов для обработки 700 тыс. строк (неприменимо)

  4. Оптимизированный перебор с делением площади T((((N/M)^2)-N)/2) - оптимизации ~O((N/M)^2) где

    • N максимальное количество линий области

    • M является количеством делений площади на любую ось. Идея состоит в том, чтобы проверять только линии, пересекающие некоторый регион (разделить область набора данных на M*M квадраты / прямоугольники). Для 700К линий лучше всего подходит 16x16 области. Измеренные времена:

      0.540 с на строки 32 КБ 1.950 с на строки 64 КБ 7.000 с на строки 128 К 27.514 с на строки 256 КБ

    расчетное время выполнения составляет 3,7 мин на линии 700 КБ (для линий длиной не более 10% всей площади). Я думаю, что это лучше, чем твои 19 минут.

  5. другое ускорение возможно с использованием нескольких процессоров / ядер

    Алгоритм полностью параллельный и для 4 процессоров / ядер 3.7min/4 -> 56s или вы можете портировать его на GPU ...

Мой оптимизированный алгоритм перебора с делением площади O((((N/M)^2)-N)/2) - оптимизации

  1. получить размер используемой площади (xmin,xmax,ymin,ymax)O(N)
  2. выберите подразделение Mлучшее, что я пытаюсь для моих случайных наборов данных 32K-256K линии было M=16

  3. цикл по всей области подразделения (равномерно разделенная область набора данных)

    создать список линий, пересекающих фактическую область подразделения, и проверить пересечение для всех линий в этом списке. Если вы не хотите дублировать пересечения, отбросьте все пересечения за пределами текущей области.

мой код (я использую BDS2006 C++ и мои собственные списки, поэтому вам нужно перенести его, чтобы он был совместим с вашим кодом)

void Twin_GLView2D::main_intersect(int M=16)
{
int ia,ib,i,j,N;
double zero=1e-6;
glview2D::_lin *l;
glview2D::_pnt p;

struct _line
    {
    double bx0,by0,bx1,by1;     // bounding rectangle
    double x0,y0,dx,dy;         // precomputed params
    } *lin,*a,*b;

struct _siz
    {
    double bx0,bx1,by0,by1;     // zone bounding rectangle
    } sz,bz;
List<_line*> zone;

// load and precompute lines
N=view.lin.num;
lin=new _line[N];
if (lin==NULL) return;
for (a=lin,l=view.lin.dat,ia=0;ia<N;ia++,a++,l++)
    {
    // line ...
    if (l->p0.p[0]<=l->p1.p[0]) { a->bx0=l->p0.p[0]; a->bx1=l->p1.p[0]; }
    else                        { a->bx0=l->p1.p[0]; a->bx1=l->p0.p[0]; }
    if (l->p0.p[1]<=l->p1.p[1]) { a->by0=l->p0.p[1]; a->by1=l->p1.p[1]; }
    else                        { a->by0=l->p1.p[1]; a->by1=l->p0.p[1]; }
    a->x0=l->p0.p[0]; a->dx=l->p1.p[0]-l->p0.p[0];
    a->y0=l->p0.p[1]; a->dy=l->p1.p[1]-l->p0.p[1];
    // global image size for zone subdivision
    if (!ia)
        {
        sz.bx0=l->p0.p[0];
        sz.by0=l->p0.p[1];
        sz.bx1=sz.bx0;
        sz.by1=sz.by0;
        }
    if (sz.bx0>l->p0.p[0]) sz.bx0=l->p0.p[0];
    if (sz.bx1<l->p0.p[0]) sz.bx1=l->p0.p[0];
    if (sz.by0>l->p0.p[1]) sz.by0=l->p0.p[1];
    if (sz.by1<l->p0.p[1]) sz.by1=l->p0.p[1];
    if (sz.bx0>l->p1.p[0]) sz.bx0=l->p1.p[0];
    if (sz.bx1<l->p1.p[0]) sz.bx1=l->p1.p[0];
    if (sz.by0>l->p1.p[1]) sz.by0=l->p1.p[1];
    if (sz.by1<l->p1.p[1]) sz.by1=l->p1.p[1];
    }
// process lines by zonal subdivision
zone.allocate(N);
view.pnt.num=0; view.pnt.allocate(view.lin.num);
sz.bx1-=sz.bx0; sz.bx1/=double(M);
sz.by1-=sz.by0; sz.by1/=double(M);
for (bz.by0=sz.by0,bz.by1=sz.by0+sz.by1,i=0;i<M;i++,bz.by0+=sz.by1,bz.by1+=sz.by1)
for (bz.bx0=sz.bx0,bz.bx1=sz.bx0+sz.bx1,j=0;j<M;j++,bz.bx0+=sz.bx1,bz.bx1+=sz.bx1)
    {
    // create list of lines for actual zone only
    zone.num=0;         // clear zone list
    for (a=lin,ia=   0;ia<N;ia++,a++)
     if ((a->bx0<=bz.bx1)&&(a->bx1>=bz.bx0))
      if ((a->by0<=bz.by1)&&(a->by1>=bz.by0))
       zone.add(a); // add line to zone list
    // check for intersection within zone only
    // O((((N/M)^2)-N)/2) - optimizations
    for (ia=   0,a=zone.dat[ia];ia<zone.num;ia++,a=zone.dat[ia])
    for (ib=ia+1,b=zone.dat[ib];ib<zone.num;ib++,b=zone.dat[ib])
        {
        // discart lines with non intersecting bound rectangles
        if (a->bx1<b->bx0) continue;
        if (a->bx0>b->bx1) continue;
        if (a->by1<b->by0) continue;
        if (a->by0>b->by1) continue;
        // 2D lines a,b intersect ?
        double x0,y0,x1,y1,t0,t1;
        // compute intersection
        t1=divide(a->dx*(a->y0-b->y0)+a->dy*(b->x0-a->x0),(a->dx*b->dy)-(b->dx*a->dy));
        x1=b->x0+(b->dx*t1);
        y1=b->y0+(b->dy*t1);
        if (fabs(a->dx)>=fabs(a->dy)) t0=divide(b->x0-a->x0+(b->dx*t1),a->dx);
        else                          t0=divide(b->y0-a->y0+(b->dy*t1),a->dy);
        x0=a->x0+(a->dx*t0);
        y0=a->y0+(a->dy*t0);
        // check if intersection exists
        if (fabs(x1-x0)>zero) continue;
        if (fabs(y1-y0)>zero) continue;
        if ((t0<0.0)||(t0>1.0)) continue;
        if ((t1<0.0)||(t1>1.0)) continue;
        // if yes add point
        p.p[0]=x0;
        p.p[1]=y0;
        p.p[2]=0.0;
        // do not add points out of zone (allmost all duplicit points removal)
        if (x0<bz.bx0) continue;
        if (x0>bz.bx1) continue;
        if (y0<bz.by0) continue;
        if (y0>bz.by1) continue;
        view.pnt.add(p);
        }
    }
view.redraw=true;
delete lin;
}

Заметки:

  1. List<T> x; такой же как T x[] с неограниченным размером
    • x.num; фактический размер x[] в т не байты!!! index = <0,x.num-1>
    • x.add(q); добавляет q к списку в конце
    • x.num=0; очищает список
    • x.allocate(N); выделить место для N элементы в списке, чтобы избежать перемещения
  2. вход List<> является view.lin... содержит баллы p0,p1 у каждого есть double p[2]... x,y
  3. выход List<> является view.pnt... содержит double p[2]... x,y

[Edit2]

Кроме того, я обнаружил, что лучшая производительность вышеуказанного алгоритма - это когда M=12+(N>>15)

  • M подсчет площадей на ось
  • N количество строк для проверки
Источник
Другие вопросы по тегам