Homec4science
Diffusion Roughness_ENAC_Project (master) Browse

R10592/SPHARM/PSD_Fragment1_diffTypesOfSmoothing.py
e6c5d1a9727emaster

Last Change
Raw File

PSD_Fragment1_diffTypesOfSmoothing.py
View Options

from matplotlib import pyplot as plt
import numpy as np
import math
import pyshtools as pysh
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import csv
import time
#Calculate the PSD of fragment with different types of smoothing and PSD of the difference between original fragments and smoothed versions
t = time.time()
#import original fragment
with open('Fragment1.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][1];
y[i] = lines[i][2];
z[i] = lines[i][3];
dim=len(lines) #number of points
#To center the object
#calculation of the centroid of the objects
Cx=num.mean(x)
Cy=num.mean(y)
Cz=num.mean(z)
#Translation of the object to have the centroid at the origin (Cx = Cy = Cz = 0)
for i in range(dim):
xc[i] = x[i] - Cx
yc[i] = y[i] - Cy
zc[i] = z[i] - Cz
d=num.zeros(dim); #distance from centroid to surface
lat=num.zeros(dim); #latitude
lon=num.zeros(dim); #longitude
#conversion in spherical coordinates
for i in range(dim):
d[i] = math.sqrt(xc[i]**2 + yc[i]**2 + zc[i]**2)
lat[i] = math.degrees(math.acos(zc[i]/d[i]))-90 #SHTools package take -90° < lat < 90°
lon[i] = math.degrees(math.atan2(yc[i], xc[i])) #and 0 < lon < 360°
####################
#smoothed fragment from Blender, factor 1
with open('Fragment1_reduced_Factor1.txt','r') as f:
linesf1 = list(csv.reader(f, delimiter = ' ', skipinitialspace = True))
xf1 =num.zeros(len(linesf1));
yf1 = num.zeros(len(linesf1));
zf1 = num.zeros(len(linesf1));
for i in range(len(linesf1)):
xf1[i] = linesf1[i][1];
yf1[i] = linesf1[i][2];
zf1[i] = linesf1[i][3];
dimf1=len(linesf1)
#calculation of the centroid of the objects
Cxf1=num.mean(xf1)
Cyf1=num.mean(yf1)
Czf1=num.mean(zf1)
#Translation of the object to have the centroid at the origin (Cx = Cy = Cz = 0)
for i in range(dimf1):
xf1c[i] = xf1[i] - Cxf1
yf1c[i] = yf1[i] - Cyf1
zf1c[i] = zf1[i] - Czf1
df1=num.zeros(dimf1); #distance from centroid to surface
latf1=num.zeros(dimf1); #latitude
lonf1=num.zeros(dimf1); #longitude
#conversion in spherical coordinates
for i in range(dimf1):
df1[i] = math.sqrt(xf1c[i]**2 + yf1c[i]**2 + zf1c[i]**2)
latf1[i] = math.degrees(math.acos(zf1c[i]/df1[i]))-90 #SHTools package take -90° < lat < 90°
lonf1[i] = math.degrees(math.atan2(yf1c[i], xf1c[i])) #and 0 < lon < 360°
####################
#smoothed fragment from Blender, factor 1
with open('Fragment1_reduced_Factor2.txt','r') as f:
linesf2 = list(csv.reader(f, delimiter = ' ', skipinitialspace = True))
xf2 =num.zeros(len(linesf2));
yf2 = num.zeros(len(linesf2));
zf2 = num.zeros(len(linesf2));
for i in range(len(linesf2)):
xf2[i] = linesf2[i][1];
yf2[i] = linesf2[i][2];
zf2[i] = linesf2[i][3];
dimf2=len(linesf2)
#calculation of the centroid of the objects
Cxf2=num.mean(xf2)
Cyf2=num.mean(yf2)
Czf2=num.mean(zf2)
#Translation of the object to have the centroid at the origin (Cx = Cy = Cz = 0)
for i in range(dimf2):
xf2c[i] = xf2[i] - Cxf2
yf2c[i] = yf2[i] - Cyf2
zf2c[i] = zf2[i] - Czf2
df2=num.zeros(dimf2); #distance from centroid to surface
latf2=num.zeros(dimf2); #latitude
lonf2=num.zeros(dimf2); #longitude
#conversion in spherical coordinates
for i in range(dimf2):
df2[i] = math.sqrt(xf2c[i]**2 + yf2c[i]**2 + zf2c[i]**2)
latf2[i] = math.degrees(math.acos(zf2c[i]/df2[i]))-90 #SHTools package take -90° < lat < 90°
lonf2[i] = math.degrees(math.atan2(yf2c[i], xf2c[i])) #and 0 < lon < 360°
####################
#cone smoothed fragment (from smoothing.cpp)
with open('Fragment1_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];
Diff1=num.zeros(dim);
Diff2=num.zeros(dim);
Diffnew=num.zeros(dim);
for i in range(len(lines)):
Diff1[i] = d[i] - df1[i] #difference of d between original and smoothed with blender with factor 1
Diff2[i] = d[i] - df2[i] #difference of d between original and smoothed with blender with factor 2
Diffnew[i] = d[i] - dnew[i] #difference of d between original and smoothed with smoothing.cpp
norm=1; #Geodesy 4-pi normalized harmonics
csphase=1; #do not apply the Condon-Shortley phase factor to the associated Legendre functions
lmax=30; #Maximum number of the spherical harmonic degree l
#Compute the spherical harmonic coefficients cilm
cilm, chi2 = pysh.expand.SHExpandLSQ(d, lat, lon, lmax, [norm, csphase]); #original
cilmf1, chi2 = pysh.expand.SHExpandLSQ(df1, latf1, lonf1, lmax, [norm, csphase]); #smoothed with blender (factor 1)
cilmf2, chi2 = pysh.expand.SHExpandLSQ(df2, latf2, lonf2, lmax, [norm, csphase]); #smoothed with blender (factor 2)
cilmdnew, chi2 = pysh.expand.SHExpandLSQ(dnew, lat, lon, lmax, [norm, csphase]); #smoothed with smoothing.cpp
cilmDiff1, chi2 = pysh.expand.SHExpandLSQ(Diff1, lat, lon, lmax, [norm, csphase]); #difference original and blender (factor 1)
cilmDiff2, chi2 = pysh.expand.SHExpandLSQ(Diff2, lat, lon, lmax, [norm, csphase]); #difference original and blender (factor 2)
cilmDiffnew, chi2 = pysh.expand.SHExpandLSQ(Diffnew, lat, lon, lmax, [norm, csphase]); #difference original and smoothing.cpp
cilmcoeff1=cilmf1-cilm; #difference applied on coefficients
cilmcoeff2=cilmf2-cilm;
cilmcoeffnew=cilmdnew-cilm;
#Calculate the power spectrum
P = pysh.spectralanalysis.spectrum(cilm);
degrees = num.arange(cilm.shape[1]);
Pf1 = pysh.spectralanalysis.spectrum(cilmf1);
degreesf1 = num.arange(cilmf1.shape[1]);
Pf2 = pysh.spectralanalysis.spectrum(cilmf2);
degreesf2 = num.arange(cilmf2.shape[1]);
Pdnew = pysh.spectralanalysis.spectrum(cilmdnew);
degreesdnew = num.arange(cilmdnew.shape[1]);
PDiff1 = pysh.spectralanalysis.spectrum(cilmDiff1);
degreesDiff1 = num.arange(cilmDiff1.shape[1]);
PDiff2 = pysh.spectralanalysis.spectrum(cilmDiff2);
degreesDiff2 = num.arange(cilmDiff2.shape[1]);
PDiffnew = pysh.spectralanalysis.spectrum(cilmDiffnew);
degreesDiffnew = num.arange(cilmDiffnew.shape[1]);
Pcoeff1 = pysh.spectralanalysis.spectrum(cilmcoeff1);
degreescoeff1 = num.arange(cilmcoeff1.shape[1]);
Pcoeff2 = pysh.spectralanalysis.spectrum(cilmcoeff2);
degreescoeff2 = num.arange(cilmcoeff2.shape[1]);
Pcoeffnew = pysh.spectralanalysis.spectrum(cilmcoeffnew);
degreescoeffnew = num.arange(cilmcoeffnew.shape[1]);
elapsed = time.time() - t
print('time elapsed [s]:')
print(elapsed)
#To export Powers
np.savetxt('Pdnew.txt', Pdnew, delimiter=' ', fmt='%s')
np.savetxt('PDiff1.txt', PDiff1, delimiter=' ', fmt='%s')
np.savetxt('PDiff2.txt', PDiff2, delimiter=' ', fmt='%s')
np.savetxt('PDiffnew.txt', PDiffnew, delimiter=' ', fmt='%s')
np.savetxt('Pf1.txt', Pf1, delimiter=' ', fmt='%s')
np.savetxt('Pf2.txt', Pf2, delimiter=' ', fmt='%s')
#Plot power spectrum
fig, ax = plt.subplots(1, 1)
plt.loglog(degrees, P, label='Original from centroid')
plt.loglog(degreescoeff1, Pcoeff1, label='coeff1')
plt.loglog(degreescoeff2, Pcoeff2, label='coeff2')
plt.loglog(degreescoeffnew, Pcoeffnew, label='coeffnew')
plt.loglog(degreesdnew, Pdnew, label='newset from centroid')
plt.loglog(degreesDiff1, PDiff1,'--', label='Difference from F1')
plt.loglog(degreesDiff2, PDiff2,'--', label='Difference from F2')
plt.loglog(degreesDiffnew, PDiffnew,'--', label='Difference from newset')
plt.loglog(degreesf1, Pf1, label='F1 from centroid')
plt.loglog(degreesf2, Pf2, label='F2 from centroid')
ax.set(xlabel='Spherical harmonic degree', ylabel='Power')
ax.grid()
plt.legend(loc="upper right")
plt.show()
c4science · Help