<jupyter_start><jupyter_text>Experiment with a single obj file<jupyter_code># Load a mesh and compute its distance field
m = gel.obj_load("../data/cochleas/shape20_pca.obj")
# View the object
viewer = gel.GLManifoldViewer()
viewer.display(m, mode='n', bg_col=[1, 1, 1])
del viewer
# Compute distance field
mdist = gel.MeshDistance(m)
s = sum([m.area(face) for face in m.faces()])
# Discretize the distance field
resolution = 64
boundary = -10
width = 20
step_size = width/resolution

output = []

for i in range(1,resolution+1):
    for j in range(1,resolution+1):
        for k in range(1,resolution+1):
            is_inside = mdist.ray_inside_test([boundary+step_size*i, boundary+step_size*j, boundary+step_size*k])
            if is_inside:
                output.append([i,j,k])

                
output = np.array(output)
np.savetxt(file+".txt", output, delimiter=",")<jupyter_output><empty_output><jupyter_text>Compute discrete distance field text files for range of cochleas.<jupyter_code># Load meshes.
cochlea_indices = range(len(cochlea_files))
cochlea_meshes = []
for cochlea_index in cochlea_indices:
    print("("+str(cochlea_index+1)+"/"+str(len(cochlea_indices))+") "+"Loading mesh from \""+cochlea_files[cochlea_index]+".obj\"...")
    cochlea_meshes.append(gel.obj_load(cochlea_folder+cochlea_files[cochlea_index]+".obj"))
print("--------------------------------------\nFinished loading meshes from "+str(len(cochlea_indices))+" files.")
# Compute distance fields.
cochlea_distance_fields = []
for cochlea_index in range(len(cochlea_meshes)):
    print("("+str(cochlea_index+1)+"/"+str(len(cochlea_meshes))+") "+"Generating distance field for \""+cochlea_files[cochlea_index]+".obj\"...")
    cochlea_distance_fields.append(gel.MeshDistance(cochlea_meshes[cochlea_index]))
print("--------------------------------------\nFinished computing distance fields for "+str(len(cochlea_meshes))+" meshes.")
# Voxelize distance fields.
resolution = 64
boundary = [-10, -10, -10]
width = 20
step_size = width/resolution
for cochlea_index in range(len(cochlea_distance_fields)):
    output = []
    print("("+str(cochlea_index+1)+"/"+str(len(cochlea_distance_fields))+") "+"Voxelizing "+cochlea_files[cochlea_index]+" to "+str(resolution)+"x"+str(resolution)+"x"+str(resolution)+" grid...")
    for i in range(1,resolution+1):
        if (i % 10 == 0):
            print("  "+str(round(i/resolution*100))+"% done...")
        for j in range(1,resolution+1):
            for k in range(1,resolution+1):
                is_inside = cochlea_distance_fields[cochlea_index].ray_inside_test([-boundary[0]-step_size*i, -boundary[1]-step_size*j, -boundary[2]-step_size*k])
                if is_inside:
                    output.append([i,j,k])
    output = np.array(output)
    output_file_name = cochlea_files[cochlea_index]+"_"+str(resolution)+".txt";
    np.savetxt(cochlea_folder+output_file_name, output, delimiter=",")
    print("  100% done. Output saved to \""+output_file_name+"\".")
print("--------------------------------------\nFinished voxelizing "+str(len(cochlea_distance_fields))+" distance fields.")<jupyter_output><empty_output><jupyter_text>Big all-in-one computation. Very heavy...<jupyter_code># Compute discrete distance field txt files for range of cochleas

mirror = False
resolution = 64
boundary = -12
width = 24
step_size = width/resolution
cochlea_indices = range(len(cochlea_files))

for file_index in cochlea_indices:
    print("("+str(file_index+1)+"/"+str(len(cochlea_indices))+") "+"Loading mesh from \""+cochlea_files[file_index]+".obj\"...")
    m = gel.obj_load(cochlea_folder+cochlea_files[file_index]+".obj")
    print("Creating distance field...")
    mdist = gel.MeshDistance(m)
    output = []
    print("Iterating over "+str(resolution)+"x"+str(resolution)+"x"+str(resolution)+" grid:")
    for i in range(1,resolution+1):
        if (i % 10 == 0):
            print("  "+str(i/resolution*100)+"% done...")
        for j in range(1,resolution+1):
            for k in range(1,resolution+1):
                if mirror:
                    is_inside = mdist.ray_inside_test([boundary+step_size*i, boundary+step_size*j, boundary+step_size*k])
                else:
                    is_inside = mdist.ray_inside_test([-boundary-step_size*i, -boundary-step_size*j, -boundary-step_size*k])
                if is_inside:
                    output.append([i,j,k])

    output = np.array(output)
    output_file_name = cochlea_folder+cochlea_files[file_index]+"_"+str(resolution)
    if mirror:
        output_file_name += "_mirrored"
    output_file_name += ".txt"
    np.savetxt(output_file_name, output, delimiter=",")
    print("  100% done.")
    print("Voxelized mesh saved to "+output_file_name+".")
    print("-------------------------------------------------------------------")
    del m, mdist, output
print("Finishes voxelizing "+str(len(cochlea_indices))+" meshes.")<jupyter_output><empty_output><jupyter_text>Special preprocessing based on the distance fieldBelow the code generates a voxel field, but instead of setting each entry to 1 it will be set to$$P(x) = \exp(-d_\alpha(x)) + \beta,$$where $\alpha$ is the clamp value and $\beta$ is the smallest value before normalisation. The preprocessing is done in matlab, here I merely save the signed distance for each voxel.<jupyter_code>def clamp(x,alpha):
    alpha = abs(alpha)
    if x > alpha:
        return alpha
    elif x < -alpha:
        return -alpha
    return x
resolution = 64
boundary = -12
width = 24
step_size = width/resolution
cochlea_indices = range(len(cochlea_files))
# cochlea_indices = range(0,1)

sign = 1;

for file_index in cochlea_indices:
    print("("+str(file_index+1)+"/"+str(len(cochlea_indices))+") ----------")
    print("  Loading mesh from \""+cochlea_files[file_index]+".obj\"...")
    m = gel.obj_load(cochlea_folder+cochlea_files[file_index]+".obj")
    print("  Computing distance field...")
    mdist = gel.MeshDistance(m)
    output = []
    print("  Iterating over "+str(resolution)+"x"+str(resolution)+"x"+str(resolution)+" grid:")
    for i in range(1,resolution+1):
        if (i % 10 == 0):
            print("    "+str(round(i/resolution*100))+"% done...")
        for j in range(1,resolution+1):
            for k in range(1,resolution+1):
                p = [-boundary-step_size*i, -boundary-step_size*j, -boundary-step_size*k]
                d = mdist.signed_distance(p)
                if (i == 1 and j == 1 and k == 1):
                    if d < 0:
                        sign = -1;
                    print("    (sign of distance = "+str(sign)+")")
                output.append([i,j,k,sign*d])
    
    output = np.array(output)
    output_file_name = cochlea_folder+cochlea_files[file_index]+"_"+str(resolution)+"_signed_distance"
    output_file_name += ".txt"
    np.savetxt(output_file_name, output, delimiter=",")
    print("    100% done.")
    print("  Voxelized distance field saved to "+output_file_name+".")
    del m, mdist, output
    sign = 1
print("Finished voxelizing "+str(len(cochlea_indices))+" meshes.")<jupyter_output><empty_output><jupyter_text>Computation of volume and surfaceVolume is computed using Halton-sequenced Monte Carlo sampling. Surface is computed by summing over the surfaces of all triangles in the mesh.<jupyter_code># n-th Corput sequence sampling function
def corput(n, base):
    q = 0
    bk = 1.0/base;
    while (n > 0):
        q  += (n % base) * bk;
        n  /= base;
        bk /= base;
    return q

# n-th Halton sequence sampling function (multidimensional Corput sampling)
def halton(n, d, shift):
    primes = (2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71) # supports up to 20 dimensions
    x = [0 for i in range(d)]
    assert(d < 20)
    for i in range(0,d):
        x[i] = corput(n,primes[i+shift])
    return x
# Halton-sequenced Monte Carlo integration given a pygel distance field mdist
from functools import reduce
def haltonmc(mdist,N,bounds,shift):
    # Compute volume here
    S = 0;
    V = reduce(lambda x,y: x*y, [bound[1]-bound[0] for bound in bounds]) # compute entire volue of the region
    
    r = np.array(list(map(lambda x: x[1]-x[0], bounds)))
    lower = np.array(list(map(lambda x: x[0], bounds)))
    upper = np.array(list(map(lambda x: x[1], bounds)))
    
    S = sum([1 if mdist.ray_inside_test(lower+r*halton(n, 3, shift)) else 0 for n in range(monte_carlo_points)])
    
    return S * V/N;

cochlea_indices = range(len(cochlea_files))

# Settings for Monte Carlo volume computation
compute_volume = True
monte_carlo_points = 50000;
halton_shift1 = 0;
halton_shift2 = 1;
bounds = [[-12, 12], [-12, 12], [-12, 12]]

# Settings for surface computation
compute_surface = False

print("Computing volume and/or surface area of " + str(len(cochlea_files)) + " meshes...")

V = []
S = []

for file_index in cochlea_indices:
    print("("+str(file_index+1)+"/"+str(len(cochlea_indices))+")")
    print("  Loading mesh from \""+cochlea_files[file_index]+".obj\"...")
    m = gel.obj_load(cochlea_folder+cochlea_files[file_index]+".obj")
    # Compute surface here
    if (compute_surface):
        print("  Integrating mesh surface...")
        s = sum([m.area(face) for face in m.faces()])
        S.append(s)
        print("    Surface = "+str(s))
    
    # Compute volume here
    if (compute_volume):
        print("  Computing volume...")
        print("    Computing distance field...")
        mdist = gel.MeshDistance(m)
        print("    Computing Monte Carlo volume integral...")
        v1 = haltonmc(mdist,monte_carlo_points,bounds,0)
        print("    Volume1 = "+str(v1))
        v2 = haltonmc(mdist,monte_carlo_points,bounds,1)
        print("    Volume2 = "+str(v2))
        v3 = haltonmc(mdist,monte_carlo_points,bounds,2)
        print("    Volume3 = "+str(v3))
        v4 = haltonmc(mdist,monte_carlo_points,bounds,3)
        print("    Volume4 = "+str(v4))
        v5 = haltonmc(mdist,monte_carlo_points,bounds,4)
        print("    Volume5 = "+str(v5))
        v6 = haltonmc(mdist,monte_carlo_points,bounds,5)
        print("    Volume6 = "+str(v6))
        v7 = haltonmc(mdist,monte_carlo_points,bounds,6)
        print("    Volume7 = "+str(v7))
        v8 = haltonmc(mdist,monte_carlo_points,bounds,7)
        print("    Volume8 = "+str(v8))
        v9 = haltonmc(mdist,monte_carlo_points,bounds,8)
        print("    Volume9 = "+str(v9))
        v10 = haltonmc(mdist,monte_carlo_points,bounds,9)
        print("    Volume10 = "+str(v10))
        V.append([v1, v2, v3, v4, v5, v6, v7, v8, v9, v10])
        del mdist
    
    del m
    
# Save to file
if (compute_volume):
    output_file_name = cochlea_folder+'cochleae_volumes.txt'
    np.savetxt(output_file_name, np.array(V), delimiter=",")
if (compute_surface):
    output_file_name = cochlea_folder+'cochleae_surfaces.txt'
    np.savetxt(output_file_name, np.array(S), delimiter=",")
print("Finished volume and/or surface computation for "+str(len(cochlea_indices))+" meshes.")<jupyter_output><empty_output><jupyter_text>Treating toy objects<jupyter_code>object_folder = "../data/objects/"
object_files = [
    "bunny",
    "duck"
]

object_files
def voxelise_mesh(resolution, object_files, object_folder):

    object_indices = range(len(object_files))
    sign = 1;
    for file_index in object_indices:
        print("("+str(file_index+1)+"/"+str(len(object_indices))+") ----------")
        print("  Loading mesh from \""+object_files[file_index]+".obj\"...")
        m = gel.obj_load(object_folder+object_files[file_index]+".obj")
        bbox = gel.bbox(m)
        print(bbox)
        width = round(1.2*(max(bbox[1])-min(bbox[0])))
        boundary = -round(1.2*abs(min(bbox[0])))
        step_size = width/resolution
        print(width, boundary, resolution, step_size)
        print("  Computing distance field...")
        mdist = gel.MeshDistance(m)
        output = []
        print("  Iterating over "+str(resolution)+"x"+str(resolution)+"x"+str(resolution)+" grid:")
        for i in range(1,resolution+1):
            if (i % 10 == 0):
                print("    "+str(round(i/resolution*100))+"% done...")
            for j in range(1,resolution+1):
                for k in range(1,resolution+1):
                    p = [boundary+step_size*i, boundary+step_size*j, boundary+step_size*k]
                    d = mdist.signed_distance(p)
                    if (i == 1 and j == 1 and k == 1):
                        if d < 0:
                            sign = -1;
                        print("    (sign of distance = "+str(sign)+")")
                    output.append([i,j,k,sign*d])
                        

        output = np.array(output)
        output_file_name = object_folder+object_files[file_index]+"_"+str(resolution)+"_signed_distance"
        output_file_name += ".txt"
        np.savetxt(output_file_name, output, delimiter=",")
        print("    100% done.")
        print("  Voxelized distance field saved to "+output_file_name+".")
        del m, mdist, output
        sign = 1
    print("Finished voxelizing "+str(len(object_indices))+" meshes.")


resolution = 64

voxelise_mesh(resolution, object_files, object_folder)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1609438814097716, 'fraction_numerical': 0.02011798517622145, 'mean_word_length': 3.090009279307145, 'pattern_counts': {'":': 0, '<': 29, '<?xml version=': 0, '>': 27, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '407102', 'license': 'MIT', 'n_tokens_mistral': 4788, 'n_tokens_neox': 4424, 'n_words': 951, 'path': 'python/GEL Distance Field Generation from OBJ Mesh.ipynb', 'repo_name': 'redzin/OptimTrans', 'size': '13222'}