3D_V2.ipynb

Kernel: Python 2 (SageMath)

#%load_ext sage
import stlwrite as stl
import cPickle as pickle
import numpy as np
from __future__ import division
import matplotlib.delaunay as triang
import pylab

/projects/sage/sage-6.10/local/lib/python2.7/site-packages/matplotlib-1.5.0-py2.7-linux-x86_64.egg/matplotlib/cbook.py:136: MatplotlibDeprecationWarning: The matplotlib.delaunay module was deprecated in version 1.4. Use matplotlib.tri.Triangulation instead.
  warnings.warn(message, mplDeprecation, stacklevel=1)

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

def create_face(x_range, y_range, z_range):
    #print('xrange : ',x_range)
    #print('yrange : ',y_range)
    #print('zrange : ',z_range)
    test_ranges = [ x_range, y_range, z_range]
    ranges = []
    for i, r in enumerate(test_ranges):
        if type(r) != np.ndarray:
            constant_index = i
        else:
            ranges.append(r)
    
    assert len(ranges) == 2
    X,Y = np.meshgrid(ranges[0], ranges[1])
    X = X.flatten()
    Y = Y.flatten()
    
    cens,edg,tri,neig = triang.delaunay(X, Y)
    faces = []
    for t in tri:
        # t[0], t[1], t[2] are the points indexes of the triangle
        t_i = [index for index in t]
        #t_i.sort()
        
        x_vals = X[t_i]
        y_vals = Y[t_i]
        #pylab.plot(X[t_i],Y[t_i])
        
        face = []
        face = [ [ x_range, x_vals[index], y_vals[index] ]  for index in range(3) ] if constant_index == 0 else face
        face = [ [ x_vals[index], y_range, y_vals[index] ]  for index in range(3) ] if constant_index == 1 else face
        face = [ [ x_vals[index], y_vals[index], z_range ]  for index in range(3) ] if constant_index == 2 else face
        faces.append( face )
    #pylab.show()

    return faces

#create model x,y grid points (assumes it's symetric about (0,0) )
WIDTH = 1#the size of one voxel in the model
SCALE = 1
FACE_POINTS = 3#How many points go into each face

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_testV2_small.stl'
with open(filename, 'wb') as fp:
    writer = stl.ASCII_STL_Writer(fp)
    #writer = stl.Binary_STL_Writer(fp)
    
    left_x  = x_grid[0]
    right_x = x_grid[-1]
    top_y   = y_grid[0]
    bot_y   = y_grid[-1]
    
    #xy_pairs = [(left_x,bot_y), (right_x,bot_y), (right_x,top_y), (left_x,top_y)]
    #height_values = [0]*4
    #writer.add_faces( create_face( xy_pairs, height_values) )
    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]
            #print("left_x : {0}\n right_x : {1}\n top_y : {2}\n bot_y : {3}".format(left_x,right_x,top_y,bot_y))
            
            x_range = np.linspace(left_x, right_x, FACE_POINTS )
            y_range = np.linspace(top_y, bot_y, FACE_POINTS )
            faces += create_face(x_range, y_range, h*SCALE) 
                
            #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
                h_range = np.linspace(bot_edge*SCALE, h*SCALE, FACE_POINTS)
                faces += create_face(x_range, bot_y, h_range)

            
            #right side
            if x_index < (num_x-1):
                bot_edge = row[x_index+1]
            else:
                bot_edge = 0
            if bot_edge != h:
                h_range = np.linspace(bot_edge*SCALE, h*SCALE, FACE_POINTS)
                faces += create_face(right_x, y_range, h_range)
            
            #top side
            if y_index != 0:
                bot_edge = data[y_index-1][x_index]
            else:
                bot_edge = 0
            if bot_edge != h:
                h_range = np.linspace(bot_edge*SCALE, h*SCALE, FACE_POINTS)
                faces += create_face(x_range, top_y, h_range)
            
            
            #left side
            if x_index != 0:
                bot_edge = row[x_index-1]
            else:
                bot_edge = 0
            if bot_edge != h:
                h_range = np.linspace(bot_edge*SCALE, h*SCALE, FACE_POINTS)
                faces += create_face(left_x, y_range, h_range)
            
            writer.add_faces( faces )