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

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