Homec4science
Diffusion Roughness_ENAC_Project (master) Browse

R10592/BOX COUNTING/BoxFragmentsAmplified.py
86c9f01a7d54master

Last Change
Raw File

BoxFragmentsAmplified.py
View Options

import numpy as num
import matplotlib.pyplot as plt
import math
import time
import csv
from mpl_toolkits.mplot3d import Axes3D
#The code calculate the number of boxes intersected for the fragments with amplified roughness
#####################################
with open('Fragment1_full.txt','r') as f:
lines = list(csv.reader(f, delimiter = ' ', skipinitialspace = True))
x =num.zeros(len(lines));
y = num.zeros(len(lines));
z = num.zeros(len(lines));
for i in range(len(lines)):
x[i] = lines[i][0];
y[i] = lines[i][1];
z[i] = lines[i][2];
#creo lista e appendo lo step 0
with open('Fragment1_full_smoothed.txt','r') as f:
lines = list(csv.reader(f, delimiter = ' ', skipinitialspace = True))
dnew =num.zeros(len(lines));
for i in range(len(lines)):
dnew[i] = lines[i][4];
#calculation of the centroid of the object
Cx=num.mean(x)
Cy=num.mean(y)
Cz=num.mean(z)
#initialization of new coordinates
xc = num.zeros(len(lines));
yc = num.zeros(len(lines));
zc = num.zeros(len(lines));
#Translation of the object to have the centroid at the origin (Cx = Cy = Cz = 0)
for i in range(len(lines)):
xc[i] = x[i] - Cx
yc[i] = y[i] - Cy
zc[i] = z[i] - Cz
d=num.zeros(len(lines)); #distance from centroid
lat=num.zeros(len(lines)); #latitude
lon=num.zeros(len(lines)); #longitude
for i in range(len(lines)):
d[i] = math.sqrt(xc[i]**2 + yc[i]**2 + zc[i]**2)
lat[i] = math.degrees(math.acos(zc[i]/d[i]))
lon[i] = math.degrees(math.atan2(yc[i], xc[i]))
dmean=num.mean(d) #mean of distances
dmin=min(d)
dmax=max(d)
dmindiff=dmean-dmin; #difference between original and smoothed
dmaxdiff=dmax-dmean;
damp=num.zeros(len(lines)) #amplification
coeff=num.zeros(len(lines))
for i in range(len(lines)):
coeff[i]=1+abs(d[i]-dnew[i])
damp[i]=d[i]*coeff[i]
xamp=num.zeros(len(lines)) #coordinates with amplification
yamp=num.zeros(len(lines))
zamp=num.zeros(len(lines))
xnew=num.zeros(len(lines)) #coordinates of smoothed fragment
ynew=num.zeros(len(lines))
znew=num.zeros(len(lines))
for i in range(len(lines)):
xamp[i] = damp[i] * math.sin(math.radians(lat[i])) * math.cos(math.radians(lon[i]))
yamp[i] = damp[i] * math.sin(math.radians(lat[i])) * math.sin(math.radians(lon[i]))
zamp[i] = damp[i] * math.cos(math.radians(lat[i]))
xnew[i] = dnew[i] * math.sin(math.radians(lat[i])) * math.cos(math.radians(lon[i]))
ynew[i] = dnew[i] * math.sin(math.radians(lat[i])) * math.sin(math.radians(lon[i]))
znew[i] = dnew[i] * math.cos(math.radians(lat[i]))
XBoxmin=1.0000001*min(min(xamp), min(yamp), min(zamp)); #Min. coordinate of the first box
XBoxmax=1.0000001*max(max(xamp), max(yamp), max(zamp)); #Max. coordinate of the first box
print('number of points:')
print(len(lines))
print('Coordinates min. and max. of the first box:')
print(XBoxmin)
print(XBoxmax)
LBox = XBoxmax - XBoxmin; #Side length of the first box
print('Length of the side of the first box:')
print(LBox)
steps=18; #Number of steps for the computation
#sf=num.zeros(steps); #scaling factors
sf=[5,7,10,15,20,25,30,40,50,60,70,80,100,120,150,200,250,300];
NBox=num.zeros(steps,dtype=int); #nr of total boxes for each scaling factor
for i in range(steps):
NBox[i]=int(sf[i]**3);
print('scaling factors:')
print(sf)
print('Nr. of total boxes:')
print(NBox)
NN=num.zeros(steps,dtype=int);
t = time.time()
for n in range(0,steps): #loop scaling factor
lbox=LBox/sf[n]
array=num.zeros(NBox[n],dtype=int)
for i in range(0,len(lines)):
ii=int( (xamp[i]-XBoxmin)/lbox)
ij=int( (yamp[i]-XBoxmin)/lbox)
ik=int( (zamp[i]-XBoxmin)/lbox)
ind=int(ii+sf[n]*ij+sf[n]*sf[n]*ik)
array[ind]=array[ind]+1
for s in range(0,NBox[n]):
if array[s]>0:
NN[n]=NN[n]+1
####################
#To export x,y,z coo.
Xamp=num.zeros((len(lines),3));
for i in range(len(lines)):
Xamp[i, 0]=xamp[i]
Xamp[i, 1] = yamp[i]
Xamp[i, 2] = zamp[i]
Xnew=num.zeros((len(lines),3));
for i in range(len(lines)):
Xnew[i, 0]=xnew[i]
Xnew[i, 1] = ynew[i]
Xnew[i, 2] = znew[i]
num.savetxt('Xamp.txt', Xamp, delimiter=' ', fmt='%s')
num.savetxt('Xnew.txt', Xnew, delimiter=' ', fmt='%s')
print('Nr of boxes intersected:')
print(NN)
elapsed = time.time() - t
print('time elapsed [s]:')
print(elapsed)
#plot amplified fragment
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(xamp, yamp, zamp, c='r', marker='o')
plt.show()
c4science · Help