Homec4science
Diffusion Roughness_ENAC_Project (master) Browse

R10592/BOX COUNTING/BoxCubeGauss.py
86c9f01a7d54master

Last Change
Raw File

BoxCubeGauss.py
View Options

import numpy as num
import matplotlib.pyplot as plt
import math
import time
import csv
from mpl_toolkits.mplot3d import Axes3D
import random
from random import randint
#The code calculate the number of boxes intersected for chosen scaling factors for cube with gaussian noise
#####################################
mm=100000; #number of points per cube side
Maxcoo=1; #max. coordinates of the cube
Mincoo=-1; #min. coordinates of the cube
#gaussian noise (different seed to have non identical cube faces)
random.seed(1)
noise1 = num.random.normal(0, 0.05, mm);
random.seed(2)
noise2 = num.random.normal(0, 0.05, mm);
random.seed(3)
noise3 = num.random.normal(0, 0.05, mm);
random.seed(4)
noise4 = num.random.normal(0, 0.05, mm);
random.seed(5)
noise5 = num.random.normal(0, 0.05, mm);
random.seed(6)
noise6 = num.random.normal(0, 0.05, mm);
#coordinates of points on face nr. 1
x1=num.zeros(mm);
y1=num.zeros(mm);
z1=num.zeros(mm);
for i in range(mm):
x1[i] = randint(1000*Mincoo, 1000*Maxcoo)/1000
y1[i] = randint(1000*Mincoo, 1000*Maxcoo)/1000
z1[i] = Mincoo+noise1[i];
#coordinates of points on face nr. 2
x2=num.zeros(mm);
y2=num.zeros(mm);
z2=num.zeros(mm);
for i in range(mm):
x2[i] = randint(1000*Mincoo, 1000*Maxcoo)/1000
y2[i] = randint(1000*Mincoo, 1000*Maxcoo)/1000
z2[i] = Maxcoo+noise2[i];
x3=num.zeros(mm);
y3=num.zeros(mm);
z3=num.zeros(mm);
for i in range(mm):
x3[i] = Mincoo+noise3[i]
y3[i] = randint(1000*Mincoo, 1000*Maxcoo)/1000
z3[i] = randint(1000*Mincoo, 1000*Maxcoo)/1000
x4=num.zeros(mm);
y4=num.zeros(mm);
z4=num.zeros(mm);
for i in range(mm):
x4[i] = Maxcoo+noise4[i]
y4[i] = randint(1000*Mincoo, 1000*Maxcoo)/1000
z4[i] = randint(1000*Mincoo, 1000*Maxcoo)/1000
x5=num.zeros(mm);
y5=num.zeros(mm);
z5=num.zeros(mm);
for i in range(mm):
x5[i] = randint(1000*Mincoo, 1000*Maxcoo)/1000
y5[i] = Mincoo+noise5[i]
z5[i] = randint(1000*Mincoo, 1000*Maxcoo)/1000
x6=num.zeros(mm);
y6=num.zeros(mm);
z6=num.zeros(mm);
for i in range(mm):
x6[i] = randint(1000*Mincoo, 1000*Maxcoo)/1000
y6[i] = Maxcoo+noise6[i]
z6[i] = randint(1000*Mincoo, 1000*Maxcoo)/1000
xc=[];
yc=[];
zc=[];
xc = num.append(xc, x1)
xc = num.append(xc, x2)
xc = num.append(xc, x3)
xc = num.append(xc, x4)
xc = num.append(xc, x5)
xc = num.append(xc, x6)
yc = num.append(yc, y1)
yc = num.append(yc, y2)
yc = num.append(yc, y3)
yc = num.append(yc, y4)
yc = num.append(yc, y5)
yc = num.append(yc, y6)
zc = num.append(zc, z1)
zc = num.append(zc, z2)
zc = num.append(zc, z3)
zc = num.append(zc, z4)
zc = num.append(zc, z5)
zc = num.append(zc, z6)
XBoxmin=1.0000001*min(min(xc), min(yc), min(zc)); #Min. coordinate of the first box
XBoxmax=1.0000001*max(max(xc), max(yc), max(zc)); #Max. coordinate of the first box
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=14; #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];
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,mm):
ii=int( (xc[i]-XBoxmin)/lbox)
ij=int( (yc[i]-XBoxmin)/lbox)
ik=int( (zc[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
print('Nr of boxes intersected:')
print(NN)
elapsed = time.time() - t
print('time elapsed [s]:')
print(elapsed)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(xc, yc, zc, c='r', marker='o')
plt.show()
c4science · Help