Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add simple test cylinder #52

Merged
merged 5 commits into from
Sep 19, 2023
Merged

Add simple test cylinder #52

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
5c7f516
Add simple test cylinder
Sep 13, 2023
585567a
Modify cylinder.py to follow repo's flake8 configuration
Sep 15, 2023
89bb101
Add commenting in postprocessing_h5py to indicate author of common fu…
Sep 15, 2023
6fde7a8
Couple more fixes to cylinder.py to satisfy mypy
Sep 15, 2023
afbd336
More reformatting cylinder.py
Sep 18, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
16 changes: 16 additions & 0 deletions postprocessing_h5py/postprocessing_common_h5py.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,9 @@

"""
This script contains a number of helper functions to create visualizations outsside of fenics.

All functions authored by David Bruneau unless specified

"""

def read_command_line():
Expand Down Expand Up @@ -137,6 +140,7 @@ def get_interface_ids(meshFile):

def get_sampling_constants(df,start_t,end_t):
'''
Author: Daniel Macdonald
T = period, in seconds,
nsamples = samples per cycle
fs = sample rate
Expand All @@ -162,6 +166,10 @@ def read_csv_file(filepath):
return df

def filter_SPI(U, W_low_cut, tag):
"""
Author: Mehdi Najafi

"""
if tag=="withmean":
U_fft = fft(U)
else:
Expand All @@ -185,6 +193,10 @@ def filter_SPI(U, W_low_cut, tag):


def filter_SPI_print(U, W_low_cut, tag):
"""
Author: Mehdi Najafi

"""
if tag=="withmean":
U_fft = fft(U)
else:
Expand All @@ -208,6 +220,10 @@ def filter_SPI_print(U, W_low_cut, tag):
return Power_25Hz/Power_0Hz

def calculate_spi(case_name, df, output_folder, meshFile,start_t,end_t,low_cut,high_cut, dvp):
"""
Author: Mehdi Najafi and David Bruneau

"""
# cut, thresh
# Get wall and fluid ids
fluidIDs, wallIDs, allIDs = get_domain_ids(meshFile)
Expand Down
177 changes: 17 additions & 160 deletions postprocessing_h5py/spectrograms.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,7 @@
except:
print("Could not import pyvista and/or vtk, install these packages or use 'RandomPoint' sampling for spectrograms.")
"""
This library contains helper functions for creating spectrograms.
This library contains helper functions for creating spectrograms. All functions authored by David Bruneau unless otherwise specified.
"""

def read_spec_config(config_file,dvp):
Expand Down Expand Up @@ -210,19 +210,8 @@ def read_spectrogram_data(case_path, mesh_name, save_deg, stride, start_t, end_t
elif sampling_method == "Spatial": # This method only works with a specified sphere
mesh, surf = postprocessing_common_pv.assemble_mesh(mesh_path)

#bounds=surf.bounds
bounds = [x_sphere-r_sphere, x_sphere+r_sphere, y_sphere-r_sphere, y_sphere+r_sphere, z_sphere-r_sphere, z_sphere+r_sphere]
#def generate_points(bounds, subdivisions=50):
# x_points=np.linspace(bounds[0], bounds[1],num=subdivisions)
# y_points=np.linspace(bounds[2], bounds[3],num=subdivisions)
# z_points=np.linspace(bounds[4], bounds[5],num=subdivisions)
# points = np.array([[x_points[0], y_points[0], z_points[0]],[x_points[0], y_points[0], z_points[1]]])
# for i in range(subdivisions):
# for j in range(subdivisions):
# for k in range(subdivisions):
# points=np.append(points,[[x_points[i], y_points[j], z_points[k]]], axis=0)
# return points[2:,:]



def generate_points(bounds, subdivisions=50):
x_points=np.linspace(bounds[0], bounds[1],num=subdivisions)
Expand Down Expand Up @@ -254,10 +243,6 @@ def generate_points(bounds, subdivisions=50):
t = time.time()
print("shaped point cloud")
print(t - start)
#coords_sphere = coords[region_ids]
#print(coords_sphere)
#print("Length of coords array")
#print(coords.shape)
tree = KDTree(coords)
_, idx = tree.query(sphere_sel.points) #find closest node to the points in the equispaced points

Expand All @@ -266,21 +251,10 @@ def generate_points(bounds, subdivisions=50):
print("queried kdtree")
print(t - start)

#print("Unsampled equispaced data")
#print(len(idx))
#print(idx)
#

# Make sure this is the correct order
equispaced_fluid_ids = [x for x in idx if x in region_ids]
#equispaced_fluid_ids=list(set(region_ids).intersection(idx))
#print("Unsampled equispaced fluid data")
#print(len(equispaced_fluid_ids))
#print(equispaced_fluid_ids)

idx_sampled = np.random.choice(equispaced_fluid_ids, n_samples)
#print("Sampled equispaced fluid data")
#print(len(idx_sampled))
#print(idx_sampled)

idx_sampled_set = set(idx_sampled)
# compare the length and print if the list contains duplicates
Expand All @@ -290,17 +264,13 @@ def generate_points(bounds, subdivisions=50):
print("No duplicates found in the list")

#sampled_point_cloud=pv.PolyData(coords[idx_sampled])
##
#plotter= pv.Plotter(off_screen=True)
#

#plotter.add_mesh(surf,
# color='red',
# show_scalar_bar=False,
# opacity=0.05)
###
#plotter.add_points(sampled_point_cloud, render_points_as_spheres=True, point_size=0.03)
#
#
#plotter.show(auto_close=False)
#plotter.show(screenshot=imageFolder + "/points"+str(idx_sampled[0])+".png")

Expand All @@ -325,32 +295,6 @@ def get_location_from_pointID(probeID,polydata):
return nearestPointLoc


#def find_nearest_points_in_radius(probe_pos,pl,num_points,polyData,radius):
# """
# probeID is actual ID of nearest point
# nearestPointLoc is xyz location of probeID
# nearestNeighbours is all the points within the radius
# nearest_points is sample of nearest_Neighbours
# NOTE: nearest_points is likely unneccesary after implementing masking
# function.
# """
# probeID = pl.FindClosestPoint(probe_pos)
# nearestPointLoc = [10000,10000,10000]
# polyData.GetPoint(probeID,nearestPointLoc)
# nearestNeighbours = vtk.vtkIdList()
# pl.FindPointsWithinRadius(radius,probe_pos,nearestNeighbours)
# num_Ids = nearestNeighbours.GetNumberOfIds()
# print('There are',num_Ids,'in radius',radius)
# idList = []
# if num_Ids >= 1:
# if num_points > num_Ids:
# num_points = num_Ids
# random_index = random.sample(range(num_Ids), num_points-1)
# idList = [nearestNeighbours.GetId(idx) for idx in random_index]
# nearest_points = [probeID] + idList
# return nearest_points


def find_points_in_sphere(cent,rad,coords):

# Calculate vector from center to each node in the mesh
Expand All @@ -371,42 +315,9 @@ def find_points_in_sphere(cent,rad,coords):
return points_in_sphere



#def vtk_read_unstructured_grid(FILE):
# '''
# Reads vtu file. Create a more general reader, particularly for h5 mesh
# going forward.
# '''
# reader = vtk.vtkXMLUnstructuredGridReader()
# reader.SetFileName(FILE)
# reader.Update()
# return reader.GetOutput()

#def vtk_mask_points(polyData, max_points=1000):
# '''
# Randomly mask points.
# SetRandomModeType(2) uses stratified spatial sampling.
# '''
# ptMask = vtk.vtkMaskPoints()
# ptMask.SetInputData(polyData)
# ptMask.RandomModeOn()
# ptMask.SetRandomModeType(2)
# ptMask.SetMaximumNumberOfPoints(max_points)
# ptMask.Update()
# return ptMask.GetOutput()

#def get_original_ids(masked_polydata, polyData):
# n = masked_polydata.GetPoints().GetNumberOfPoints()
# original_ids = np.zeros(n)
# pll = vtk.vtkPointLocator()
# pll.SetDataSet(polyData)
# for pt in range(n):
# location = masked_polydata.GetPoints().GetPoint(pt)
# original_ids[pt] = pll.FindClosestPoint(location)
# return original_ids

def shift_bit_length(x):
'''
Author: Dan MacDonald
round up to nearest pwr of 2
https://stackoverflow.com/questions/14267555/find-the-smallest-power-of-2-greater-than-n-in-python
'''
Expand All @@ -433,6 +344,7 @@ def get_psd(dfNearest,fsamp,scaling="density"):

def get_spectrogram(dfNearest,fsamp,nWindow,overlapFrac,window,start_t,end_t, scaling='spectrum', interpolate = False):
'''
Author: Daniel Macdonald
Calculates spectrogram
input dfNearest is a pandas df of shape (num_points, num_timesteps)
fsamp is sampling frequency
Expand Down Expand Up @@ -475,20 +387,8 @@ def get_spectrogram(dfNearest,fsamp,nWindow,overlapFrac,window,start_t,end_t, sc
Pxx_mean[Pxx_mean<0] = 1e-16
return Pxx_mean, freqs, bins

def spectrogram_scaling(Pxx_mean,thresh_percent):
Pxx_scaled = np.log(Pxx_mean)
max_val = np.max(Pxx_scaled)
min_val = np.min(Pxx_scaled)
thresh = max_val-(max_val-min_val)*thresh_percent
print('Pxx_scaled max', max_val)
print('Pxx_scaled max', min_val)
print('Pxx threshold', thresh)

Pxx_threshold_indices = Pxx_scaled < thresh
Pxx_scaled[Pxx_threshold_indices] = thresh
return Pxx_scaled, max_val, min_val, thresh

def spectrogram_scaling_old(Pxx_mean,lower_thresh):
def spectrogram_scaling(Pxx_mean,lower_thresh):
# Author: Daniel Macdonald
Pxx_scaled = np.log(Pxx_mean)
max_val = np.max(Pxx_scaled)
min_val = np.min(Pxx_scaled)
Expand Down Expand Up @@ -523,25 +423,26 @@ def butter_bandpass(lowcut, highcut, fs, order=5, btype='band'):
return b, a

def butter_bandpass_filter(data, lowcut=25.0, highcut=15000.0, fs=2500.0, order=5, btype='band'):
# Author: Daniel Macdonald
b, a = butter_bandpass(lowcut, highcut, fs, order=order,btype=btype)
y = filtfilt(b, a, data)
return y

def filter_time_data(df,fs,lowcut=25.0,highcut=15000.0,order=6,btype='highpass'):
# Author: Daniel Macdonald
df_filtered = df.copy()
for row in range(df.shape[0]):
df_filtered.iloc[row] = butter_bandpass_filter(df.iloc[row],lowcut=lowcut,highcut=highcut,fs=fs,order=order,btype=btype)
return df_filtered

def compute_average_spectrogram(df, fs, nWindow,overlapFrac,window,start_t,end_t,thresh, scaling="spectrum", filter_data=False,thresh_method="new"):
# Author: Daniel Macdonald
if filter_data == True:
df = filter_time_data(df,fs)

Pxx_mean, freqs, bins = get_spectrogram(df,fs,nWindow,overlapFrac,window,start_t,end_t, scaling) # mode of the spectrogram
if thresh_method == "new":
if thresh_method == "old":
Pxx_scaled, max_val, min_val, lower_thresh = spectrogram_scaling(Pxx_mean,thresh)
elif thresh_method == "old":
Pxx_scaled, max_val, min_val, lower_thresh = spectrogram_scaling_old(Pxx_mean,thresh)
elif thresh_method == "log_only":
Pxx_scaled = np.log(Pxx_mean)
max_val = np.max(Pxx_scaled)
Expand All @@ -557,8 +458,6 @@ def compute_average_spectrogram(df, fs, nWindow,overlapFrac,window,start_t,end_t
return bins, freqs, Pxx_scaled, max_val, min_val, lower_thresh

def plot_spectrogram(fig1,ax1,bins,freqs,Pxx,ylim_,title=None,convert_a=0.0,convert_b=0.0,x_label=None,color_range=None):
#plt.figure(figsize=(14,7)) #fig size same as before


if color_range == None:
im = ax1.pcolormesh(bins, freqs, Pxx, shading = 'gouraud')#, vmin = -30, vmax = -4) # look up in matplotlib to add colorbar
Expand Down Expand Up @@ -586,37 +485,8 @@ def time_convert(x):
ax2.set_xlabel(x_label)


def plot_spectrogram_1ax(bins,freqs,Pxx,ylim_,title=None,path=None,convert_a=1.0,convert_b=0.0,x_label=None,color_range=None):
#plt.figure(figsize=(14,7)) #fig size same as before

bins = bins*convert_a+convert_b
fig1, ax1 = plt.subplots()
if color_range == None:
im = ax1.pcolormesh(bins, freqs, Pxx, shading = 'gouraud')#, vmin = -30, vmax = -4) # look up in matplotlib to add colorbar
else:
im = ax1.pcolormesh(bins, freqs, Pxx, shading = 'gouraud', vmin = color_range[0], vmax = color_range[1]) # look up in matplotlib to add colorbar


fig1.colorbar(im, ax=ax1)
#fig1.set_size_inches(8, 6) #fig1.set_size_inches(10, 7)

ax1.set_xlabel(x_label)
ax1.set_ylabel('Frequency [Hz]')

ax1.set_ylim([0,ylim_]) # new
ax1.set_title('{}'.format(title),y=1.08)

if path != None:
fig1.savefig(path)
path_csv = path.replace(".png",".csv")
#freqs_txt = np.array2string(freqs, precision=2, separator=',',)
data_csv = np.append(freqs[np.newaxis].T,Pxx, axis=1)
bins_txt = "freq (Hz)/bins (s)," + np.array2string(bins, precision=2, separator=',',).replace("[","").replace("]","")
np.savetxt(path_csv, data_csv,header=bins_txt, delimiter=",")


def chromagram_from_spectrogram(Pxx,fs,n_fft,n_chroma=24,norm=True):

# Author: Daniel Macdonald
# Calculate chroma filterbank
chromafb = chroma_filterbank(
sr=fs,
Expand All @@ -639,6 +509,7 @@ def chromagram_from_spectrogram(Pxx,fs,n_fft,n_chroma=24,norm=True):
return chroma

def calc_chroma_entropy(chroma,n_chroma):
# Author: Daniel Macdonald
chroma_entropy = -np.sum(chroma * np.log(chroma), axis=0) / np.log(n_chroma)
#print(np.sum(chroma,axis=0))
#chroma_entropy = entropy(normed_power) / np.log(n_chroma)
Expand Down Expand Up @@ -672,6 +543,7 @@ def plot_chromagram(fig1,ax1,bins,chroma,title=None,path=None,convert_a=1.0,conv


def get_sampling_constants(df,start_t,end_t):
# Author: Daniel Macdonald
'''
T = period, in seconds,
nsamples = samples per cycle
Expand All @@ -682,23 +554,8 @@ def get_sampling_constants(df,start_t,end_t):
fs = nsamples/T
return T, nsamples, fs

#def sample_polydata(df, polydata_geo, num_points=2000):
# polydata_masked = vtk_mask_points(polydata_geo, num_points)
# #print(polydata_masked)
# '''
# Find the indices of the df that we'll use, filter df
# '''
# sampled_points_og_ids = get_original_ids(polydata_masked, polydata_geo)
# df_sampled = df.filter(sampled_points_og_ids.astype(int), axis = 0)
# '''
# Reset index here allows us to easily search the masked polydata
# '''
# df_sampled.reset_index(level=0, inplace=True)
# df_sampled.index.names = ['NewIds']
# df_sampled.drop('Ids', 1, inplace = True)
# return df_sampled, polydata_masked

# The following functions are taken from librosa to generate chroma filterbank ############################################################

# The remaining functions are taken from librosa to generate chroma filterbank

def octs_to_hz(octs, tuning=0.0, bins_per_octave=12):
"""Convert octaves numbers to frequencies.
Expand Down
Loading
Loading