Пересечение трехмерной сетки точек движущейся сферой

@MeXaL · Регистрация: 10.09.2023

Студворк — интернет-сервис помощи студентам

Дана трехмерная сетка из точек, выровненных по осям X, Y, Z с равномерным шагом ∆s начиная с заданной опорной точки O. Любая точка P в этой сетке может быть определена как: P = O + (ix, iy, iz)∆s, де ix, iy и iz - целочисленные индексы, принадлежащие открытым диапазонам [0, nx), [0, ny) и [0, nz) соответственно. Сетка пересекается движущейся сферой радиуса R. Движение центра сферы определяется 3d-кривой f(t), t ∈ [0, 1]. Для упрощения реализации дана выборка f(t) с шагом ∆t (0 < ∆t ≪ 1) и получаем последовательность 3d-точек f(0), f(∆t), ... , f(1). Каждая пара последовательных точек в этой последовательности может рассматриваться как начальная и конечная точки линейного движения сферы. Точки, пересекающиеся с движущейся сферой, считаются удаленными.
Задание: Реализовать функцию, которая принимает входные параметры (∆s, O, nx, ny, nz, R, f(t), ∆t), моделирует удаление точек сетки, которые пересекаются с линейными перемещениями сферы, и выводит все оставшиеся точки, видимые сверху.

Я пока мало знаю о компьютерной графике и алгоритмах, которые в ней используются. Вот что я пытался сделать:

C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
// Bresenham's algorithm for drawing circle
std::vector<Point3> GetCircleBoundary(const Point3& center, double radius)
{
    std::vector<Point3> points;
 
    int x = 0;
    int y = radius;
    int delta = 1 - 2 * radius;
    int error = 0;
    while (y >= x)
    {
        points.emplace_back(center.x() + x, center.y() + y, center.z());
        points.emplace_back(center.x() + x, center.y() - y, center.z());
        points.emplace_back(center.x() - x, center.y() + y, center.z());
        points.emplace_back(center.x() - x, center.y() - y, center.z());
        points.emplace_back(center.x() + y, center.y() + x, center.z());
        points.emplace_back(center.x() + y, center.y() - x, center.z());
        points.emplace_back(center.x() - y, center.y() + x, center.z());
        points.emplace_back(center.x() - y, center.y() - x, center.z());
        error = 2 * (delta + y) - 1;
        if ((delta < 0) && (error <= 0))
        {
            delta += 2 * ++x + 1;
            continue;
        }
        if ((delta > 0) && (error > 0))
        {
            delta -= 2 * --y + 1;
            continue;
        }
        delta += 2 * (++x - --y);
    }
    return points;
}
 
 
/// @param refPoint reference point O of the cloud, which is a point with the minimum values along
/// all coordinate axes
/// @param nx number of points in cloud along x axis
/// @param ny number of points in cloud along y axis
/// @param nz number of points in cloud along z axis
/// @param deltaS distance between neighboring cloud points along x, y and z axis
/// @param sphereRadius radius R of the sphere
/// @param curve 3d curve that defines trajectory of the sphere
/// @param deltaT step size for 3d curve parameter
void Process(
    const Point3 refPoint,
    const int nx,
    const int ny,
    const int nz,
    const double deltaS,
    const double sphereRadius,
    const Curve& curve,
    const double deltaT)
{
    int maxZ = refPoint.z() + nz * deltaS;
    std::vector<std::vector<std::optional<Point3>>> topSkin(ny, std::vector<std::optional<Point3>>(nx));    
 
    for (int y = 0; y < ny; ++y)
        for (int x = 0; x < nx; ++x)
            topSkin[y][x] = Point3(refPoint.x() + x * deltaS, refPoint.y() + y * deltaS, maxZ);
 
    Point3 sphereCenter;
    for (double t = curve.GetBeginParameter(); t <= curve.GetEndParameter(); t += deltaT)
    {
        sphereCenter = curve.Evaluate(t);
 
        if (sphereCenter.z() + sphereRadius <= maxZ)
            continue;
 
        double r{};
        std::vector<Point3> points;
        for (double z = maxZ; z <= sphereCenter.z() + sphereRadius; z += deltaS)
        {
            r = std::sqrt(sphereRadius * sphereRadius - (z - sphereCenter.z()) * (z - sphereCenter.z()));
            points = GetCircleBoundary(Point3(sphereCenter.x(), sphereCenter.y(), z), r);
 
            for (const auto& point : points)
            {
                std::size_t y = (point.y() - refPoint.y()) / deltaS;
                std::size_t x = (point.x() - refPoint.x()) / deltaS;
                if (2 * sphereCenter.z() - point.z() > refPoint.z())
                    topSkin[y][x] = Point3(point.x(), point.y(), 2 * sphereCenter.z() - point.z());
                else
                    topSkin[y][x] = {};
            }
        }
    }
 
        // WriteToFile(topSkin);
}

Тут для каждого момента времени я нахожу каждую точку сферы, которая 'вылазит' выше верхней поверхности сетки, и говорю, что тогда сверху будет видна точка, симметричная этой относительно центра сферы.
Направьте на путь истинный, если можно, как это сделать получше?

@zayats80888 · 01.03.2024, 12:12

Сообщение от MeXaL

Я пока мало знаю о компьютерной графике и алгоритмах, которые в ней используются.

По-моему эта задача не имеет отношения к графики, она про геометрию.

Сообщение от MeXaL

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

Во-первых, судя по использованию алгоритма Брезенхема, "каждая точка сферы" далеко не каждая, а только с целочисленными координатами, поэтому такой подход не корректен.
Во-вторых, в два последовательных момента времени сфера может и не пересекать точку сетки, а вот в промежутке между ними по ходу двежения вполне может. Т.е. работать нужно не с точками кривой пути, а с отрезками.

Сообщение от MeXaL

Направьте на путь истинный, если можно, как это сделать получше?

Я бы решал так:
1) Все вычисления делал бы в системе координат сетки, чтоб координаты точек сетки были целочисленными с началом в опорной точке.
Перевёл бы в эту систему радиус сферы и координаты кривой пути.
2) Состояние сетки(наличие/отсутствие) моделировал бы массивом флагов (можно даже битовым массивом, если размеры большие), типа:

C++
1
2
3
4
5
6
7
8
class Grid {
public:
  Grid(int nx, int ny, int nz);
  bool hasPoint(int ix, int iy, int iz) const;
  void removePoint(int ix, int iy, int iz);
private:
  std::vector<bool> m_state;
};

3) Для каждого отрезка пути вычислял бы AABB, раширял бы его на R (т.е. вычислял бы AABB капсулы, "построенной" сферой при движении вдоль этого отрезка). Для каждой "не удаленной" (Grid::hasPoint) точки сетки, поподающей внутрь этой коробки (три цикла в диапазонах целочисленных "координат" AABB) вычислял бы расстояние (квадрат расстояния) до отрезка пути, и если оно меньше радиуса, "удалял" бы точку сетки (Grid::removePoint).
4) В конце останется только пройтись по "столбцам сетки сверху вниз", чтоб сформировать массив точек, видимых сверху.

@MeXaL · 01.03.2024, 18:52 **[ТС]**

Та да, я понял, что не пойми что придумал, и что надо с отрезками работать, а не с точками.

То есть, если я правильно понял, вы предлагаете проходиться вообще по всем (неудаленным) точкам сетки на каждом шаге и определять расстояние от каждой из них до отрезка, пройденного сферой.
Это то, что я в самом начале и думал сделать, но это должно быть очень затратно, наверно. Может есть какой-то лучше способ? Октодеревья это в ту сторону?

@zayats80888 · 01.03.2024, 21:08

Сообщение от MeXaL

То есть, если я правильно понял, вы предлагаете проходиться вообще по всем (неудаленным) точкам сетки на каждом шаге

Не по всем, а только по тем, которые потенциально могут быть затронуты (для этого и вычисляем AABB).

@MeXaL · 03.03.2024, 11:53 **[ТС]**

Извините, я что-то запутался, что вы подразумеваете под AABB?

@zayats80888 · 03.03.2024, 12:19

Сообщение от MeXaL

Извините, я что-то запутался, что вы подразумеваете под AABB?

Axis Aligned Bounding Box

@MeXaL · 03.03.2024, 14:38 **[ТС]**

Что я не так делаю?

C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
class PointCloud
{
public:
    PointCloud(int nx, int ny, int nz)
        : nx(nx), ny(ny), nz(nz)
        , state(nx * ny * nz, false) {}
 
    bool hasPoint(int ix, int iy, int iz) const
    {
        return !state[ix + iy * ny + iz * nx * ny];
    }
 
    void removePoint(int ix, int iy, int iz)
    {
        state[ix + iy * ny + iz * nx * ny] = true;
    }
 
    int x() const noexcept { return nx; }
    int y() const noexcept { return nx; }
    int z() const noexcept { return nx; }
 
private:
    std::vector<bool> state;
    int nx{}, ny{}, nz{};
};
 
struct Rect
{
    Point3 point;
    Point3 size;
};
 
Rect GetAABB(const Point3& src, const Point3& dest, double sphereRadius)
{
    Point3 point;
    for (int i = 0; i < 3; ++i)
        point[i] = (src[i] < dest[i] ? src[i]  : dest[i]) - sphereRadius;
 
    Point3 size;
    for (int i = 0; i < 3; ++i)
        size[i] = std::abs(src[i] - dest[i]) + 2 * sphereRadius;
 
    return Rect{point, size};
}
 
Point3 CrossProduct(const Point3& a, const Point3& b) noexcept
{
    geo::Point3D result;
    result[0] = a[1] * b[2] - a[2] * b[1];
    result[1] = a[2] * b[0] - a[0] * b[2];
    result[2] = a[0] * b[1] - a[1] * b[0];
    return result;
}
 
double DotProduct(const Point3& a, const Point3& b) noexcept
{
    return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
}
 
// returns squared distance from point to vector (dest - src)
double Distance(const Point3& src, const Point3& dest, const Point3& point)
{
    if (DotProduct(dest - src, point - src) < 0)
    {
        double d1 = (point - src) .Length2();
        double d2 = (point - dest).Length2();
        return d1 < d2 ? d1 : d2;
    }
    return (CrossProduct(dest - src, point - src) / ~(dest - src)).Length2();
}
 
 
 
void Calculate(
    const Point3& refPoint,
    const int nx,
    const int ny,
    const int nz,
    const double deltaS,
    const double sphereRadius,
    const Curve& curve,
    const double deltaT,
    const std::filesystem::path& outputFileName)
{
    io::TestOutput to(outputFileName);
 
    PointCloud pointCloud(nx, ny, nz);
 
    double t = curve.GetBeginParameter();
    Point3 src = curve.Evaluate(t) - refPoint;
    t += deltaT;
 
    Point3 dest;
    double distance{};
 
    for (; t <= curve.GetEndParameter(); t += deltaT)
    {
        dest = curve.Evaluate(t) - refPoint;
        auto [point, size] = GetAABB(src, dest, sphereRadius);
 
        for (int iz = point.z(); iz < point.z() + size.z(); ++iz)
            for (int iy = point.y(); iy < point.y() + size.y(); ++iy)
                for (int ix = point.x(); ix < point.x() + size.x(); ++ix)
                {
                    if (pointCloud.hasPoint(ix, iy, iz))
                    {
                        distance = Distance(src, dest, Point3{ix * deltaS, iy * deltaS, iz * deltaS});
                        if (distance <= sphereRadius * sphereRadius)
                            pointCloud.removePoint(ix, iy, iz);
                    }
                }
 
        src = dest;
    }
 
 
    // output resulting points
    for (int iy = 0; iy < ny; ++iy)
        for (int ix = 0; ix < nx; ++ix)
            for (int iz = nz - 1; iz >= 0; --iz)
            {
                if (pointCloud.hasPoint(ix, iy, iz))
                {
                    to.Write(Point3{ix * deltaS, iy * deltaS, iz * deltaS} + refPoint);
                    break;
                }
            }
}

@zayats80888 · 03.03.2024, 15:11

Сообщение от MeXaL

Что я не так делаю?

1) Все ваши вычисления не в системе координат сетки, вы забыли про масштаб (1 / deltaS), поэтому некорректно целочисленные координаты AABB сопоставлять с координатами точек сетки (последние в вашей системе координат не целочисленные)
2) Ваша функция нахождения расстояния от точки до отрезка некорректная, псевдокод:

Code
1
2
3
4
5
6
7
8
9
10
// squared distance between point c and segment ab
double Distance(Point a, Point b, Point c)
{
    Point ab = b – a, ac = c – a, bc = c – b;
    double e = Dot (ac, ab);
    if (e <= 0.0f) return Dot(ac, ac);
    double f = Dot (ab, ab);
    if (e >= f) return Dot(bc, bc);
    return Dot(ac, ac) – e * e / f
}

3) С AABB вы тоже намудрили (имхо, удобнее представление не {min, size}, а {min, max})

@MeXaL · 03.03.2024, 19:35 **[ТС]**

Уже все стало лучше

C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
class PointCloud
{
public:
    PointCloud(int nx, int ny, int nz)
        : nx(nx), ny(ny), nz(nz)
        , m_state(nx * ny * nz, false) {}
 
    bool hasPoint(int ix, int iy, int iz) const
    {
        return !m_state[ix + iy * ny + iz * nx * ny];
    }
 
    void removePoint(int ix, int iy, int iz)
    {
        m_state[ix + iy * ny + iz * nx * ny] = true;
    }
 
    int x() const noexcept { return nx; }
    int y() const noexcept { return ny; }
    int z() const noexcept { return nz; }
 
private:
    std::vector<bool> m_state;
    int nx{}, ny{}, nz{};
};
 
struct BoundingBox
{
    Point3 min;
    Point3 max;
};
 
BoundingBox GetBoundingBox(const Point3& a, const Point3& b, double sphereRadius, const Point3& bounds)
{
    Point3 min, max;
    for (int i = 0; i < 3; ++i)
    {
        min[i] = std::min(a[i], b[i]) - sphereRadius;
        if (min[i] < 0)
            min[i] = 0;
        else if (min[i] >= bounds[i])
            min[i] = bounds[i] - 1;
 
        max[i] = std::max(a[i], b[i]) + sphereRadius;
        if (max[i] < 0)
            max[i] = 0;
        else if (max[i] >= bounds[i])
            max[i] = bounds[i] - 1;
    }
    return BoundingBox{min, max};
}
 
double DotProduct(const Point3& a, const Point3& b) noexcept
{
    return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
}
 
// returns squared distance from c to vector ab
double Distance(const Point3& a, const Point3& b, const Point3& c)
{
    Point3 ab = b - a, ac = c - a, bc = c - b;
    double e = DotProduct(ac, ab);
    if (e <= 0.0f) return DotProduct(ac, ac);
    double f = DotProduct(ab, ab);
    if (e >= f) return DotProduct(bc, bc);
    return DotProduct(ac, ac) - e * e / f;
}
 
class CoordConverter
{
public:
    CoordConverter(const Point3 refPoint = {}, double factor = 1.0)
        : m_refPoint(refPoint)
        , m_factor(factor)
    {
        if (factor == 0.0)
            throw std::runtime_error("Factor is equal to 0");
    }
 
    const Point3& getRefPoint() const noexcept { return m_refPoint; }
    double getFactor() const noexcept { return m_factor; }
 
    void setRefPoint(const Point3& point) noexcept { m_refPoint = point; }
    void setFactor(double factor)
    {
        if (factor == 0.0)
            throw std::runtime_error("Factor is equal to 0");
        m_factor = factor;
    }
 
    Point3 rel(const Point3& point) const noexcept
    {
        return (point - m_refPoint) / m_factor;
    }
 
    double rel(double distance) const noexcept
    {
        return distance / m_factor;
    }
 
    Point3 abs(const Point3& point) const noexcept
    {
        return point * m_factor + m_refPoint;
    }
 
    double abs(double distance) const noexcept
    {
        return distance * m_factor;
    }
 
private:
    Point3 m_refPoint;
    double m_factor = 1.0;
};
 
 
void Calculate(
    const Point3& refPoint,
    const int nx,
    const int ny,
    const int nz,
    const double deltaS,
    const double sphereRadius,
    const Curve& curve,
    const double deltaT,
    const std::filesystem::path& outputFileName)
{
    io::TestOutput to(outputFileName);
 
    PointCloud pointCloud(nx, ny, nz);
 
    CoordConverter conv(refPoint, deltaS);
    double relSphereRadius = conv.rel(sphereRadius);
 
    double t = curve.GetBeginParameter();
    Point3 a = conv.rel(curve.Evaluate(t));
    t += deltaT;
 
    Point3 b;
    double distance{};
 
    for (; t <= curve.GetEndParameter(); t += deltaT)
    {
        b = conv.rel(curve.Evaluate(t));            
        
        auto [min, max] = GetBoundingBox(a, b, relSphereRadius, Point3{nx, ny, nz});
 
        for (int iz = min.z(); iz <= max.z(); ++iz)
            for (int iy = min.y(); iy <= max.y(); ++iy)
                for (int ix = min.x(); ix <= max.x(); ++ix)
                {
                    if (pointCloud.hasPoint(ix, iy, iz))
                    {
                        distance = Distance(a, b, Point3{ix, iy, iz});
                        if (distance <= std::pow(relSphereRadius, 2))
                            pointCloud.removePoint(ix, iy, iz);
                    }
                }
 
        a = b;
    }
 
 
 
    // TODO: output resulting points
    for (int iy = 0; iy < ny; ++iy)
        for (int ix = 0; ix < nx; ++ix)
            for (int iz = nz - 1; iz >= 0; --iz)
            {
                if (pointCloud.hasPoint(ix, iy, iz))
                {
                    to.Write(conv.abs(Point3{ix, iy, iz}));
                    break;
                }
            }
}

Но иногда вот такое возникает (3 полоски, а должна быть одна)

@zayats80888 · 03.03.2024, 21:00

Сообщение от MeXaL

Но иногда вот такое возникает (3 полоски, а должна быть одна)

Трудно сказать почему, но похожее бывает при выходе за пределы в одномерном представлении многомерного массива.
Попробуйте ассертов наставить в коде.

Кликните здесь для просмотра всего текста

А еще мне не нравиться, как вы по "диапазону" AABB итерируете, лучше преобразовать min и max в int перед циклом, причем для min использовать ceil, а для max floor, так будет точнее.

Добавлено через 2 минуты
И еще вопрос, а какие значения nx, ny и nz? Переполнения в выражении

Сообщение от MeXaL

nx * ny * nz

случайно нет?

@MeXaL · 03.03.2024, 21:25 **[ТС]**

Спасибо вам большое за оказанную помощь!! Дело было действительно в переполнении.

@Igor3D · 04.03.2024, 14:33

Не стоит связываться с 3-мерным массивом (иногда называют lattice) без необходимости. Здесь можно обойтись 2-мерным. И удаляются все точки внутри (хотя бы одной) сферы - а они могут быть и не видны сверху. В общем, слишком прямолинейно

Не по теме:

Зато как владеет "современным С++" :)

@MeXaL · 04.03.2024, 15:12 **[ТС]**

Спасибо за замечание. Изучил C++, теперь буду изучать векторную алгебру и компьютерную графику. Планирую пойти в геймдев или разработку CAD-, CAM- систем.

@zayats80888 · 04.03.2024, 15:20

Сообщение от Igor3D

Здесь можно обойтись 2-мерным.

Можно, но тогда придется повторять расчеты по траектории до тех пор, пока массив не перестанет меняться. Или же использовать представление "2-мерный массив дырок", но тогда может случится перерасход памяти в сравнении с 3-мерным массивом битов. Все зависит от траектории и размеров сферы относительно сетки. Тут как с сортировкой, где-то лучше "расческой", а где-то "вставками" эффективней, зависит от входных данных.

Не по теме:

Сообщение от Igor3D

Зато как владеет "современным С++" :)

Что там современного увидели? Structured binding уже лет 7 как ввели, noexcept уже лет 13. Если это убрать(ну и filesystem тоже), то в с++98 скомпилируется.

@Igor3D · 04.03.2024, 16:53

Сообщение от zayats80888

Можно, но тогда придется повторять расчеты по траектории до тех пор, пока массив не перестанет меняться.

Этого не понял

Сообщение от zayats80888

Или же использовать представление "2-мерный массив дырок", но тогда может случится перерасход памяти в сравнении с 3-мерным массивом битов.

Если требуется только "вид сверху" - вообще не вопрос. Если же в общем виде нужен как бы кусок сыра изрытый червями/кротами - то и хранить для каждого элемента 2-мерного массива пару (контейнер пар) индексов "с какого по какой" удалено.

Не по теме:

Сообщение от zayats80888

Что там современного увидели?

К сожалению, увлечение всеми std::новшествами часто идет в ущерб содержательной части :(

@zayats80888 · 04.03.2024, 17:02

Сообщение от Igor3D

Этого не понял

Ну это в случае, если храним только "слой, видимый сверху", а не разреженное представления 3-мерного массива.
В этом случае мы игнорируем "дырки" под этим слоем, но при самопересечении траетории эти дырки могут "вскрыться", а чтобы их "увидеть", придется повторить расчет. Поэтому рассчет придеться повторять до тех пор, пока видимый слой не перестанет меняться.

@Igor3D · 04.03.2024, 18:52

Сообщение от zayats80888

Ну это в случае, если храним только "слой, видимый сверху", а не разреженное представления 3-мерного массива.
В этом случае мы игнорируем "дырки" под этим слоем, но при самопересечении траетории эти дырки могут "вскрыться", а чтобы их "увидеть", придется повторить расчет. Поэтому рассчет придеться повторять до тех пор, пока видимый слой не перестанет меняться.

Теперь понял. Да, придется хранить "дырки в глубине" для возможного объединения. С др стороны, битовый массив может оказаться не очень удобным (напр при отрисовке того же "вида сверху")

@MeXaL · 07.03.2024, 13:18 **[ТС]**

Друзья, помогите еще раз, времени мало осталось, а у меня не выходит. Мне сказали оптимизировать решение по памяти. Я подумал, что можно по-другому делать: хранить массив 2x2 z-координат видимых сверху точек. И тогда для каждого отрезка пути нужно находить часть bounding-box'а, которая 'вылазит' выше чем сетка, и для каждой попадающей внутрь сферы точки искать симметричную ей.

Вот это я написал. Подскажите где тут что не так, на тестах не все точки попадают.
Тут я для каждой точки ищу где в этот момент был центр сферы на отрезке, нахожу z-координату симметричную относительно центра и говорю, что ее сверху будет видно.

C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
class PointCloud
{
public:
    static constexpr int NOT_VISIBLE = -1;
 
    PointCloud(int nx, int ny, int nz)
            : nx(nx), ny(ny), nz(nz)
            , m_state(ny, std::vector<int>(nx, nz - 1))
    {}
 
    bool IsVisibleFromAbove(int ix, int iy) const
    {
        return m_state[iy][ix] != NOT_VISIBLE;
    }
 
    int GetVisibleFromAbove(int ix, int iy) const
    {
        return m_state[iy][ix];
    }
 
    void SetVisibleFromAbove(int ix, int iy, int iz)
    {
        if (iz < 0)
            m_state[iy][ix] = NOT_VISIBLE;
        else if (iz < m_state[iy][ix])
            m_state[iy][ix] = iz;
    }
 
    int x() const noexcept { return nx; }
    int y() const noexcept { return ny; }
    int z() const noexcept { return nz; }
 
private:
    int nx{}, ny{}, nz{};
    std::vector<std::vector<int>> m_state;
};
 
struct BoundingBox
{
    geo::Point3<int> min;
    geo::Point3<int> max;
 
    BoundingBox(const geo::Point3<int>& min = {}, const geo::Point3<int>& max = {}) noexcept
            : min(min)
            , max(max)
    {}
 
    bool IntersectPlane(int z) const noexcept
    {
        return z >= min.z() && z <= max.z();
    }
 
    void CutZ(int z) noexcept
    {
        min.z(z);
    }
 
    void CutXY(const geo::Point3<int>& bounds) noexcept
    {
        for (int i = 0; i < 2; ++i)
        {
            if (min[i] < 0)
                min[i] = 0;
            else if (min[i] >= bounds[i])
                min[i] = bounds[i] - 1;
 
            if (max[i] < 0)
                max[i] = 0;
            else if (max[i] >= bounds[i])
                max[i] = bounds[i] - 1;
        }
    }
};
 
// Find bounding box of the capsule built by the sphere as it moves along the segment ab
BoundingBox GetBoundingBox(const geo::Point3D& a, const geo::Point3D& b, double sphereRadius)
{
    geo::Point3<int> min, max;
    for (int i = 0; i < 3; ++i)
    {
        min[i] = std::ceil(std::min(a[i], b[i]) - sphereRadius);
        max[i] = std::floor(std::max(a[i], b[i]) + sphereRadius);
    }
    return BoundingBox{min, max};
}
 
double DotProduct(const geo::Point3D& a, const geo::Point3D& b) noexcept
{
    return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
}
 
// Calculate squared distance from point c to segment ab
double Distance(const geo::Point3D& a, const geo::Point3D& b, const geo::Point3D& c)
{
    geo::Point3D ab = b - a, ac = c - a, bc = c - b;
    double e = DotProduct(ac, ab);
    if (e <= 0.0f) return ac.Length2();
    double f = ab.Length2();
    if (e >= f) return bc.Length2();
    return ac.Length2() - e * e / f;
}
 
geo::Point3D FindCenter(const geo::Point3D& a, const geo::Point3D& b, const geo::Point3D& c)
{
    geo::Point3D ab = b - a, ac = c - a, bc = c - b;
    double e = DotProduct(ac, ab);
    if (e <= 0.0f) return a;
    double f = ab.Length2();
    if (e >= f) return b;
    return a + ab * e / f;
}
 
double FindSymmetric(double point_z, double center_z) noexcept
{
    return 2 * center_z - point_z;
}
 
class CoordConverter
{
public:
    CoordConverter(const geo::Point3D& refPoint = {}, double factor = 1.0)
        : m_refPoint(refPoint)
    {
        SetFactor(factor);
    }
 
    const geo::Point3D& GetRefPoint() const noexcept { return m_refPoint; }
    double GetFactor() const noexcept { return m_factor; }
 
    void SetRefPoint(const geo::Point3D& point) noexcept { m_refPoint = point; }
    void SetFactor(double factor)
    {
        CheckFactor(factor);
        m_factor = factor;
    }
 
    geo::Point3D ToRel(const geo::Point3D& point) const noexcept
    {
        return (point - m_refPoint) / m_factor;
    }
 
    double ToRel(double distance) const noexcept
    {
        return distance / m_factor;
    }
 
    geo::Point3D ToAbs(const geo::Point3D& point) const noexcept
    {
        return point * m_factor + m_refPoint;
    }
 
    double ToAbs(double distance) const noexcept
    {
        return distance * m_factor;
    }
 
private:
    geo::Point3D m_refPoint;
    double m_factor = 1.0;
 
    static void CheckFactor(double factor)
    {
        if (factor == 0.0)
            throw std::runtime_error("Factor is equal to 0");
    }
};
 
 
void Calculate(
    const geo::Point3D& refPoint,
    const int nx,
    const int ny,
    const int nz,
    const double deltaS,
    const double sphereRadius,
    const geo::Curve& curve,
    const double deltaT,
    const std::filesystem::path& outputFileName)
{
    io::TestOutput to(outputFileName);
 
    PointCloud pointCloud(nx, ny, nz);
 
    CoordConverter conv(refPoint, deltaS);
    double relSphereRadius = conv.ToRel(sphereRadius);
 
    double t = curve.GetBeginParameter();
    geo::Point3D a = conv.ToRel(curve.Evaluate(t));
    t += deltaT;
    geo::Point3D b;
    geo::Point3<int> bounds(nx, ny, nz);
    for (; t <= curve.GetEndParameter(); t += deltaT)
    {
        b = conv.ToRel(curve.Evaluate(t));          
        
        BoundingBox box = GetBoundingBox(a, b, relSphereRadius);
 
        if (box.IntersectPlane(bounds.z() - 1))
        {
            box.CutZ(bounds.z() - 1);
            box.CutXY(bounds);
            for (int iz = box.min.z(); iz <= box.max.z(); ++iz)
                for (int iy = box.min.y(); iy <= box.max.y(); ++iy)
                    for (int ix = box.min.x(); ix <= box.max.x(); ++ix)
                    {
                        geo::Point3D point(static_cast<double>(ix), static_cast<double>(iy), static_cast<double>(iz));
                        if (Distance(a, b, point) <= relSphereRadius * relSphereRadius)
                        {
                            geo::Point3D center = FindCenter(a, b, point);
                            double iz = FindSymmetric(point.z(), center.z());
                            pointCloud.SetVisibleFromAbove(ix, iy, iz);
                        }
                    }
        }
 
        a = b;
    }
 
    for (int iy = 0; iy < ny; ++iy)
        for (int ix = 0; ix < nx; ++ix)
        {
            if (pointCloud.IsVisibleFromAbove(ix, iy))
            {
                int iz = pointCloud.GetVisibleFromAbove(ix, iy);
                geo::Point3D point(static_cast<double>(ix), static_cast<double>(iy), static_cast<double>(iz));
                to.Write(conv.ToAbs(point));
            }
        }
}

@zayats80888 · 07.03.2024, 14:22

Сообщение от MeXaL

Я подумал, что можно по-другому делать: хранить массив 2x2 z-координат видимых сверху точек.

Я об этом и писал, как раз про то, что весь расчет придется повторять, пока массив не перестанет меняться.

Сообщение от MeXaL

И тогда для каждого отрезка пути нужно находить часть bounding-box'а, которая 'вылазит' выше чем сетка, и для каждой попадающей внутрь сферы точки искать симметричную ей.

Что-то как-то мудрено. В первом приближении: для каждой точки слоя, попадающей внутрь AABB, пока она лежит внутри капсулы(дистанция до отрезка, как и прежде) уменьшаем её z координату на 1.

Добавлено через 17 минут
MeXaL, псевдокод(может так будет понятнее)

Code
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
for(bool isChanged = true; isChanged;)
{
  isChanged = false;
  foreach(segment in path)
  {
    box = aabb(segment, radius, gridBounds)
    foreach(y in box.min.y..box.max.y)
    {
      foreach(x in box.min.x..box.max.x)
      {
        z = grid[x][y];
        if (z in box.min.z..box.max.z)
        {
          while (sqDistance(segment, point(x, y, z) <= radius * radius))
          {
            --z;
            isChanged = true;
            grid[x][y] = z;
          }
        }
      }
    }
  }
}

@MeXaL · 07.03.2024, 14:46 **[ТС]**

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

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

C++
1
2
3
4
5
6
7
8
9
std::pair<double, geo::Point3D> Distance(const geo::Point3D& a, const geo::Point3D& b, const geo::Point3D& c) noexcept
{
    geo::Point3D ab = b - a, ac = c - a, bc = c - b;
    double ac_ab = DotProduct(ac, ab);
    if (ac_ab <= 0.0f) return { ac.Length2(), a };
    double ab2 = ab.Length2();
    if (ac_ab >= ab2) return { bc.Length2(), b };
    return { ac.Length2() - ac_ab * ac_ab / ab2, a + ab * ac_ab / ab2 };
}

А потом искать симметричную z-координату

Вот это вот в целом почти работает, но кое-где дырки остаются

C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
class TopView
{
public:
    static constexpr int NOT_VISIBLE = -1;
 
    TopView(int nx, int ny, int nz)
            : nx(nx), ny(ny), nz(nz)
            , m_state(ny, std::vector<int>(nx, nz - 1))
    {}
 
    bool IsVisibleFromAbove(int ix, int iy) const
    {
        return m_state[iy][ix] != NOT_VISIBLE;
    }
 
    int GetVisibleFromAbove(int ix, int iy) const
    {
        return m_state[iy][ix];
    }
 
    void SetVisibleFromAbove(int ix, int iy, int iz)
    {
        if (iz < 0)
            m_state[iy][ix] = NOT_VISIBLE;
        else if (iz < m_state[iy][ix])
            m_state[iy][ix] = iz;
    }
 
    int x() const noexcept { return nx; }
    int y() const noexcept { return ny; }
    int z() const noexcept { return nz; }
 
private:
    int nx{}, ny{}, nz{};
    std::vector<std::vector<int>> m_state;
};
 
struct BoundingBox
{
    geo::Point3<int> min;
    geo::Point3<int> max;
 
    BoundingBox(const geo::Point3<int>& min = {}, const geo::Point3<int>& max = {}) noexcept
            : min(min)
            , max(max)
    {}
 
    bool IntersectPlane(int z) const noexcept
    {
        return z >= min.z() && z <= max.z();
    }
 
    void CutZ(int z) noexcept
    {
        min.z(z);
    }
 
    void CutXY(const geo::Point3<int>& bounds) noexcept
    {
        for (int i = 0; i < 2; ++i)
        {
            if (min[i] < 0)
                min[i] = 0;
            else if (min[i] >= bounds[i])
                min[i] = bounds[i] - 1;
 
            if (max[i] < 0)
                max[i] = 0;
            else if (max[i] >= bounds[i])
                max[i] = bounds[i] - 1;
        }
    }
};
 
// Find bounding box of the capsule built by the sphere as it moves along the segment ab
BoundingBox GetBoundingBox(const geo::Point3D& a, const geo::Point3D& b, double sphereRadius) noexcept
{
    geo::Point3<int> min, max;
    for (int i = 0; i < 3; ++i)
    {
        min[i] = std::ceil(std::min(a[i], b[i]) - sphereRadius);
        max[i] = std::floor(std::max(a[i], b[i]) + sphereRadius);
    }
    return BoundingBox(min, max);
}
 
double DotProduct(const geo::Point3D& a, const geo::Point3D& b) noexcept
{
    return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
}
 
std::pair<double, geo::Point3D> Distance(const geo::Point3D& a, const geo::Point3D& b, const geo::Point3D& c) noexcept
{
    geo::Point3D ab = b - a, ac = c - a, bc = c - b;
    double ac_ab = DotProduct(ac, ab);
    if (ac_ab <= 0.0f) return { ac.Length2(), a };
    double ab2 = ab.Length2();
    if (ac_ab >= ab2) return { bc.Length2(), b };
    return { ac.Length2() - ac_ab * ac_ab / ab2, a + ab * ac_ab / ab2 };
}
 
double FindSymmetricZ(double pointZ, double centerZ) noexcept
{
    return 2 * centerZ - pointZ;
}
 
geo::Point3D GetPoint3D(const geo::Point3<int>& point) noexcept
{
    return geo::Point3D(static_cast<double>(point.x()), static_cast<double>(point.y()), static_cast<double>(point.z()));
}
 
geo::Point3D GetPoint3D(int x, int y, int z) noexcept
{
    return geo::Point3D(static_cast<double>(x), static_cast<double>(y), static_cast<double>(z));
}
 
class CoordConverter
{
public:
    CoordConverter(const geo::Point3D& refPoint = {}, double factor = 1.0)
        : m_refPoint(refPoint)
    {
        SetFactor(factor);
    }
 
    const geo::Point3D& GetRefPoint() const noexcept { return m_refPoint; }
    double GetFactor() const noexcept { return m_factor; }
 
    void SetRefPoint(const geo::Point3D& point) noexcept { m_refPoint = point; }
    void SetFactor(double factor)
    {
        CheckFactor(factor);
        m_factor = factor;
    }
 
    geo::Point3D ToRel(const geo::Point3D& point) const noexcept
    {
        return (point - m_refPoint) / m_factor;
    }
 
    double ToRel(double distance) const noexcept
    {
        return distance / m_factor;
    }
 
    geo::Point3D ToAbs(const geo::Point3D& point) const noexcept
    {
        return point * m_factor + m_refPoint;
    }
 
    double ToAbs(double distance) const noexcept
    {
        return distance * m_factor;
    }
 
private:
    geo::Point3D m_refPoint;
    double m_factor = 1.0;
 
    static void CheckFactor(double factor)
    {
        if (factor == 0.0)
            throw std::runtime_error("Factor is equal to 0");
    }
};
 
 
/// Calculate the result and output it to the outputFileName
///
/// @param refPoint reference point O of the cloud, which is a point with the minimum values along
/// all coordinate axes
/// @param nx number of points in cloud along x axis
/// @param ny number of points in cloud along y axis
/// @param nz number of points in cloud along z axis
/// @param deltaS distance between neighboring cloud points along x, y and z axis
/// @param sphereRadius radius R of the sphere
/// @param curve 3d curve that defines trajectory of the sphere
/// @param deltaT step size for 3d curve parameter
/// @param outputFileName name of the output file with result
void Calculate(
    const geo::Point3D& refPoint,
    const int nx,
    const int ny,
    const int nz,
    const double deltaS,
    const double sphereRadius,
    const geo::Curve& curve,
    const double deltaT,
    const std::filesystem::path& outputFileName)
{
    io::TestOutput to(outputFileName);
 
    TopView topView(nx, ny, nz);
 
    CoordConverter conv(refPoint, deltaS);
    double relSphereRadius = conv.ToRel(sphereRadius);
 
    double t = curve.GetBeginParameter();
    geo::Point3D a = conv.ToRel(curve.Evaluate(t)), b;
    t += deltaT;
 
    for (; t <= curve.GetEndParameter(); t += deltaT)
    {
        b = conv.ToRel(curve.Evaluate(t));          
        
        BoundingBox box = GetBoundingBox(a, b, relSphereRadius);
        if (box.IntersectPlane(nz - 1))
        {
            box.CutZ(nz - 1);
            box.CutXY(geo::Point3<int>(nx, ny, nz));
            for (int iz = box.min.z(); iz <= box.max.z(); ++iz)
                for (int iy = box.min.y(); iy <= box.max.y(); ++iy)
                    for (int ix = box.min.x(); ix <= box.max.x(); ++ix)
                    {
                        geo::Point3D point = GetPoint3D(ix, iy, iz);
                        auto [distance2, sphereCenter] = Distance(a, b, point);
                        if (distance2 <= relSphereRadius * relSphereRadius)
                        {
                            double iz = FindSymmetricZ(point.z(), sphereCenter.z());
                            topView.SetVisibleFromAbove(ix, iy, iz - 1);
                        }
                    }
        }
 
        a = b;
    }
 
    for (int iy = 0; iy < ny; ++iy)
        for (int ix = 0; ix < nx; ++ix)
        {
            if (topView.IsVisibleFromAbove(ix, iy))
            {
                int iz = topView.GetVisibleFromAbove(ix, iy);
                geo::Point3D point = GetPoint3D(ix, iy, iz);
                to.Write(conv.ToAbs(point));
            }
        }
}

Новые блоги и статьи Все статьи Все блоги /
Новый ноутбук volvo 07.12.2025 Всем привет. По скидке в "черную пятницу" взял себе новый ноутбук Lenovo ThinkBook 16 G7 на Амазоне: Ryzen 5 7533HS 64 Gb DDR5 1Tb NVMe 16" Full HD Display Win11 Pro	Музыка, написанная Искусственным Интеллектом volvo 04.12.2025 Всем привет. Некоторое время назад меня заинтересовало, что уже умеет ИИ в плане написания музыки для песен, и, собственно, исполнения этих самых песен. Стихов у нас много, уже вышли 4 книги, еще 3. . .	От async/await к виртуальным потокам в Python IndentationError 23.11.2025 Армин Ронахер поставил под сомнение async/ await. Создатель Flask заявляет: цветные функции - провал, виртуальные потоки - решение. Не threading-динозавры, а новое поколение лёгких потоков. Откат?. . .	Поиск "дружественных имён" СОМ портов Argus19 22.11.2025 Поиск "дружественных имён" СОМ портов На странице: https:/ / norseev. ru/ 2018/ 01/ 04/ comportlist_windows/ нашёл схожую тему. Там приведён код на С++, который показывает только имена СОМ портов, типа,. . .	Сколько Государство потратило денег на меня, обеспечивая инсулином. Programma_Boinc 20.11.2025 Сколько Государство потратило денег на меня, обеспечивая инсулином. Вот решила сделать интересный приблизительный подсчет, сколько государство потратило на меня денег на покупку инсулинов. . . .
Ломающие изменения в C#.NStar Alpha Etyuhibosecyu 20.11.2025 Уже можно не только тестировать, но и пользоваться C#. NStar - писать оконные приложения, содержащие надписи, кнопки, текстовые поля и даже изображения, например, моя игра "Три в ряд" написана на этом. . .	Мысли в слух kumehtar 18.11.2025 Кстати, совсем недавно имел разговор на тему медитаций с людьми. И обнаружил, что они вообще не понимают что такое медитация и зачем она нужна. Самые базовые вещи. Для них это - когда просто люди. . .	Создание Single Page Application на фреймах krapotkin 16.11.2025 Статья исключительно для начинающих. Подходы оригинальностью не блещут. В век Веб все очень привыкли к дизайну Single-Page-Application . Быстренько разберем подход "на фреймах". Мы делаем одну. . .	Фото: Daniel Greenwood kumehtar 13.11.2025	Расскажи мне о Мире, бродяга kumehtar 12.11.2025 — Расскажи мне о Мире, бродяга, Ты же видел моря и метели. Как сменялись короны и стяги, Как эпохи стрелою летели. - Этот мир — это крылья и горы, Снег и пламя, любовь и тревоги, И бескрайние. . .

@MeXaL 2 / 2 / 0 Регистрация: 10.09.2023 Сообщений: 70
	01.03.2024, 18:52 [ТС]
	Та да, я понял, что не пойми что придумал, и что надо с отрезками работать, а не с точками. То есть, если я правильно понял, вы предлагаете проходиться вообще по всем (неудаленным) точкам сетки на каждом шаге и определять расстояние от каждой из них до отрезка, пройденного сферой. Это то, что я в самом начале и думал сделать, но это должно быть очень затратно, наверно. Может есть какой-то лучше способ? Октодеревья это в ту сторону? 0

@MeXaL 2 / 2 / 0 Регистрация: 10.09.2023 Сообщений: 70
	03.03.2024, 11:53 [ТС]
	Извините, я что-то запутался, что вы подразумеваете под AABB? 0

@MeXaL 2 / 2 / 0 Регистрация: 10.09.2023 Сообщений: 70
	03.03.2024, 21:25 [ТС]
	Спасибо вам большое за оказанную помощь!! Дело было действительно в переполнении. 1

@Igor3D 1958 / 814 / 114 Регистрация: 01.10.2012 Сообщений: 4,743 Записей в блоге: 2
	04.03.2024, 14:33
	Не стоит связываться с 3-мерным массивом (иногда называют lattice) без необходимости. Здесь можно обойтись 2-мерным. И удаляются все точки внутри (хотя бы одной) сферы - а они могут быть и не видны сверху. В общем, слишком прямолинейно Не по теме: Зато как владеет "современным С++" :) 1

@MeXaL 2 / 2 / 0 Регистрация: 10.09.2023 Сообщений: 70
	04.03.2024, 15:12 [ТС]
	Спасибо за замечание. Изучил C++, теперь буду изучать векторную алгебру и компьютерную графику. Планирую пойти в геймдев или разработку CAD-, CAM- систем. 0

Пересечение трехмерной сетки точек движущейся сферой

Решение