https://github.com/jonasmb/curlnoisejittering
Tip revision: 42ba8b5ee132682b7804c6f47a7e1b91d6d73c9d authored by Jonàs Martínez on 26 February 2024, 12:21:49 UTC
Add Worley noise image 250x250 pixels
Add Worley noise image 250x250 pixels
Tip revision: 42ba8b5
psa_wrapper.py
import subprocess
import os
import numpy
import shutil
import matplotlib
import matplotlib.pyplot as plt
def read_psa_raw_data(filename):
field_1 = []
field_2 = []
with open(filename, 'r') as f:
for line in f:
fields = line.strip('\n').split(' ')
assert len(fields) == 2
field_1.append(float(fields[0]))
field_2.append(float(fields[1]))
return (field_1, field_2)
def psa_wrapper(samples, temp_post_fix = "", points_id=0):
psabin_path = os.path.join(os.path.dirname(__file__), 'psa_mod/build/psabin')
temp_path = os.path.join(os.path.dirname(__file__), 'temp' + temp_post_fix)
if not os.path.exists(temp_path):
os.makedirs(temp_path)
print(os.path.curdir, os.path.exists(psabin_path))
assert os.path.exists(psabin_path) or os.path.exists(psabin_path + '.exe')
points_name = 'points' + str(points_id) + '_'
points_path = os.path.join(temp_path, points_name)
for i in range(len(samples)):
points = samples[i]
points_psa_txt = points_path + str(i) + '.txt'
with open(points_psa_txt, 'w') as f:
f.write(str(len(points)) + '\n')
for point in points:
f.write(str(point[0]) + ' ' + str(point[1]) + '\n')
average = len(samples) > 1
flags = ' '
if average:
flags +='--avg'
flags += ' --raw --ani --rp --rdf --spectral '
if average:
path_ani = os.path.join(os.path.dirname(__file__), 'avg_ani.txt')
else:
path_ani = os.path.join(os.path.dirname(__file__), points_name + '0_ani.txt')
if os.path.exists(path_ani):
os.remove(path_ani)
if average:
path_rdf = os.path.join(os.path.dirname(__file__), 'avg_rdf.txt')
else:
path_rdf = os.path.join(os.path.dirname(__file__), points_name + '0_rdf.txt')
if os.path.exists(path_rdf):
os.remove(path_rdf)
if average:
path_rp = os.path.join(os.path.dirname(__file__), 'avg_rp.txt')
else:
path_rp = os.path.join(os.path.dirname(__file__), points_name + '0_rp.txt')
if os.path.exists(path_rp):
os.remove(path_rp)
if average:
path_spectral = os.path.join(os.path.dirname(__file__), 'avg_spectral.txt')
else:
path_spectral = os.path.join(os.path.dirname(__file__), points_name + '0_spectral.txt')
if os.path.exists(path_spectral):
os.remove(path_spectral)
if average:
cmd_str = psabin_path + flags + points_path + '*.txt'
else:
cmd_str = psabin_path + flags + points_path + '0.txt'
print(subprocess.run(cmd_str, shell=True))
# read back data from files
with open(path_spectral, 'r') as f:
lines = f.readlines()
assert len(lines) == 2
effnyquist = float(lines[0].strip('\n'))
oscillations = float(lines[1].strip('\n'))
print("effnyquist = " + str(effnyquist))
print("oscillations = " + str(oscillations))
rp_raw_data = read_psa_raw_data(path_rp)
ani_raw_data = read_psa_raw_data(path_ani)
rdf_raw_data = read_psa_raw_data(path_rdf)
os.remove(path_ani)
os.remove(path_rdf)
os.remove(path_rp)
os.remove(path_spectral)
if not average:
#os.remove(os.path.join(os.path.dirname(__file__), points_name + '0_spec.png'))
os.remove(points_path + '0.txt')
return (
effnyquist,
oscillations,
rp_raw_data,
ani_raw_data,
rdf_raw_data,
)
def plot_data(data, xlabel, ylabel, name, to_pgf=False):
if to_pgf:
matplotlib.use("pgf") # for latex (see https://ctan.org/pkg/pgf)
matplotlib.rcParams.update({
"pgf.texsystem": "pdflatex",
'font.family': 'serif',
'text.usetex': True,
'pgf.rcfonts': False,
})
fig, ax = plt.subplots()
ax.plot(data[0], data[1], c='k')
ax.set(xlabel=xlabel, ylabel=ylabel)
ax.set_xlim(left=0.0)
if to_pgf:
fig.savefig(filename + ".pgf")
else:
plt.show()
if __name__ == '__main__':
numpy.random.seed(42)
samples = []
num_points = 128
num_samples = 1
for i in range(num_samples):
samples.append(numpy.random.random_sample((num_points, 2)))
print(len(samples), samples[0].shape)
(effnyquist, oscillations, data_rp, data_ani, data_rdf) = psa_wrapper(samples)
plot_data(data_rp, 'Frequency', 'Radial Power', 'plot_rp')
plot_data(data_ani, 'Frequency', 'Anisotropy', 'plot_ani')
plot_data(data_rdf, 'Distance', 'RDF', 'plot_rdf')
# heuristic measure of anisotropy derived from raw data PSA
print("anisotropy heuristic mean =" + str(numpy.mean(data_ani[1][1:])))