# Demo entry 6623978

tubes

Submitted by anonymous on Jun 12, 2017 at 14:40
Language: Python 3. Code size: 2.4 kB.

```import numpy as np
import csv
import codecs
import glob

tube_num = 0
with open('box-list.csv', 'r') as csvfile1: #reading the list of box files created by ls -1 stack5 > box-list.csv
for i, row in enumerate(filelist):
boxname = "stack5/" + row[0]
#print i
with open(boxname, 'r') as csvfile2:
px = np.zeros(300)
py = np.zeros(300)
for j,line in enumerate(boxfile): #reading points one by one in each file
px[j] = int(line[0])
py[j] = int(line[1])
particle_num = j

vx = np.zeros(particle_num)
vy = np.zeros(particle_num)
direction = np.zeros(particle_num)
for k in range(particle_num): #here I get the vector that connect each 2 consecutive points
tube_num = tube_num + 1
vx[k] = px[k+1] - px[k]
vy[k] = py[k+1] - py[k]
for k in range(particle_num-1):# here I get the dot product of each two consecutive vectors
direction[k] =  (vx[k] * vx[k+1] + vy[k] * vy[k+1])/(np.sqrt(vx[k]**2 + vy[k]**2)* np.sqrt(vx[k+1]**2 + vy[k+1]**2))
#				print direction[k]
temp1 = 0
boxnew = "stack-tubes/" + boxname[7:60]#output name
nstart = 0
nend = 0
with open(boxnew, 'w') as outfile:
for k in range(particle_num-1):
temp2 = int(direction[k] > 0.95) #check if 2 vectors are parallel, i.e. on the same tube
turn = temp2 - temp1#have we "turned" or not! 1: start of new filament, 0: same filament, -1: end of the old filament
if turn == 1:
pstart = k
#print "start = %d" %(pstart)
if turn == -1:
pend = k + 1
#print  "end = %d" %(pend)
if pend - pstart > 1:
pointstr = "%d\t%d\n" %(px[pstart],py[pstart])
outfile.write(pointstr)
pointstr = "%d\t%d\n" %(px[pend],py[pend])
outfile.write(pointstr)
nend = nend + 1
nstart = nstart + 1
#often the file ends but the filament is there, so turn!= -1, for that:
if k == particle_num-2 and turn == 0 and nstart==nend+1:
pend = k + 2
#print  "endf = %d" %(k)
pointstr = "%d\t%d\n" %(px[k],py[k])
outfile.write(pointstr)
nend = nend + 1
temp1 = temp2
#del px,py,vx, vy, direction
if nstart!=nend:#make sure all filaments start and end!
print "warning %s %d" %(boxnew,nstart-nend)

```

This snippet took 0.01 seconds to highlight.

Back to the Entry List or Home.