In my blog post about cutting a turbine, I manually extracted the coordinates from the GCODE file, deleted extra coordinates, and moved rows around in Excel in order to synchronize the coordinates. This was laborious time consuming work, and it would be great if there was a program that automated this. This post is a work in progress hoping to develop an algorithm to do just that. Here it is on Github: https://github.com/agb253/3DHotwire
My previous turbine 3D model was quite large. Let’s only take the first four layers of the GCODE to help understand the coordinate synchronization algorithm.
Here is an algorithm to find the line going from a point on layer one to a point on layer 4 that minimizes the distance from the line to the surface of the model. It does this by finding the minimum mean square distance between between the line and coordinates on layer 2 and layer 3.
Here are the first four layers in our GCODE. The red points are examples of coordinates that should be synchronized.
The synchronization program will:
-
- Read in the GCODE file.
- Remove everything except the coordinates.
- Separate each layer into its own list. P1 and P2 are the top and bottom perimeters. The interior perimeters will be called iP1 through iPN from top to bottom
- Starting with the first location in the gcode for both P1[1] and P2[1] (in this example, let’s suppose we only have one part), find the distance of each point on iP1 to the line.
- Store the closest point to the line as C-iPN (closest interior perimeter point).
- Repeat step 4 and step 5 for every iPN (interior layers: iP1, iP2,…,iPN), storing the closest point on each iP to the line created in step 4.
- Calculate the sum of the distances of C-iP’s to the line. Store this number as the objective function.
- Remain on the same point on P1 or P1[1], but now draw a line to the next coordinate on P2 or P2[2]. (OBJ[1,1]= the objective function for P1[1] to P2[1], OBJ[1,2]=P1[1] to P2[2])
- Repeat step 4 through 8 for every coordinate in P2[1 through n].
- Find the lowest OBJ[1,j]; that is the P2[j] coordinate that should be synchronized with P1[1].
- For every point on P1[1 through N], repeat step 4-10.
You can download the GCODE file for the turbine blade at this link:
3DHotwire/Turbine_blade-all_GCODE.txt.gcode
Here is the synchronization code written in Python. I’m not a professional programmer, so excuse the mess. It works!! …kinda… and it is very slow since it is checking every line created by the combination of coordinates from P1[1-N] to P2[1-N], and finding the distance from every coordinate to the every line on every interior layer. To speed things up, I only used one out of every 20 interior layers.
import re from tkinter import Tk from tkinter.filedialog import askopenfilename import time Tk().withdraw() filename = askopenfilename()#prompt user to select the gcode file i=1 coords= [] coordinates = [] with open(filename) as gcode: for line in gcode: line = line.strip() coord = re.findall(r'[XY].?\d+.\d+', line) if coord: a = coord[0] b = coord[1] aa = float(a[1:]) bb = float(b[1:]) c = [aa,bb,i] coords.append(c) elif "layer_num" in line: coordinates.append(coords) coords = [] i=i+1 ####calculate the shortest distance between a point and line(defined by two points) import numpy as np P1 = coordinates[1] #first layer is not coordinates[0] as that is blank P2 = coordinates[len(coordinates)-1]#last layer j = 0 #top profile coordinate iterator k = 0 #bottom profile iterator l = 2 #interior layer iterator. Starts at 2 because first layer is empty [], and the first layer is P1 m = 0 #interior coordinate iterator min_dist = 99999999.9 #some large distance for each point in a layer sum_min_dist = 0 lowest_dist = 9999999.9#some large distance for total of layer's minimum distance start = time.time() def distance(p, q, rs): x = p-q return np.linalg.norm( np.outer(np.dot(rs-q, x)/np.dot(x, x), x)+q-rs, axis=1) for sets in P1: for points in P2: while l < len(coordinates)-1: #don't count the first layer as it is blank [] and don't use P1, start at coordinates[2] to last coordinate -1 while m < len(coordinates[l]): p = np.array(P1[j]) #coord of top profile q = np.array(P2[k]) #coord of bottom profile rs = np.array(coordinates[l][m])#interior layer at the m-th point d = distance(p, q, rs) if min_dist > d: #first time around min_dist should be large number, so this will be true min_dist = d #distance to the closest coordinate [m] on a layer [l] #print("min_dist " +str(min_dist)) m = m + 1 #next coordinate on the same layer [l] #end result is min_dist is equal to the distance to the closest coordinate sum_min_dist = sum_min_dist + min_dist #sum of the closest coordinate distances on each layer l = l +20 #next layer. You can skip layers by making l+ some number of skipped layers min_dist = 9999999.9 #initialize at some large distance m = 0 #which coordinate on P2 has the lowest sum_min_dist? #print(sum_min_dist) if sum_min_dist < lowest_dist: #if the sum of the closest distances for all layers is less than lowest_dist = sum_min_dist #the previous lowest sum of closest distance, then it may be a synch point synch_coord = P2[k] #so keep that new lowest distance and coordinate in memory sum_min_dist = 0 #reset the sum of the closest coordinates and try the next coordinate on P2 k = k + 1 #try the next coordinate on P2 l = 2 print(str(P1[j]) +" and " +str(synch_coord)) end = time.time() print(end - start) j = j + 1 #next coordinate on P1 lowest_dist = 9999999.9 #some large number k = 0
I was given the advice to use Anaconda3 for python compiler, write the code in a series of separate functions for easier parallelization, and try the “Numba” just in time compiler. It uses the “@jit” decorator. Whatever that means… Anyway, here is that code. Unfortunately, it runs slower than it did without Numba (68.15 sec per synchronized coordinate without Numba, and 84.65 sec with Numba).
import re from tkinter import Tk from tkinter.filedialog import askopenfilename import time import numpy as np from numba import jit Tk().withdraw() filename = askopenfilename()#prompt user to select the gcode file i=1 coords= [] coordinates = [] with open(filename) as gcode: for line in gcode: line = line.strip() coord = re.findall(r'[XY].?\d+.\d+', line) if coord: a = coord[0] b = coord[1] aa = float(a[1:]) bb = float(b[1:]) c = [aa,bb,i] coords.append(c) elif "layer_num" in line: coordinates.append(coords) coords = [] i=i+1 ####calculate the shortest distance between a point and line(defined by two points) P1 = coordinates[1] #first layer is not coordinates[0] as that is blank P2 = coordinates[len(coordinates)-1]#last layer start = time.time() @jit def distance(p, q, rs): x = p-q return np.linalg.norm( np.outer(np.dot(rs-q, x)/np.dot(x, x), x)+q-rs, axis=1) @jit def loops3(l,j,k): m = 0 #interior coordinate iterator min_dist = 99999999.9 #some large distance for each point in a layer while m < len(coordinates[l]): p = np.array(P1[j]) #coord of top profile q = np.array(P2[k]) #coord of bottom profile rs = np.array(coordinates[l][m])#interior layer at the m-th point d = distance(p, q, rs) if min_dist > d: #first time around min_dist should be large number, so this will be true min_dist = d #distance to the closest coordinate [m] on a layer [l] #print("min_dist " +str(min_dist)) m = m + 1 #next coordinate on the same layer [l] return(min_dist) @jit def loops2(j,k): l = 2 #interior layer iterator. Starts at 2 because first layer is empty [], and the first layer is P1 sum_min_dist = 0 while l < len(coordinates)-1: #don't count the first layer as it is blank [] and don't use P1, start at coordinates[2] to last coordinate -1 #end result is min_dist is equal to the distance to the closest coordinate min_dist = loops3(l,j,k) sum_min_dist= sum_min_dist + min_dist #sum of the closest coordinate distances on each layer l = l +20 #next layer. You can skip layers by making l+ some number of skipped layers return(sum_min_dist) @jit def loops(j): k = 0 #bottom profile iterator lowest_dist = 9999999.9#some large distance for total of layer's minimum distance while k < len(P2): synchro_dist=loops2(j,k) if synchro_dist < lowest_dist: #if the sum of the closest distances for all layers is less than lowest_dist = synchro_dist #the previous lowest sum of closest distance, then it may be a synch point synch_coord = P2[k] #current lowest coordinate on P2 through k coordinates on P2 k = k + 1 #try the next coordinate on P2 return(synch_coord) j = 0 #top profile coordinate iterator while j < len(P1): synchro_coords = loops(j) print(str(P1[j]) +str(synchro_coords)) j = j + 1 end = time.time() print(end - start)
I came up with a new algorithm that runs faster! It moves along both P1 and P2 depending on which movement has a better minimum synchronization distance.
- I updated the code to use the long algorithm for the first coordinate only, i.e., check P1[1] versus all P2[1 through n], to find the first set of synchronized coordinates.
- The P2 list is reordered starting from the synchronized coordinate.
- We check the synchronization distance for three cases: the next point on P2 to the next point on P1, the next point on P1 to the current point on P2, and the next point on P2 to the current point on P1. Whichever has the best synchronization distance, use that combination as a synchronization coordinate.
I had a problem with the algorithm getting stuck in a local minimum, and the same coordinate would keep repeating. In an attempt to fix this, I added code to check if the last three coordinates on either P1 or P2 are repeats. If repeats on P1, then step to the next one or two on P1 (depending on which gives you the better synchronization distance), and vice versa for repeats on P2.
#Read in P2 to a new P2* starting at location k n=0 PP2=[] lP = len(P2)-kopt # #print(lP) while n < lP: PP2.append(P2[kopt+n]) n=n+1 n=0 while n < kopt: PP2.append(P2[n]) n=n+1 P2=PP2 Q=0 R=1 w=[1,2,3] #w tells you how many times you've stepped the P1 side in a row x=[1,2,3] #x tells you how many times you've stepped the P2 side in a row U=[] while Q<len(P1)-2: while R<len(P2)-1: T=loops2(Q,R+1) #iterate along P2 U=loops2(Q+1,R) #iterate along P1 V=loops2(Q+1,R+1)#iterate along both P2 and P1 W=loops2(Q,R+2) X=loops2(Q+2,R) if x[0]==x[1]==x[2]: #check if your last four points on P2 are the same T=loops2(Q,R+1) W=loops2(Q,R+2) if T < W: x.pop(0) x.append(P2[R+1]) w.pop(0) w.append(P1[Q]) file.write((str(P1[Q]) +" ; " +str(P2[R+1]) +'\n')) R = R+1 else: x.pop(0) x.append(P2[R+2]) w.pop(0) w.append(P1[Q]) file.write((str(P1[Q]) +" ; " +str(P2[R+2]) +'\n')) R = R+2 elif w[0]==w[1]==w[2]:#check if your last four points on P1 are the same X=loops2(Q+2,R) U=loops2(Q+1,R) if U<X: x.pop(0) x.append(P2[R]) w.pop(0) w.append(P1[Q+1]) file.write((str(P1[Q+1]) +" ; " +str(P2[R])+'\n')) Q = Q+1 else: x.pop(0) x.append(P2[R]) w.pop(0) w.append(P1[Q+2]) file.write((str(P1[Q+2]) +" ; " +str(P2[R])+'\n')) Q = Q+2 else: if U<V and U<T: x.pop(0) x.append(P2[R]) w.pop(0) w.append(P1[Q+1]) file.write((str(P1[Q+1]) +" ; " +str(P2[R])+'\n')) Q=Q+1 elif T<U and T<V: x.pop(0) x.append(P2[R+1]) w.pop(0) w.append(P1[Q]) file.write((str(P1[Q]) +" ; " +str(P2[R+1])+'\n')) R=R+1 else: x.pop(0) x.append(P2[R+1]) w.pop(0) w.append(P1[Q+1]) file.write((str(P1[Q+1]) +" ; " +str(P2[R+1])+'\n')) Q = Q+1 R = R+1
So how did it do? Here is a slice of the back end of a Lamborghini Gallardo, and the GCODE used to find synchronization coordinates.
You can see the coordinates that have been synchronized in the picture below. Left (P1) and Right (P2) start at slightly different locations. At some points, the coordinates look fairly well synchronized, but further on, they are way out of synch. Sometimes P1 is too far along its path, and sometimes it’s P2.
