Merge pull request #51 from dbruneau-mie/master

Add domain-based spectrogram sampling

postprocessing_h5py/create_spectrogram_from_components.py

@@ -69,14 +69,18 @@ def create_spectrogram_composite(case_path, dvp, n_samples, thresh_val, max_plot
     fig2, ax2_1 = plt.subplots()
     fig2.set_size_inches(7.5, 5) #fig1.set_size_inches(10, 7)
     title = "threshold Pxx = {}".format(thresh_val)
-    path_to_fig = os.path.join(imageFolder, fullname)
-    spec.plot_spectrogram(fig2,ax2_1,bins,freqs,Pxx,ylim,title=title,path=path_to_fig,x_label="Time (s)",color_range=[thresh_val,max_plot])
-    fig2.savefig(path_to_fig)
+    path_to_spec = os.path.join(imageFolder, fullname)
+    spec.plot_spectrogram(fig2,ax2_1,bins,freqs,Pxx,ylim,title=title,x_label="Time (s)",color_range=[thresh_val,max_plot])
+    fig2.savefig(path_to_spec)
+    data_csv = np.append(freqs[np.newaxis].T,Pxx, axis=1)
+    path_csv = path_to_spec.replace(".png",".csv")
+    np.savetxt(path_csv, data_csv,header=bins_txt, delimiter=",")
+
 
 
 if __name__ == '__main__':
     # Load in case-specific parameters
-    case_path, mesh_name, save_deg, stride,  start_t, end_t, lowcut, ylim, r_sphere, x_sphere, y_sphere, z_sphere, dvp, _, _, interface_only, sampling_method, component, _, point_id = spec.read_command_line_spec()
+    case_path, mesh_name, save_deg, stride,  start_t, end_t, lowcut, ylim, _,_,_, r_sphere, x_sphere, y_sphere, z_sphere, dvp, _, _, interface_only, sampling_method, component, _, point_id = spec.read_command_line_spec()
 
     # Read fixed spectrogram parameters from config file
     config_file = os.path.join(os.path.dirname(os.path.realpath(__file__)),"Spectrogram.config")

postprocessing_h5py/create_spectrograms_chromagrams.py

@@ -105,7 +105,7 @@ def create_spectrogram_composite(case_name, dvp, df, start_t, end_t,
                                                                                         thresh_method="old")
     bins = bins+start_t # Need to shift bins so that spectrogram timing is correct
     spec.plot_spectrogram(fig1,ax2,bins,freqs,Pxx,ylim,color_range=[thresh_val,max_plot])
-    
+
 
     # Chromagram ------------------------------------------------------------
     n_fft = spec.shift_bit_length(int(df.shape[1]/nWindow))*2 
@@ -190,11 +190,37 @@ def create_spectrogram_composite(case_name, dvp, df, start_t, end_t,
     fig2, ax2_1 = plt.subplots()
     fig2.set_size_inches(7.5, 5) #fig1.set_size_inches(10, 7)
     title = "Pxx max = {:.2e}, Pxx min = {:.2e}, threshold Pxx = {}".format(max_val, min_val, lower_thresh)
+    spec.plot_spectrogram(fig2,ax2_1,bins,freqs,Pxx,ylim,title=title,x_label="Time (s)",color_range=[thresh_val,max_plot])
+    # Save data to files (spectrogram, chromagram, SBI)
+
     fullname = dvp+"_"+case_name + '_'+str(nWindow)+'_windows_'+'_'+"thresh"+str(thresh_val)+"_spectrogram"
-    path_to_fig = os.path.join(imageFolder, fullname + '.png')
-    spec.plot_spectrogram(fig2,ax2_1,bins,freqs,Pxx,ylim,title=title,path=path_to_fig,x_label="Time (s)",color_range=[thresh_val,max_plot])
-    fig2.savefig(path_to_fig)
+    path_to_spec = os.path.join(imageFolder, fullname + '.png')
+    fig2.savefig(path_to_spec)
+    path_csv = path_to_spec.replace(".png",".csv")
+    #freqs_txt = np.array2string(freqs, precision=2, separator=',',)
+    data_csv = np.append(freqs[np.newaxis].T,Pxx, axis=1)
+    bins_txt =  np.array2string(bins, max_line_width=10000, precision=2, separator=',',).replace("[","").replace("]","")
+    np.savetxt(path_csv, data_csv,header=bins_txt, delimiter=",")
+
+    # Save data to files (spectrogram, chromagram, SBI)
+    fullname = dvp+"_"+case_name + '_'+str(nWindow)+'_windows_'+'_chromagram'
+    path_to_chroma = os.path.join(imageFolder, fullname + '.png')
+    path_csv = path_to_chroma.replace(".png",".csv")
+    chroma_y = np.linspace(0,1,chroma.shape[0])
+    #freqs_txt = np.array2string(freqs, precision=2, separator=',',)
+    data_csv = np.append(chroma_y[np.newaxis].T,chroma, axis=1)
+    bins_txt =  np.array2string(bins_raw, max_line_width=10000, precision=2, separator=',',).replace("[","").replace("]","")
+    np.savetxt(path_csv, data_csv,header=bins_txt, delimiter=",")
 
+    fullname = dvp+"_"+case_name + '_'+str(nWindow)+'_windows_'+'_SBI'
+    path_to_SBI = os.path.join(imageFolder, fullname + '.png')
+    path_csv = path_to_SBI.replace(".png",".csv")
+    #freqs_txt = np.array2string(freqs, precision=2, separator=',',)
+    print(bins)
+    print(chroma_entropy)
+    data_csv = np.array([bins,chroma_entropy]).T
+    np.savetxt(path_csv, data_csv,header="t (s), SBI", delimiter=",")
+    #bins_raw,chroma
 
 def sonify_point(case_name, dvp, df, start_t, end_t, overlapFrac, lowcut, imageFolder):
 
@@ -220,7 +246,7 @@ def sonify_point(case_name, dvp, df, start_t, end_t, overlapFrac, lowcut, imageF
 
 if __name__ == '__main__':
     # Load in case-specific parameters
-    case_path, mesh_name, save_deg, stride,  start_t, end_t, lowcut, ylim, r_sphere, x_sphere, y_sphere, z_sphere, dvp, _, _, interface_only, sampling_method, component, _, point_id = spec.read_command_line_spec()
+    case_path, mesh_name, save_deg, stride,  start_t, end_t, lowcut, ylim, sampling_region, fluid_sampling_domain_ID, solid_sampling_domain_ID, r_sphere, x_sphere, y_sphere, z_sphere, dvp, _, _, interface_only, sampling_method, component, _, point_id = spec.read_command_line_spec()
 
     # Read fixed spectrogram parameters from config file
     config_file = os.path.join(os.path.dirname(os.path.realpath(__file__)),"Spectrogram.config")
@@ -235,6 +261,9 @@ def sonify_point(case_name, dvp, df, start_t, end_t, overlapFrac, lowcut, imageF
                                                                                                            end_t, 
                                                                                                            n_samples, 
                                                                                                            ylim, 
+                                                                                                           sampling_region,
+                                                                                                           fluid_sampling_domain_ID,
+                                                                                                           solid_sampling_domain_ID,
                                                                                                            r_sphere, 
                                                                                                            x_sphere, 
                                                                                                            y_sphere, 

postprocessing_h5py/determine_modal_frequencies.py

@@ -33,8 +33,10 @@
 
 def determine_modal_frequencies(case_path, dvp, n_samples, thresh_val, max_plot):
 
-    idy_time=12 # time index for determining peak freqs
-    thresh_peak=-40 # Power threshold in db for determining peaks
+
+
+    #idy_time=12 # time index for determining peak freqs
+    thresh_peak=-41 # Power threshold in db for determining peaks
     peak_idx_min_gap=4 # minimum spacing between peaks, in terms of frequency index (4 for surgical cases)
     # Get viz path
     visualization_path = postprocessing_common_h5py.get_visualization_path(case_path)    
@@ -45,6 +47,7 @@ def determine_modal_frequencies(case_path, dvp, n_samples, thresh_val, max_plot)
 
     # create spec name
     fullname=os.path.basename(x_csv_files[0]).replace("d_combined_","d_overlayLines_").replace(".csv",".png")
+    fullname=os.path.basename(x_csv_files[0]).replace("d_combined_","d_overlayLines_").replace(".csv",'')
 
     # get bins from x csv file header
     bins_file = open(x_csv_files[0], "r")
@@ -61,50 +64,60 @@ def determine_modal_frequencies(case_path, dvp, n_samples, thresh_val, max_plot)
 
     # Frequencies are the first column of the data
     freqs=csv_x_data[:,0]
-
-    # Average the component10
-    Pxx = csv_x_data[:,1:]
-    peaks_indices = find_peaks_cwt(Pxx[:,idy_time], np.arange(3,6) ,noise_perc=40) # gap_thresh=10
-    print(peaks_indices)
-    #print(Pxx[peaks_indices,idy_time])
-    #print(peaks_indices)
-    peak_idx_prev = 0
-    peak_idx_filtered = []
-    for peak_idx in peaks_indices:
-        if Pxx[peak_idx,idy_time] > thresh_peak:
-            if (peak_idx - peak_idx_prev >= peak_idx_min_gap):
-                peak_idx_filtered.append(peak_idx)
-            else:
-                if Pxx[peak_idx,idy_time] >  Pxx[peak_idx_prev,idy_time]:
-                    del peak_idx_filtered[-1]
+    Pxx = csv_x_data[:,1:]    
+    n_ts = len(Pxx[0,:])
+
+    print(n_ts)
+    for idy_time in range(n_ts):
+        # Average the component10
+        peaks_indices = find_peaks_cwt(Pxx[:,idy_time], np.arange(3,6) ,noise_perc=40) # gap_thresh=10
+        #print(peaks_indices)
+        #print(Pxx[peaks_indices,idy_time])
+        #print(peaks_indices)
+        peak_idx_prev = 0
+        peak_idx_filtered = []
+        for peak_idx in peaks_indices:
+            if Pxx[peak_idx,idy_time] > thresh_peak:
+                if (peak_idx - peak_idx_prev >= peak_idx_min_gap):
                     peak_idx_filtered.append(peak_idx)
-                #else:
-                #    peak_idx_filtered.append(peak_idx_prev)
-
-        peak_idx_prev = peak_idx
-
-    modal_freqs = freqs[peak_idx_filtered]
-
-    # create separate spectrogram figure
-    fig2, ax2_1 = plt.subplots()
-    fig2.set_size_inches(7.5, 5) #fig1.set_size_inches(10, 7)
-    title = "threshold Pxx = {}".format(thresh_val)
-    path_to_fig = os.path.join(imageFolder, fullname)
-
-    spec.plot_spectrogram(fig2,ax2_1,bins,freqs,Pxx,ylim,title=title,path=path_to_fig,x_label="Time (s)",color_range=[thresh_val,max_plot])
-    for i in range(0,4):
-        plt.axhline(y = modal_freqs[i], color = 'r', linestyle = '-', label = "{:.2f} Hz".format(modal_freqs[i]))
+                else:
+                    if Pxx[peak_idx,idy_time] >  Pxx[peak_idx_prev,idy_time]:
+                        del peak_idx_filtered[-1]
+                        peak_idx_filtered.append(peak_idx)
+                    #else:
+                    #    peak_idx_filtered.append(peak_idx_prev)
+
+            peak_idx_prev = peak_idx
+
+        modal_freqs = freqs[peak_idx_filtered]
+
+        # create separate spectrogram figure
+        fig2, ax2_1 = plt.subplots()
+        fig2.set_size_inches(7.5, 5) #fig1.set_size_inches(10, 7)
+        title = "threshold Pxx = {}".format(thresh_val)
+        file_name_ts = fullname + "_" + str(idy_time) + ".png"
+        path_to_fig = os.path.join(imageFolder, file_name_ts)
+
+
+        spec.plot_spectrogram(fig2,ax2_1,bins,freqs,Pxx,ylim,title=title,x_label="Time (s)",color_range=[thresh_val,max_plot])
+        for i in range(0,4):
+            try:
+                plt.axhline(y = modal_freqs[i], color = 'r', linestyle = '-', label = "{:.2f} Hz".format(modal_freqs[i]))
+                print(os.path.basename(case_path) + " Mode {}: {:.2f} Hz".format(modal_freqs[i],idy))
+
+            except:
+                print("final mode found")
+
+        #print(os.path.basename(case_path) + " Mode 0: {:.2f} Hz, Mode 1: {:.2f} Hz, Mode 2: {:.2f} Hz, Mode 3: {:.2f} Hz, ts {}".format(modal_freqs[0],modal_freqs[1],modal_freqs[2],modal_freqs[3],idy))
+        plt.axvline(x = bins[idy_time], color = 'b', linestyle = '-')
+        plt.legend()
+        fig2.savefig(path_to_fig)
 
-    print(os.path.basename(case_path) + " Mode 0: {:.2f} Hz, Mode 1: {:.2f} Hz, Mode 2: {:.2f} Hz, Mode 3: {:.2f} Hz,".format(modal_freqs[0],modal_freqs[1],modal_freqs[2],modal_freqs[3]))
-
-    plt.axvline(x = bins[idy_time], color = 'b', linestyle = '-')
-    plt.legend()
-    fig2.savefig(path_to_fig)
 
 
 if __name__ == '__main__':
     # Load in case-specific parameters
-    case_path, mesh_name, save_deg, stride,  start_t, end_t, lowcut, ylim, r_sphere, x_sphere, y_sphere, z_sphere, dvp, _, _, interface_only, sampling_method, component, _, point_id = spec.read_command_line_spec()
+    case_path, mesh_name, save_deg, stride,  start_t, end_t, lowcut, ylim, _,_,_, r_sphere, x_sphere, y_sphere, z_sphere, dvp, _, _, interface_only, sampling_method, component, _, point_id = spec.read_command_line_spec()
 
     # Read fixed spectrogram parameters from config file
     config_file = os.path.join(os.path.dirname(os.path.realpath(__file__)),"Spectrogram.config")

postprocessing_h5py/postprocessing_common_h5py.py

@@ -100,6 +100,35 @@ def get_domain_ids(meshFile):
     allIDs = np.unique(allTopology) 
     return fluidIDs, wallIDs, allIDs
 
+def get_domain_ids_specified_region(meshFile,fluid_sampling_domain_ID,solid_sampling_domain_ID):
+    # This function obtains the topology for the fluid, solid, and all elements of the input mesh
+    vectorData = h5py.File(meshFile,"r")
+    domainsLoc = 'domains/values'
+    domains = vectorData[domainsLoc][:] # Open domain array
+    id_wall = (domains==solid_sampling_domain_ID).nonzero() # domain = 2 is the solid
+    id_fluid = (domains==fluid_sampling_domain_ID).nonzero() # domain = 1 is the fluid
+
+    topologyLoc = 'domains/topology'
+    allTopology = vectorData[topologyLoc][:,:] 
+    wallTopology=allTopology[id_wall,:] 
+    fluidTopology=allTopology[id_fluid,:]
+
+    wallIDs = np.unique(wallTopology) # find the unique node ids in the wall topology, sorted in ascending order
+    fluidIDs = np.unique(fluidTopology) # find the unique node ids in the fluid topology, sorted in ascending order
+    allIDs = np.unique(allTopology) 
+
+    return fluidIDs, wallIDs, allIDs
+
+def get_domain_ids(meshFile):
+    # This function obtains a list of the node IDs for the fluid, solid, and all elements of the input mesh
+
+    # Get topology of fluid, solid and whole mesh
+    fluidTopology, wallTopology, allTopology = get_domain_topology(meshFile)
+    wallIDs = np.unique(wallTopology) # find the unique node ids in the wall topology, sorted in ascending order
+    fluidIDs = np.unique(fluidTopology) # find the unique node ids in the fluid topology, sorted in ascending order
+    allIDs = np.unique(allTopology) 
+    return fluidIDs, wallIDs, allIDs
+
 def get_interface_ids(meshFile):
     fluidIDs, wallIDs, allIDs = get_domain_ids(meshFile)
     interfaceIDs_set = set(fluidIDs) - (set(fluidIDs) - set(wallIDs))