3D_V1.ipynb

Kernel: Python 2 (SageMath)

#%load_ext sage
import stlwrite as stl
import cPickle as pickle
import numpy as np
from __future__ import division

with open("../data/10^2.pkl") as f:
    data = pickle.load(f)["iterations"]
    data += 1

def create_face(xy_pairs, height_values):
    #points should be defined clockwise as seen from the outside of the model
    face = []
    for index, value in enumerate(xy_pairs):
        face.append( value + (height_values[index],) )
    return face

#create model x,y grid points (assumes it's symetric about (0,0) )
WIDTH = 1#the size of one voxel (x,y) in the model
SCALE = 1#the height multiplier

num_y, num_x = data.shape

x_edge = (num_x//2)*WIDTH + WIDTH/2
y_edge = (num_y//2)*WIDTH + WIDTH/2

x_grid = np.linspace(-x_edge, x_edge, num_x+2)
y_grid = np.linspace(y_edge, -y_edge, num_y+2)

filename = '../3D_files/edges_testV1_small-1.stl'
with open(filename, 'wb') as fp:
    writer = stl.ASCII_STL_Writer(fp)
    
    for y_index, row in enumerate(data):
        for x_index, value in enumerate(row):
            faces = []
            #create lid at height h*scale
            h = value
            left_x  = x_grid[x_index]
            right_x = x_grid[x_index+1]
            top_y   = y_grid[y_index]
            bot_y   = y_grid[y_index+1]
            
            xy_pairs = [(left_x,bot_y), (left_x,top_y), (right_x,top_y), (right_x,bot_y)]
            height_values = [h*SCALE]*4
            faces.append( create_face( xy_pairs, height_values) )
            
            #add sides           
            #bottom side
            if y_index < (num_y - 1):
                #get bottom height from that neighbour
                bot_edge = data[y_index+1][x_index]                    
            else:
                #we are the last row
                #so bottom height is zero
                bot_edge = 0
            if bot_edge != h:#Don't add face if neighbour is same height
                xy_pairs      = [(left_x,bot_y),(left_x,bot_y),(right_x,bot_y),(right_x,bot_y)]
                height_values = [bot_edge*SCALE] + [h*SCALE]*2 +  [bot_edge*SCALE]
                faces.append( create_face( xy_pairs, height_values) )
                
            
            #right side
            if x_index < (num_x-1):
                bot_edge = row[x_index+1]
            else:
                bot_edge = 0
            if bot_edge != h:
                xy_pairs      = [(right_x,bot_y),(right_x,bot_y),(right_x,top_y),(right_x,top_y)]
                height_values =  [bot_edge*SCALE] + [h*SCALE]*2 +  [bot_edge*SCALE]
                faces.append( create_face( xy_pairs, height_values) )
            
            #top side
            if y_index != 0:
                bot_edge = data[y_index-1][x_index]
            else:
                bot_edge = 0
            if bot_edge != h:
                xy_pairs      = [(right_x,top_y),(right_x,top_y),(left_x,top_y),(left_x, top_y)]
                height_values =  [bot_edge*SCALE] + [h*SCALE]*2 +  [bot_edge*SCALE]
                faces.append( create_face( xy_pairs, height_values) )
            
            
            #left side
            if x_index != 0:
                bot_edge = row[x_index-1]
            else:
                bot_edge = 0
            if bot_edge != h:
                xy_pairs      = [(left_x,top_y), (left_x,top_y),(left_x,bot_y),(left_x,bot_y)]
                height_values = [bot_edge*SCALE] + [h*SCALE]*2 +  [bot_edge*SCALE]
                faces.append( create_face( xy_pairs, height_values) )
            
            writer.add_faces( faces )