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
	filelist = csv.reader(csvfile1,delimiter="\t")
	for i, row in enumerate(filelist):
		boxname = "stack5/" + row[0]
		#print i
		with open(boxname, 'r') as csvfile2:
			boxfile = csv.reader(csvfile2,delimiter="\t")
			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.

Delete this entry (admin only).