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

// 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);
}

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
            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;
};
 
 
TopView topView(nx, ny, nz);
bool isChanged = true;
while (isChanged)
{
    isChanged = false;
    for (; t < curve.GetEndParameter(); t += deltaT)
    {
        b = conv.ToRel(curve.Evaluate(t));
        auto [min, max] = GetBoundingBox(a, b, relSphereRadius, bounds);
        for (int y = min.y(); y <= max.y(); ++y)
        {
            for (int x = min.x(); x <= max.x(); ++x)
            {
                int z = topView.GetVisibleFromAbove(x, y);
                if (z >= min.z() && z <= max.z())
                {
                    while (Distance(a, b, geo::Point3D(x, y, z)) <= relSphereRadius * relSphereRadius)
                    {
                        --z;
                        isChanged = true;
                        topView.SetVisibleFromAbove(x, y, z);
                    }
                }
            }
        }
    }
}