initial commit

README.md

@@ -0,0 +1,8 @@
+This code package implements the density integration methodology for BOS outlined in:
+
+Rajendran, L. K., Zhang, J., Bane, S., & Vlachos, P. (2020). Uncertainty-based weighted least squares density integration for background-oriented schlieren. Experiments in Fluids. 
+
+Please cite the above paper if you use this code package for your work.
+
+sample_script.py is a sample python script that loads the sample dataset in 'sample-data.mat' and calls the density integration + uncertainty quantification function to perform the calculations.
+It saves the result to 'sample-result.mat' and a figure to 'sample-result.png'

loadmat_functions.py

@@ -0,0 +1,35 @@
+# The following functions convert an object to a struct so that it can be saved to a mat file
+import scipy.io as sio
+
+def loadmat(filename):
+    '''
+	this function should be called instead of direct spio.loadmat
+	as it cures the problem of not properly recovering python dictionaries
+	from mat files. It calls the function check keys to cure all entries
+	which are still mat-objects
+	'''
+    data = sio.loadmat(filename, struct_as_record=False, squeeze_me=True)
+    return _check_keys(data)
+
+def _check_keys(dict):
+    '''
+	checks if entries in dictionary are mat-objects. If yes
+	todict is called to change them to nested dictionaries
+	'''
+    for key in dict:
+        if isinstance(dict[key], sio.matlab.mio5_params.mat_struct):
+            dict[key] = _todict(dict[key])
+    return dict
+
+def _todict(matobj):
+    '''
+	A recursive function which constructs from matobjects nested dictionaries
+	'''
+    dict = {}
+    for strg in matobj._fieldnames:
+        elem = matobj.__dict__[strg]
+        if isinstance(elem, sio.matlab.mio5_params.mat_struct):
+            dict[strg] = _todict(elem)
+        else:
+            dict[strg] = elem
+    return dict

sample_script.py

@@ -0,0 +1,113 @@
+#!/usr/bin/env python2
+# -*- coding: utf-8 -*-
+
+import numpy as np
+import matplotlib
+# matplotlib.use('Agg')
+import matplotlib.pyplot as plt
+import sys
+import os
+from scipy.io import savemat
+from scipy.io import loadmat
+import timeit
+
+from DensityIntegrationUncertaintyQuantification import Density_integration_Poisson_uncertainty
+from DensityIntegrationUncertaintyQuantification import Density_integration_WLS_uncertainty
+from DensityIntegrationUncertaintyQuantification import Density_integration_WLS_uncertainty_weighted_average
+
+# Load the data: 
+# all loaded variables is in standard physical units: X, Y: [m], rho_x, rho_y: [kg/m^4] 
+data = loadmat('sample-data.mat', squeeze_me=True)
+
+# set the density uncertainty at the boundary points [kg/m^3]
+sigma_rho_dirichlet = 0.01
+
+# calculate density and uncertainty (Poisson)
+rho_poisson, sigma_rho_poisson = Density_integration_Poisson_uncertainty(data['X'], data['Y'], data['mask'],
+                                                         data['rho_x'], data['rho_y'],
+                                                         data['dirichlet_label'], data['rho_dirichlet'],
+                                                         uncertainty_quantification=True,
+                                                         sigma_grad_x=data['sigma_rho_x'], sigma_grad_y=data['sigma_rho_y'],
+                                                         sigma_dirichlet=sigma_rho_dirichlet)
+
+# calculate density and uncertainty (WLS)
+rho_wls, sigma_rho_wls = Density_integration_WLS_uncertainty(data['X'], data['Y'], data['mask'],
+                                                         data['rho_x'], data['rho_y'],
+                                                         data['dirichlet_label'], data['rho_dirichlet'],
+                                                         uncertainty_quantification=True,
+                                                         sigma_grad_x=data['sigma_rho_x'], sigma_grad_y=data['sigma_rho_y'],
+                                                         sigma_dirichlet=sigma_rho_dirichlet)
+
+# save the results to file
+savemat(file_name='sample-result.mat', mdict={'X': data['X'], 'Y': data['Y'], 'rho_poisson': rho_poisson, 'sigma_rho_poisson': sigma_rho_poisson,
+                                    'rho_wls': rho_wls, 'sigma_rho_wls': sigma_rho_wls}, long_field_names=True)
+
+# Plot the results
+fig1 = plt.figure(1, figsize=(12,8))
+plt.figure(1)
+ax1 = fig1.add_subplot(3,2,1)
+ax2 = fig1.add_subplot(3,2,2)
+ax3 = fig1.add_subplot(3,2,3)
+ax4 = fig1.add_subplot(3,2,4)
+ax5 = fig1.add_subplot(3,2,5)
+ax6 = fig1.add_subplot(3,2,6)
+
+# plot x gradient
+plt.axes(ax1)
+plt.pcolor(data['X'], data['Y'], data['rho_x'], vmin=-60, vmax=60)
+plt.colorbar()
+plt.title('rho_x')
+
+# plot y gradient
+plt.axes(ax2)
+plt.pcolor(data['X'], data['Y'], data['rho_y'], vmin=-60, vmax=60)
+plt.colorbar()
+plt.title('rho_y')
+
+# plot density (poisson)
+plt.axes(ax3)
+plt.pcolor(data['X'], data['Y'], rho_poisson, vmin=1.2, vmax=1.5)
+plt.colorbar()
+plt.title('rho, Poisson')
+
+# plot density uncertainty (poisson)
+plt.axes(ax4)
+plt.pcolor(data['X'], data['Y'], sigma_rho_poisson, vmin=0.0, vmax=0.01)
+plt.colorbar()
+plt.title('sigma rho, Poisson')
+
+# plot density (wls)
+plt.axes(ax5)
+plt.pcolor(data['X'], data['Y'], rho_wls, vmin=1.2, vmax=1.5)
+plt.colorbar()
+plt.title('rho, WLS')
+
+# plot density uncertainty (wls)
+plt.axes(ax6)
+plt.pcolor(data['X'], data['Y'], sigma_rho_wls, vmin=0.0, vmax=0.01)
+plt.colorbar()
+plt.title('sigma rho, WLS')
+
+plt.tight_layout()
+
+# save plot to file
+plt.savefig('sample-result.png')
+plt.close()
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+