How to make a mesh surface that looks like cubes stacked in the shape of a sphere?

Question

I want to make a mesh object that looks like the outside surface of cubes stacked in the shape of a sphere. I thought I saw this once, but I can't find it.

The final size will be big (up to 100x100x100 cubes) so I think making a million cubes and merging them (modifier Boolean Union) is going to take quite a while, especially if I have to do it differently for every frame.

As a first step, I've stacked cubes in a hollow spherical shell.

What I really want is just a single mesh object of the outside surface. But for 100x100x100 that's a million cubes to Union. There must be an easier way? 2) if not, how can I select say 20,000 cubes in a "shell" and Union them at once?

Through brute-force application of logic & math & coffee, I guess one could define an algorithm to generate the mesh directly, but I somehow think that I saw this done before in an elegant way. "cubification"? "cubation"?

Here's an example of a small shell that takes about ten seconds to generate. In the picture I've removed the top.

import bpy
import numpy as np

nx, ny, nz = 21, 21, 21
x0, y0, z0 = 10, 10, 10
r1, r2 = 6.0, 7.4

scale = 0.5

x, y, z = np.arange(nx), np.arange(ny), np.arange(nz)

rsq = (x[None, None, :]-x0)**2 + (y[None, :, None]-y0)**2 + (z[:, None, None]-z0)**2
sphere_1 = rsq < r1**2
sphere_2 = rsq < r2**2
shell = sphere_2 * (-sphere_1)    # I love python!

X, Y, Z = np.meshgrid(x, y, z)

ipoints = zip( X.flatten(), 
               Y.flatten(), 
               Z.flatten() )

points = zip(scale*(X.flatten()-x0), 
             scale*(Y.flatten()-y0), 
             scale*(Z.flatten()-z0) ) 

points = [thing for thing in points]  # convert zip object to list (python 3!)
ipoints = [thing for thing in ipoints]

points_shell = [thing for i,thing in enumerate(points)
                 if shell[ipoints[i]]]

for point in points_shell:
    bpy.ops.mesh.primitive_cube_add(radius=0.5*scale, location=point)

Answer 1

I think you are on the right track here, but the slowest part is the creation of all the primitives. (reason - each new object's name must be checked for name collisions, and this gets progressively slower as more objects are added, it's not something we can prevent, via python)

Instead use from_pydata and generate the cube data as a function of the their center x,y,z and then extend the verts and face indices that way. Lateron you can do a post processing step which removes all the redundant (inner) faces

As related question / answer here: Extruding problem shows how to construct the pydata for cubes in one massive list

Here this is rather fast, but doesn't remove redundant geometry -- but it is miles faster than adding primitives: https://gist.github.com/zeffii/7ef17b4932e61ee12221

Answer 2

One possible way to do this would be to use the 'Remesh' modifier. It has an option to turn the selected object into blocks, though specifying the exact amount of blocks may be a bit difficult.

Here I have a UV Sphere with a remesh modifier set to 'Blocks'. I also added a Subsurf modifier before it so that the surface is smoother for the remesh:

You can use the 'Octree Depth' property to increase the resolution, but it will get increasingly slow when increasing this number.

Answer 3

Just to be nice:

import bpy
import numpy as np


def basic_cube(scale=None, offset=None, idx_offset=None):
    verts = [[-1.0, -1.0, -1.0], [-1.0, -1.0, 1.0], [-1.0, 1.0, -1.0], [-1.0, 1.0, 1.0], [1.0, -1.0, -1.0], [1.0, -1.0, 1.0], [1.0, 1.0, -1.0], [1.0, 1.0, 1.0]]
    edges = []
    faces = [[0, 1, 3, 2], [2, 3, 7, 6], [6, 7, 5, 4], [4, 5, 1, 0], [2, 6, 4, 0], [7, 3, 1, 5]]
    if idx_offset is not None:
        npf = np.array(faces) + (idx_offset * 8)
        faces = npf.tolist()
    if scale is not None:
        npv = np.array(verts) * scale
        if offset is not None:
            npv += offset
            verts = npv.tolist()
        
    return verts, edges, faces
    

def link_mesh(verts, edges, faces, name='name', ob=None):
    """Generate and link a new object from pydata.
    If object already exists replace its data
    with a new mesh and delete the old mesh."""
    if ob is None:
        mesh = bpy.data.meshes.new(name)
        mesh.from_pydata(verts, edges, faces)
        mesh.update()
        mesh_ob = bpy.data.objects.new(name, mesh)
        bpy.context.collection.objects.link(mesh_ob)
        return mesh_ob

    mesh_ob = ob
    old = ob.data
    mesh = bpy.data.meshes.new(name)
    mesh.from_pydata(verts, edges, faces)
    mesh.update()
    mesh_ob.data = mesh
    bpy.data.meshes.remove(old)
    return mesh_ob


def cubify():

    nx, ny, nz = 21, 21, 21
    x0, y0, z0 = 10, 10, 10
    r1, r2 = 6.0, 7.4

    scale = 0.5

    x, y, z = np.arange(nx), np.arange(ny), np.arange(nz)

    rsq = (x[None, None, :]-x0)**2 + (y[None, :, None]-y0)**2 + (z[:, None, None]-z0)**2
    sphere_1 = rsq < r1**2
    sphere_2 = rsq < r2**2
    shell = sphere_2 * (~sphere_1)    # I love python!

    X, Y, Z = np.meshgrid(x, y, z)

    ipoints = zip( X.ravel(), 
                   Y.ravel(), 
                   Z.ravel() )

    points = zip(scale*(X.ravel()-x0), 
                 scale*(Y.ravel()-y0), 
                 scale*(Z.ravel()-z0) ) 

    points = [thing for thing in points]  # convert zip object to list (python 3!)
    ipoints = [thing for thing in ipoints]

    points_shell = [thing for i,thing in enumerate(points)
                     if shell[ipoints[i]]]

    verts = []
    faces = []
    edges = []
    count = 0
    for point in points_shell:
        print(point)    
        v, e, f = basic_cube(scale=scale * 0.5, offset=np.array(point), idx_offset=count)
        verts += v
        faces += f   
        count += 1
        
    ob = link_mesh(verts, edges, faces, name='cubify', ob=None)