Пересечение линии и параллелепипеда

[avc]

Задачка: Есть бесконечная линия заданная как Point3d и Vector3d и есть габариты некоего объекта заданные как Extents3d. Пересекает ли линия габариты?
Вроде простая задачка, но что-то у меня выходит слишком сложное решение. Может есть готовый быстрый алгоритм на чистой векторной алгебре?
На всякий случай поясню зачем мне это надо: хочу избавится от вызова капризной внешней функции acedNEntSelPEx для поиска объектов под точкой. Но прежде чем искать реальные пересечения линии взгляда со всеми объектами чертежа, хотелось бы по быстрому отбросить объекты, у которых даже GeometricExtents не попадает на линию.
Вот велосипед, который я изобрел:

Код - C# [Выбрать]

1.     /// <summary>
2.     /// Луч в направлении direction пересекает габариты. на самом деле строятся другие габариты в плоскости Plane(picked, direction).
3.     /// При пустом direction вызывате IsInside
4.     /// </summary>
5.     public static bool IsIntersect(this Point3d picked, Vector3d direction, Extents3d ext)
6.     {
7.       if (direction == new Vector3d()) // считаем что все происходит в плоскости XY
8.         return picked.IsInside(ext);
9.       // все вершины прямоугольника заданного Extents3d
10.       Point3d[] extPoints = new Point3d[]{
11.       ext.MinPoint,
12.       new Point3d(ext.MaxPoint.X, ext.MinPoint.Y, ext.MinPoint.Z),
13.       new Point3d(ext.MinPoint.X, ext.MaxPoint.Y, ext.MinPoint.Z),
14.       new Point3d(ext.MinPoint.X, ext.MinPoint.Y, ext.MaxPoint.Z),
15.       new Point3d(ext.MaxPoint.X, ext.MaxPoint.Y, ext.MinPoint.Z),
16.       new Point3d(ext.MaxPoint.X, ext.MinPoint.Y, ext.MaxPoint.Z),
17.       new Point3d(ext.MinPoint.X, ext.MaxPoint.Y, ext.MaxPoint.Z),
18.       ext.MaxPoint
19.       };
20.       // процирование всех вершин на плоскость перпендикулярную direction и поиск габаритов этого множества точек
21.       double minX = double.PositiveInfinity, minY = double.PositiveInfinity, minZ = double.PositiveInfinity;
22.       double maxX = double.NegativeInfinity, maxY = double.NegativeInfinity, maxZ = double.NegativeInfinity;
23.       using (Plane pl = new Plane(picked, direction))
24.         foreach (Point3d p in extPoints)
25.         {
26.           Point3d proj = p.Project(pl, direction);
27.           if (proj.X < minX) minX = proj.X;
28.           else if (proj.X > maxX) maxX = proj.X;
29.           if (proj.Y < minY) minY = proj.Y;
30.           else if (proj.Y > maxY) maxY = proj.Y;
31.           if (proj.Z < minZ) minZ = proj.Z;
32.           else if (proj.Z > maxZ) maxZ = proj.Z;
33.         }
34.       // построение нового Extents3d по габаритам спроецированных точек (заведомо больше чем заданные!)
35.       Extents3d projExt = new Extents3d(new Point3d(minX, minY, minZ), new Point3d(maxX, maxY, maxZ));
36.       return picked.IsInside(projExt);
37.     }
38.  
39.     /// <summary>
40.     /// Gets a value indicating whether the specified point is inside the extents.
41.     /// </summary>
42.     /// <param name="pt">The instance to which the method applies.</param>
43.     /// <param name="extents">The extents 3d supposed to contain the point.</param>
44.     /// <returns>true if the point is inside the extents; otherwise, false.</returns>
45.     public static bool IsInside(this Point3d pt, Extents3d extents, double tolerance = 0.0000001)
46.     {
47.       return
48.         (pt.X > extents.MinPoint.X - tolerance) &&
49.         (pt.Y > extents.MinPoint.Y - tolerance) &&
50.         (pt.Z > extents.MinPoint.Z - tolerance) &&
51.         (pt.X < extents.MaxPoint.X + tolerance) &&
52.         (pt.Y < extents.MaxPoint.Y + tolerance) &&
53.         (pt.Z < extents.MaxPoint.Z + tolerance);
54.     }

Код рабочий, но наверно не быстрей, чем сразу вызывать brep.GetLineContainment. Особенно напрягает необходимость создавать и диспозить Plane - любой DisposableWrapper у меня вызывает подозрения

[Александр Ривилис]

Extents3d - это параллелепипед, а не прямоугольник.

[avc]

Extents3d - это параллелепипед, а не прямоугольник.

Долго искал, где ж я так опечатался Исправил

[avc]

Кстати насчет прямоугольников. Вариант поиска пересечения линии с шестью ортогональными прямоугольниками-гранями этого параллелепипеда я тоже рассматривал, но даже не смог подступиться. Гуглятся школьные задачки с описанием линии формулами и т.п.

[Александр Ривилис]

Думаю, что твой вариант самый быстрый, но не самый точный.
Кстати, для определения нахождения точки внутри параллелепипеда есть класс BoundBlock3d с методом Contains

[avc]

Спасибо за весомое мнение 
Буду тестировать, что производительней: сразу искать пересечения, проверять Extents3d, или же (по вашей наводке) попробовать метод BoundBlock3d.IsDisjoint c линией описанной как BoundBlock3d с двумя нулевыми размерами...

[Дима_]

Не проверял - но по моему так должно работать (возможно, на больших расстояниях и малых отклонениях не будет хватать точности - тогда вектор надо приводить не к 1, а к максимально используемому значению - но это вопрос "притирки-наладки" - идея в коде):

Код - C# [Выбрать]

1.         private static Vector3d UnitVector(Vector3d v)
2.         {
3.             var arr = v.ToArray();
4.             var mn = arr[v.LargestElement];
5.             return new Vector3d(arr.Select(itm => itm / mn).ToArray());
6.         }
7.  
8.         public static bool IsIntersect(Point3d pt, Vector3d vct, Extents3d ext)
9.         {
10.             Point3d[] extPoints = new Point3d[]{
11.                                                   ext.MinPoint,
12.                                                   new Point3d(ext.MaxPoint.X, ext.MinPoint.Y, ext.MinPoint.Z),
13.                                                   new Point3d(ext.MinPoint.X, ext.MaxPoint.Y, ext.MinPoint.Z),
14.                                                   new Point3d(ext.MinPoint.X, ext.MinPoint.Y, ext.MaxPoint.Z),
15.                                                   new Point3d(ext.MaxPoint.X, ext.MaxPoint.Y, ext.MinPoint.Z),
16.                                                   new Point3d(ext.MaxPoint.X, ext.MinPoint.Y, ext.MaxPoint.Z),
17.                                                   new Point3d(ext.MinPoint.X, ext.MaxPoint.Y, ext.MaxPoint.Z),
18.                                                   ext.MaxPoint
19.                                                 };
20.             var v = UnitVector(vct);
21.             var vcs = extPoints.Select(ptt => UnitVector(pt.GetVectorTo(ptt))).ToList();
22.             return vcs.Any(vc => v.X >= vc.X && v.Y >= vc.Y && v.Z >= vc.Z)&&
23.                    vcs.Any(vc => v.X <= vc.X && v.Y <= vc.Y && v.Z <= vc.Z);
24.         }

[Александр Ривилис]

Дима_,
А зачем нужен метод UnitVector если есть стандартный Vector3d.GetNormal() ?

[Дима_]

Дима_,
А зачем нужен метод UnitVector если есть стандартный Vector3d.GetNormal() ?

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

[Дима_]

Стоп - а GetNormal - это кажется не то? Мне единичный вектор нужен - сейчас посмотрю...

[Александр Ривилис]

Стоп - а GetNormal - это кажется не то? Мне единичный вектор нужен - сейчас посмотрю...

Посмотри-посмотри.

[avc]

Бррр. Хоть убей не усваиваю суть кода без типизации и на лямбдах. Но непременно попробую as is.  Спасибо.
GetNormal это как ни странно не нормаль какая-нибудь, а как раз единичный вектор, хотя по названию никак не догадаться, да.
Пока успел потестировать скорость с проверкой Extents и без нее. На моем велосипеде проверка не ускоряет, а затормаживает работу от полутора (на блоках) до 12 раз (на солидах). Похоже IntersectWith и GetLineContainment и так очень быстрые. Тысячи объектов за доли секунды. Это под отладчиком. Жаль только IntersectWith реализован только у половины Entity

[Дима_]

Посмотри-посмотри.

Как ни странно - но нужен "тот" единичный вектор - который возвращает UnitVector. Во вложении пример который "не ловит" GetNormal.

[Александр Ривилис]

Как ни странно - но нужен "тот" единичный вектор - который возвращает UnitVector.

Еще раз посмотрел его код. Он возвращает не единичный вектор в смысле аналитической геометрии.

[Дима_]

Хоть убей не усваиваю суть кода без типизации и на лямбдах

Надо привыкать мыслить абстрактно (не обращая внимание на синтаксис - хотя он и иногда раздражает). Суть простая - смотрим есть ли среди 8 единичных векторов (от исходной точки - до 8 угловых точек), такие, у которых все "координаты" больше и все меньше (то есть - чтоб был и тот и другой), чем у нашего вектора.