New function to generate CLUMP with uniform radius of sphere

[adrienluyckx]

Recently, I had the need to make flexible CLUMP composed of spheres of the same size with no overlap. In this manner they could break and have no problem of new volume created in the case of breakage.

I thus started from your GenerateClump_Euclidean_3D.py framework and expanded it to a "uniform_radius" version where spheres are placed from the center to the edges of the shape (allowing for small reduction in the radius of the sphere to fill better the volume).

Some comment might be useless but the main idea are there, you could maybe add a new function similar to the package.

Here is the code,
`def GenerateClump_Euclidean_3D_uniform_radius(inputGeom, N, rMin, div, **kwargs):
""" Function to generate clumps using the Euclidean map for voxelated, 3D particles
:param inputGeom: Directory of stl file, used to generate spheres
:param N: Number of spheres to be generated.
:param rMin: Minimum allowed radius: When this radius is met, the generation procedure stops even before N spheres are generated.
:param div: Division number along the shortest edge of the AABB during voxelisation (resolution). If not given, div=50 (default value in iso2mesh).
:param kwargs: Can contain either of the optional variables "output", "outputVTK", "visualise", "rMax_ratio".
- File name for output of the clump in .txt form. If not assigned, a .txt output file is not created.

 :Keyword Arguments:
     * *output*: txt file with centroids and radii, with format: [x,y,z,r] (optional)*

     * *outputVTK*: vtk file containing [x,y,z,r] which can be visualised on Paraview using sphere Glyphs (optional)*

     * *visualise*: whether to plot the clump and mesh (optional)*.

     * *rMax_ratio*: limit the maximum sphere radius to a ratio of the radius of the largest inscribed circle (optional)*.

:return: mesh: Trimesh object containing all relevant parameters of input polyhedron. The user can get:
                mesh.vertices
                mesh.faces
                mesh.centroid
                mesh.volume
                mesh.moment_inertia
                mesh.principal_inertia_components
                mesh.principal_inertia_vectors

        clump: Clump object containing all relevant clump parameters
                clump.positions		:	M-by-3 matrix containing the position of each generated sphere.
                clump.radii			:	M-by-1 vector containing the radius of each generated sphere
                clump.minRadius		:	Minimum generated sphere (might differ from rmin)
                clump.maxRadius		:	Maximum generated sphere
                clump.numSpheres	:	Total number of spheres
"""

if isinstance(inputGeom, str):
    if inputGeom.endswith(".stl"):
        mesh = trimesh.load_mesh(inputGeom)  # this will be used for voxalization
        F = mesh.faces
        P = mesh.vertices

        # Calculate extreme coordinates & centroid of the AABB of the particle
        minX, minY, minZ = np.min(P[:, 0]), np.min(P[:, 1]), np.min(P[:, 2])
        maxX, maxY, maxZ = np.max(P[:, 0]), np.max(P[:, 1]), np.max(P[:, 2])
        aveX, aveY, aveZ = np.mean((minX, maxX)), np.mean((minY, maxY)), np.mean((minZ, maxZ))

        # Center the particle to the centroid of its AABB
        P[:, 0] -= aveX
        P[:, 1] -= aveY
        P[:, 2] -= aveZ

        # Voxalization
        min_AABB = np.min(
            (np.abs(maxX - minX), np.abs(maxY - minY), np.abs(maxZ - minZ)))  # find the shortest length of axes
        voxel_size = min_AABB / div  # determine the voxel size

        img_temp = mesh.voxelized(pitch=voxel_size, method="subdivide").fill()  # voxalize

        intersection = np.pad(np.array(img_temp.matrix), ((2, 2), (2, 2), (2, 2)), mode='constant')  # pad the array with 2 voxels

        # I skipped the part "Ensure the voxel size is the same in all 3 directions -> Might be an overkill, but still".
        # Maybe add it later - Utku

    elif inputGeom.endswith(".mat"):

        # Load the .mat file
        data = loadmat(inputGeom)
        vox_keys = list(data.keys())

        # MATLAB's load function sometimes loads additional meta information starting with '__'
        # So, let's filter out such keys
        vox_keys = [key for key in vox_keys if not key.startswith('__')]

        vox = data[vox_keys[0]]

        img = vox['img'][0][0]
        voxel_size = vox['voxel_size'][0][0]

        opt = 2  # see vol2mesh function in iso2mesh
        isovalues = []  # see vol2mesh function in iso2mesh

        # Convert voxel data to a mesh using trimesh
        vox_mesh = trimesh.voxel.Voxel(img).as_mesh()
        P = vox_mesh.vertices
        F = vox_mesh.faces

        P = P * voxel_size[0][0]

        # Create a structure equivalent in Python
        FV = {'vertices': P, 'faces': F}

        RBP, _ = RigidBodyParameters.RBP(FV)
    else:
        raise ValueError("Not recognised inputGeom format.")
elif isinstance(type(inputGeom), type):  # check if it is a class ("struct").
    try:
        P = inputGeom.Vertices
    except AttributeError:
        P = inputGeom.vertices

    try:
        F = inputGeom.Faces
    except AttributeError:
        F = inputGeom.faces

    RBP, _ = RigidBodyParameters.RBP(F, P)

else:
    raise ValueError("Not recognised inputGeom format.")

rMax_ratio = kwargs.get('rMax_ratio')
if rMax_ratio is None or rMax_ratio<=0:
	rMax_ratio=1

################################################################################################
#                                   Main Body of the Function                                  #
################################################################################################

clump = Clump()  # instantiate Clump object for later use

# Dimensions of the new image
halfSize = np.array(intersection.shape) / 2


edtImage = distance_transform_edt(intersection)
rInscribed = np.max(edtImage)

radius = rInscribed * rMax_ratio
if radius * voxel_size < rMin:
    raise ValueError("Chosen radius smaller than rMin")

radius_phys = rMin
radius_vox = radius_phys / voxel_size


valid_centers = np.argwhere(edtImage >= radius_vox)
centroid_vox = np.round(halfSize).astype(int)


dists = np.linalg.norm(valid_centers - centroid_vox, axis=1)
i0 = np.argmin(dists)

centers = [valid_centers[i0]]


valid_set = set(map(tuple, valid_centers))
used_set  = set()

centers = [valid_centers[i0]]
used_set.add(tuple(centers[0]))

radius_factor = 0.9  # reduce by 10% each time we fail
radius_phys_current = radius_phys  # start with rMin
radius_vox_current = radius_phys_current / voxel_size


k = 1
# Precompute shell offsets
tol = 1.0
R = int(np.ceil(2*radius_vox_current + tol))
shell_offsets = np.array([[i,j,l] for i in range(-R,R+1) 
                                    for j in range(-R,R+1)
                                    for l in range(-R,R+1)
                                    if abs(np.linalg.norm([i,j,l]) - 2*radius_vox_current) <= tol])
from scipy.spatial import cKDTree
while k < N:
    candidates = np.unique(np.vstack([c + shell_offsets for c in centers]), axis=0)
    # Filter valid centers
    mask_valid = np.array([tuple(c) in valid_set and tuple(c) not in used_set for c in candidates])
    candidates = candidates[mask_valid]

    if len(candidates) == 0:
        radius_phys_current *= radius_factor
        if radius_phys_current < rMin * 0.5:
            print(f"Cannot place more spheres. Stopped at {k} spheres")
            break
        radius_vox_current = radius_phys_current / voxel_size
        # Update shell for new radius
        R = int(np.ceil(2*radius_vox_current + tol))
        shell_offsets = np.array([[i,j,l] for i in range(-R,R+1) 
                                            for j in range(-R,R+1)
                                            for l in range(-R,R+1)
                                            if abs(np.linalg.norm([i,j,l]) - 2*radius_vox_current) <= tol])
        continue

    # Use KDTree to reject overlapping spheres
    tree = cKDTree(np.array(centers))
    distances, _ = tree.query(candidates, k=1)
    candidates = candidates[distances >= 2*radius_vox_current - 1e-6]

    if len(candidates) == 0:
        radius_phys_current *= radius_factor
        if radius_phys_current < rMin * 0.5:
            print(f"Cannot place more spheres. Stopped at {k} spheres")
            break
        radius_vox_current = radius_phys_current / voxel_size
        # Update shell
        R = int(np.ceil(2*radius_vox_current + tol))
        shell_offsets = np.array([[i,j,l] for i in range(-R,R+1) 
                                            for j in range(-R,R+1)
                                            for l in range(-R,R+1)
                                            if abs(np.linalg.norm([i,j,l]) - 2*radius_vox_current) <= tol])
        continue

    # Pick candidate closest to centroid
    dists_to_centroid = np.linalg.norm(candidates - centroid_vox, axis=1)
    next_center = candidates[np.argmin(dists_to_centroid)]

    centers.append(next_center)
    used_set.add(tuple(next_center))
    k += 1


for c in centers:
    xyzC = c - halfSize + 1
    clump.positions = np.vstack((clump.positions, xyzC * voxel_size))
    clump.radii = np.vstack((clump.radii, radius_phys))


clump.minRadius = np.min(clump.radii)
clump.maxRadius = np.max(clump.radii)
clump.numSpheres = len(clump.radii)

output = kwargs.get('output')
if output is not None:
    np.savetxt(output, np.asarray(np.hstack((clump.positions, clump.radii))),
               delimiter=",")  # In PyCharm this line seems to have an error but it does not. Known issue.

outputVTK = kwargs.get('outputVTK')
if outputVTK is not None:
    clump_to_VTK(clump, filepath=outputVTK)

visualise = kwargs.get('visualise')
if visualise is not None and visualise:
    clump_plotter_pyvista(clump)

return mesh, clump

`