<jupyter_start><jupyter_text>Make a (slightly corrupt) spherical topographic map of Mercury and my cat!<jupyter_code>import numpy as np
import math
from PIL import Image
Image.MAX_IMAGE_PIXELS = None # working with very high def images, so override this lol
def scale_data(img_path = mercury_demo, min_thresh = 0.0, max_thresh = None, protrusion = 5.0):
    ''' 
    A function that loads the image flux data and scales the flux values to fall between 0 and a desired protrusion value.
    
    Parameters
    ----------
    img_path: String, optional
        image file location (default mercury_demo, a very high resolution topographic map of Mercury)
    min_thresh: float, optional
        minimum desired flux threshhold value (default 0.0)
    max_thresh: float, optional
        maximum desired flux threshold value (default None)
    protrusion: float, optional
        max height value to scale image fluxes to
        
    Returns
    -------
    img2D: ndarray
        2D array of image with scaled flux values
    '''
    image = Image.open(img_path)
    img1D = np.array(image.getdata())
    
    if max_thresh == None:
        max_height = np.max(img1D)
    else:
        max_height = max_thresh
        img1D[img1D > max_height] = max_height # set all values abolve max to max
        
    if min_thresh != 0: 
        min_height = min_thresh
    else:
        min_height = np.min(img1D)
        img1D = img1D - min_height # overwrite, subtract floor from all values
        
    if (max_height == min_height): # try to avoid dividing by zero
        height_scalar = 1
    else:
        height_scalar = protrusion / abs(max_height - min_height)
    
    img1D = img1D * height_scalar # now on a scale from 0 to 5 mm
    
    img2D = np.reshape(img1D, (image.size[1], image.size[0])) # finally shape the array into an image

    return img2D
def define_shape(resolution = 100, radius = 100.0): 
    ''' 
    A function that builds the spherical sampling mesh. It begins with an icosahedron and divides the vertices of each face into a number of smaller triangles specivied by the resolution input. 
    
    Parameters
    ----------
    resolution: int, optional
        how many facets along one edge of the original icosahedron (default 100 for a balance of runtime vs resolution)
    radius: float, optional
        radius of base sphere in mm before the image data is overlayed (default 100 mm for now) 
        
    Returns
    -------
    icos_mesh: ndarray
        [5 x 2n x n x 3] array simulating the shape of the net of the icosahedron filled in with its populated faces. At each index in [5 x 2n x n] is the 1D array of spherical coordinates of each point [r, 0, φ]
        the bulk of the heavily populated icosahedron, including everything other than the very top and bottom. 
    icos_apex: array
        1D array of sphere's apex [r, 0, φ] (Note: I realize I could have only passed the radius and simplified calculations down the line by setting z = radius, but that will be for an update)
    icos_base: array
        1D array of sphere's base [r, 0, φ] (Note: same as apex)
    '''    
    
    # DEFINE SPHERICAL COORDS OF ICOSAHEDRON VERTICES
    # sphere coords are very easy to deal with at this point because I can just override the radius with height later
    # ORDER: (r, 0, φ) defintions: 0 is right-hand sweep from x to y, φ is angle down from upwards z
    # (I made these coordinates - three equally spaced phi values, two sets of five equally spaced staggered theta values)
    v0  = np.array([radius,   0, 180]) # base
    v1  = np.array([radius,   0, 120])
    v2  = np.array([radius, 324,  60])
    v3  = np.array([radius,  72, 120])
    v4  = np.array([radius,  36,  60])
    v5  = np.array([radius, 144, 120])
    v6  = np.array([radius, 108,  60])
    v7  = np.array([radius, 216, 120])
    v8  = np.array([radius, 180,  60])
    v9  = np.array([radius, 288, 120])
    v10 = np.array([radius, 252,  60])
    v11 = np.array([radius,   0,   0]) # apex
        
    icos_apex = v11 # top of icos
    icos_base = v0 # bottom of icos
    
    # 10 horizontal pairs, bottom to top left to right: (order read off from a grid-like net of an icosahedron)
    band = [[v1, v3], [v2, v4], [v3, v5], [v4, v6], [v5, v7], [v6, v8], [v7, v9], [v8, v10], [v9, v1], [v10, v2]] 
    # hot dang I think I can populate the sphere without even using the top or botom of the icos in calculations (?)
    
    n = resolution # sub triangle density per icos vertex, n is a good shorthand variable of number density  
    
    # all the steps will be the same, so calculate it once: 
    theta_step = (v3[1] - v1[1]) / n # calculate one step in theta here (distance / #steps)
    phi_step = 180 / (3*n) # woohoo symmetry
    
    big_squares = np.zeros(10, dtype = np.ndarray) # store (n x n) local population squares in here, concatenate later
    
    for i in range(len(band)): # there are 10 sub bands (horizontal vectors) in the list that defines the icos middle band  
        #   A . B
        #    .   .   <== my parallelogram :)
        #     C . D
    
        A = band[i][0] # left point
        B = band[i][1] # right point  
        
        pop = np.zeros((n, n), dtype = np.ndarray) # make a temp 2D grid to be populated, local population
          
        pop[0][0] = A # set grandparent/base node to be first entry
        
        # calculate the diagonal entries:
        for ii in range(n): # recall order (r, 0, φ)
            theta = (A[1] + theta_step * ii) % 360 # starting (left) angle (A's 0) with n equal spacings to the next point's 0.
            phi = (A[2]) % 180   # both A and B have the same altitude/same φ, so it doesn't matter which one I choose
            pop[ii][ii] = [1, theta, phi] # fill in the diagonal with the points equally spaced on the horizontal line
            
        for ii in range(n): 
            parent = pop[ii][ii] # the rest of the entries in the mesh will be calculated from the horizontal parent nodes
            
            # children above: parent index ii has ii children above it. 
            # node 0 is the grandparent node. it is on a vertex of the icos and is at the top left of matrix. 
            for jj in range (0, ii):
                
                # THIS THETA AND THE THETA BELOW ARE WHAT I SUSPECT TO BE THE CAUSE OF THE BUG THAT MAKES BIG PETALS ON MY SPHERE, BOTH ON THE TOP AND BOTTOM
                # I have spent so long trying to fix it. it's some sort of issue of me not scaling their spacings depending on their height in the face. 
                # I have more ideas that I'll try eventually but not now unfortunately :(
                theta = (parent[1] - abs(jj - ii) * 0.5 * theta_step) % 360 # distance from diag corresponds to half step
                phi = (parent[2] - phi_step * (ii - jj)) % 180 # distance from diagonal corresponds to one phi step
                pop[jj][ii] = [1, theta, phi] # first index used to be ii - jj but idt that's right
                
            # children below: parent index ii has (n - 1) - ii children below it. 
            # the sum of parent's index plus number of children is (n - 1) !
            for jj in range(ii + 1, n): # here we skip the last parent because it has no lower kids, being at (n,n)
                
                # BUG LOCATION TO FIX LATER
                theta = (parent[1] + abs(jj - ii) * 0.5 * theta_step) % 360 
                phi = (parent[2] + phi_step * (jj - ii)) % 180
                pop[jj][ii] = [1, theta, phi] 
             
        big_squares[i] = pop # fill big squares temp array with our ten local pop arrays                
       
    icos_mesh = np.zeros(5, dtype = np.ndarray) # store the five big columns of the grid array, each column made of two big squares
    
    for i in range (5): # fill in the big grid array in chunks
        top = big_squares[2 * i + 1] # odd indices in big_squares correspond to top row of big grid array
        bottom = big_squares[2 * i] # even indices in big_squares correspond to bottom row of grid
        icos_mesh[i] = np.concatenate((top, bottom), axis = 0, dtype = np.ndarray) # build 2n tall column
    
    return icos_mesh, icos_apex, icos_base # the top and bottom of the icos, and the bulk of the shape
    # I found it interesting that the total number of points in the populated icosahedron is 10n^2 + 2 (with n being the resolution)
            
def data_overlay(sampling_data, icos_mesh, icos_apex, icos_base): 
    ''' 
    A function that updates the sphere mesh by sampling the image data at the cartesian-converted coordinates of the sphere and using their scaled fluxes to override the radii of the sphere mesh coordinates. 
    
    Parameters
    ----------
    sampling_data: ndarray
        2D array of image with scaled flux values to sample from
    icos_mesh: ndarray
        [5 x 2n x n x 3] array simulating the shape of the net of the icosahedron filled in with mesh points of the populated faces
    icos_apex: array
        1D array of the sphere's apex [r, 0, φ]
    icos_base: array
        1D array of sphere's base [r, 0, φ]
        
    Returns
    -------
    mesh_data: ndarray
        icos_mesh input with updated radii sampled from the sampling data
    mesh_apex: array
        icos_apex updated with radius sampled from the sampling data
    mesh_base: array
        icos_base updated with radius sampled from the sampling data
    '''
    
    x_pixels = len(sampling_data[0])
    y_pixels = len(sampling_data)
    print("sampling data dimensions:  x:", x_pixels, " y:", y_pixels)
    print("For the best quality 3D model, consider choosing a mesh resolution of about 10-20% the x-resolution of your image.")

    x_scale = x_pixels / 360 # for theta
    y_scale = y_pixels / 180 # for phi
    
    sphere_radius = icos_apex[0] # all radii should be the same right now
    
    mesh_data = np.zeros((len(icos_mesh), len(icos_mesh[0]), len(icos_mesh[0][0])), dtype = np.ndarray) 
    
    # go vertex by vertex of triangle, ask the dataset what height it should be at, and push vertex to that height        
    for ii in range (len(icos_mesh)): # which of the 5 big columns
        for jj in range (len((icos_mesh[0]))): # which row of a big column
            for kk in range (len((icos_mesh[0][0]))): # which entry (real column)
                                
                theta = icos_mesh[ii][jj][kk][1] 
                phi = icos_mesh[ii][jj][kk][2]
                
                x = int(math.floor(theta * x_scale) - 1) # turn theta and phi into an index to get a pixel from
                y = int(math.floor(phi * y_scale) - 1)
                
                # AUDREY in the future you can just add the sampling radii to the point bc all points are already at sphere radius.
                radius = sampling_data[y][x] + sphere_radius # we calculated this earlier, add to base radius of sphere
                
                sphere_coord = [radius, theta, phi]
                cart_coord = spherical_to_cartesian(sphere_coord)
                
                # fill in pixel grid
                mesh_data[ii][jj][kk] = cart_coord
    
    # AUDREY FIX THESE AFTER EVERYTHING RUNS
    # now do the apex and the base    
    apex_x = int(math.floor(icos_apex[1] * x_scale) - 1) # turn it into an index to get a pixel from
    apex_y = int(math.floor(icos_apex[2] * y_scale) - 1)
                
    apex_radius = sampling_data[apex_y][apex_x] + sphere_radius # we calculated this earlier
                
    apex_sphere = [apex_radius, icos_apex[1], icos_apex[2]]
    mesh_apex = spherical_to_cartesian(apex_sphere)
    
    base_x = int(math.floor(icos_base[1] * x_scale) - 1) # turn it into an index to get a pixel from
    base_y = int(math.floor(icos_base[2] * y_scale) - 1)
                
    base_radius = sampling_data[base_y][base_x] + sphere_radius # we calculated this earlier
                
    base_sphere = [base_radius, icos_base[1], icos_base[2]]
    mesh_base = spherical_to_cartesian(base_sphere)
    
    return mesh_data, mesh_apex, mesh_base                                
# spherical to cartesian coordinate converter

def spherical_to_cartesian(sphere_coord): 
    ''' 
    A function that converts a spherical coordinate to a 3D cartesian coordinate with trig.
    
    Parameters
    ----------
    sphere_coord: array
        [r, 0, φ] array of spherical coordinate to be converted 
        
    Returns
    -------
    cart_coord: array
        [x, y, z] array of cartesian coordinate
    '''
    radius = sphere_coord[0]
    theta = sphere_coord[1] * math.pi / 180
    phi = sphere_coord[2] * math.pi / 180
    
    x = radius * math.sin(phi) * math.cos(theta)
    y = radius * math.sin(phi) * math.sin(theta)
    z = radius * math.cos(phi)
    
    cart_coord = [x, y, z]
    
    return cart_coord
def stitch_triangles(A, B, C, D, lines): 
    ''' 
    A void function that defines two right-handed triangular facets for a parallelogram with points A, B, C, and D, and writes the corresponding lines of the .STL file for each triangle. The lines of the .STL file are then appended to a specified list. 
    
    Parameters
    ----------
    A: array
        3D cartesian coordinate defining the upper-left point of a parallelogram
    B: array
        3D cartesian coordinate defining the upper-right point of a parallelogram
    C: array
        3D cartesian coordinate defining the lower-left point of a parallelogram
    D: array
        3D cartesian coordinate defining the lower-right point of a parallelogram
    lines: list
        list to append the lines of the .STL file for each triangle to
        
    Returns
    -------
    None
    '''
    #   A . B
    #    .   .   <== my parallelogram :)
    #     C . D
    
    # keywords for writing an STL file, with right amount of spacing
    facet = '  facet' # two spaces before
    outerloop = '    outer loop' # four spaces before, two words
    vertex = '      vertex ' # six spaces before, one after
    endloop = '    endloop' # four spaces before, one word
    endfacet = '  endfacet' # two spaces before, one word
    
    lines.extend((facet, outerloop, vertex + str(A)[1:-1], vertex + str(C)[1:-1], vertex + str(D)[1:-1], endloop, endfacet))
    lines.extend((facet, outerloop, vertex + str(A)[1:-1], vertex + str(D)[1:-1], vertex + str(B)[1:-1], endloop, endfacet))
def make_file(lines, filename): # lines of the file
    ''' 
    A void function that writes the .STL file from the individually written list of lines, and saves it with a specified file name. 
    
    Parameters
    ----------
    lines: list
        list of lines to make up the file
    filename: String
        name of the .STL file (passed from constructor with default '3D_map')
        
    Returns
    -------
    None
    '''
    
    with open(filename + ".stl", 'w') as f:
        f.write('\n'.join(lines))
def constructor(mesh, apex, base, filename = '3D_map'):
    ''' 
    A void function that sends paralellograms made of coincident points to the triangle stitcher to result in a closed solid defined by 20*n^2 facets and writes a .STL file..
    
    Parameters
    ----------
    mesh: ndarray
        spherical mesh with radii sampled from the sampling data
    apex: array
        apex of the spherical mesh with radius sampled from the sampling data
    base: array
        base of the spherical mesh with radius sampled from the sampling data
    filename: String, optional
        name of the .STL file (default '3D_map')
        
    Returns
    -------
    None
    '''
    n = len(mesh[0][0]) # resolution of the mesh
    
    lines = [] # initiate list to hold lines of the .STL file
    lines.append('solid') # first line of the .STL file is "solid", initiating the definition of a solid
    
    for ii in range (len(mesh)): # 5 big columns
        for jj in range (0, 2*n - 1): # jj from 0 to 2*n - 2
            for kk in range (0, n - 1): # kk from 0 to n - 2
            
                # GENERAL CASE (the three other points are to the right and below)
                A = mesh[ii][jj][kk]
                B = mesh[ii][jj][kk + 1]
                C = mesh[ii][jj + 1][kk]
                D = mesh[ii][jj + 1][kk + 1]
                
                stitch_triangles(A, B, C, D, lines) # stitch triangles every time we name a NEW set of points A, B, C, D
        
            kk = n - 1 # manually set kk for the weird right column cases
            A = mesh[ii][jj][kk] # update who A is every time we name a different set of indices
            
            # CASE: OPEN TOP CORNER (I think these uppercase notes only make sense if I draw you a diagram of how I'm storing them)
            if (jj < n): # last col of segment, and in the part that sticks above the next segment
                
                # SUBCASE: VERTEX CASE
                if (jj == 0): # every segment's top right most entry must connect to the icos vertex 
                    B = apex # AUDREY swap the apex and base if the file looks werid !!
                
                # SUBCASE: "GENERAL" OPEN TOP CORNER CASE
                else: 
                    B = mesh[(ii + 1) % 5][0][n - jj] 
                
                # last two points calculated the same way for all OPEN CORNER
                C = mesh[ii][jj + 1][n - 1]
                D = mesh[(ii + 1) % 5][0][n - 1 - jj]
            
            # CASE: SEGMENT SPANNING
            else: # (means jj >= n)
                B = mesh[(ii + 1) % 5][jj - n][0]
                C = mesh[ii][jj + 1][n - 1]
                D = mesh[(ii + 1) % 5][jj - n + 1][0]
        
            stitch_triangles(A, B, C, D, lines) # stitch triangles for new points
            
        # CASE: BOTTOM ROW
        jj = 2*n - 1 # manually set jj for the weird bottom row cases
        
        # SUBCASE: VERTEX CASE
        kk = 0 # manually set kk for vertex case
        
        A = mesh[ii][jj][kk] # new kk, new A
        B = mesh[ii][2*n - 1][1] # FUTURE AUDREY: make it work with n = 1 fool
        C = base # make one triangle that connects with the vertex
        D = mesh[(ii + 1) % 5][2*n - 1][0]
        
        stitch_triangles(A, B, C, D, lines)
        
        # SUBCASE: "GENERAL" OPEN BOTTOM CORNER CASE
        # FUTURE AUDREY NOTE: there's a bug here that makes the petal shaped holes in the bottom of the sphere. 
        # it just means that two of your points are defined incorrectly but you gotta figure that out. 
        # currently two points are in the right spot, and the other two points are sent to the top of the sphere. 
        for kk in range (1, n - 1): # will go from 1 to n - 2
            A = mesh[ii][jj][kk] # new kk, new A
            B = mesh[ii][jj][kk + 1]
            C = mesh[ii][n - kk][kk]
            D = mesh[ii][n - kk - 1][kk] 
            
            stitch_triangles(A, B, C, D, lines)
        
        # SUBCASE: BOTTOM RIGHT "WEIRD FLAT TRIANGLE" CASE
        kk = n - 1 # manually set kk for the bottom right column case
            
        A = mesh[ii][jj][kk] # new kk, new A    
        B = mesh[(ii + 1) % 5][n - 1][0] # requires three points (B, D, C) to be in a line
        C = mesh[(ii + 1) % 5][n][0]
        D = mesh[(ii + 1) % 5][n + 1][0]
        
        stitch_triangles(A, B, C, D, lines)
       
    lines.append('endsolid') # last line of the STL file, means that the solid is fully defined and closed by the listed facets.
    
    make_file(lines, filename) # writes the file and downloads it to the same folder as this python script<jupyter_output><empty_output><jupyter_text>- - - - - - - - - - - - - - - - - - - - - - RENDER HERE!- - - - - - - - - - - - - - - - - - - - - -<jupyter_code>mercury_demo = 'pics/mercury1.png' # awesome hi-res topographic pic of mercury
moon_demo = 'pics/moonbump4k.jpg'  # topographic moon
ringo_cat = 'pics/ringo.png' # this is my friend's fat cat named <NAME>
image = scale_data()
# define points of a mesh of a certain radius and resolution. this base shape can be kept and used for any image. 
# raising the resolution to something like 500 will create an awesome model but will take a few minutes to render. 
# (it will also be a massive file, about 1.2GB because it is 20*n^2 = 5 million triangles)
icos_mesh, icos_apex, icos_base = define_shape(resolution = 200)
# drape the image height data onto the mesh
mesh_data, mesh_apex, mesh_base = data_overlay(image, icos_mesh, icos_apex, icos_base)
# construct and write an STL file
constructor(mesh = mesh_data, apex = mesh_apex, base = mesh_base, filename = 'mercury0')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08895675780492535, 'fraction_numerical': 0.018276129532674035, 'mean_word_length': 2.119933454325469, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '238300', 'license': 'MIT', 'n_tokens_mistral': 6631, 'n_tokens_neox': 6104, 'n_words': 2580, 'path': 'star lamp.ipynb', 'repo_name': 'shadowbehindthebread/star_lamp', 'size': '20633'}