Homec4science
Diffusion Roughness_ENAC_Project (master) Browse

R10592/BOX COUNTING/BoxSphereGauss.py
86c9f01a7d54master

Last Change
Raw File

BoxSphereGauss.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 a sphere with gaussian noise
random.seed(10)
m = 10000; # number of points of the object (sphere in this case)
d = num.ones(m)
lat = num.zeros(m); # latitude in degrees, initialization
lon = num.zeros(m); # longitude in degrees, initialization
noise = num.random.normal(0, 0.05, m);
print(min(noise))
print(max(noise))
for i in range(m):
lat[i] = randint(-90, 90) # generate random latitudes and longitudes angles (theta and phi)
lon[i] = randint(-180, 180)
d[i] = d[i] + noise[i]
x = num.zeros(m); # initialization of the cartesian coordinates
y = num.zeros(m);
z = num.zeros(m);
for i in range(m): # conversion from spherical coordinates into cartesian coo.
x[i] = d[i] * math.sin(math.radians(lon[i])) * math.cos(math.radians(lat[i]))
y[i] = d[i] * math.sin(math.radians(lon[i])) * math.sin(math.radians(lat[i]))
z[i] = d[i] * math.cos(math.radians(lon[i]))
XBoxmin=1.0000001*min(min(x), min(y), min(z)); #Min. coordinate of the first box
XBoxmax=1.0000001*max(max(x), max(y), max(z)); #Max. coordinate of the first box
print('number of points:')
print(m)
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,m):
ii=int( (x[i]-XBoxmin)/lbox)
ij=int( (y[i]-XBoxmin)/lbox)
ik=int( (z[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 coordinates of sphere with noise
XGauss=num.zeros((m,3));
for i in range(m):
XGauss[i, 0] = x[i]
XGauss[i, 1] = y[i]
XGauss[i, 2] = z[i]
num.savetxt('XGauss.txt', XGauss, delimiter=' ', fmt='%s')
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(x, y, z, c='r', marker='o')
plt.show()
c4science · Help