Initial Upload

Initial Upload. Can be made faster with some precomputation, but speed is acceptable relative to reported results.

Turbulence_Sim_v1_python/integrals_spatial_corr.py

@@ -0,0 +1,85 @@
+import numpy as np
+import scipy.integrate as integrate
+from scipy.special import jv
+import os
+from math import gamma
+
+def genTiltPSD(s_max, spacing, N, **kwargs):
+	D = 0.1
+	D_r0 = 20
+	wavelength = kwargs.get('wavelength', 500 * (10 ** (-9)))
+	L = 7000
+	delta0 = L * wavelength / (2 * D)
+	delta0 *= 1  # Adjusting factor
+	# Constants in [Channan, 1992]
+	c1 = 2 * ((24 / 5) * gamma(6 / 5)) ** (5 / 6)
+	c2 = 4 * c1 / np.pi * (gamma(11 / 6)) ** 2
+	s_arr = np.arange(0, s_max, spacing)
+	I0_arr = np.float32(s_arr * 0)
+	I2_arr = np.float32(s_arr * 0)
+	for i in range(len(s_arr)):
+		I0_arr[i] = I0(s_arr[i])
+		I2_arr[i] = I2(s_arr[i])
+	i, j = np.int32(N / 2), np.int32(N / 2)
+	[x, y] = np.meshgrid(np.arange(1, N + 0.01, 1), np.arange(1, N + 0.01, 1))
+	s = np.sqrt((x - i) ** 2 + (y - j) ** 2)
+	s *= spacing
+	C0 = (In_m(s, spacing*100, I0_arr) + In_m(s, spacing*100, I2_arr)) / I0(0)
+	C0[i, j] = 1
+	C0_scaled = C0 * I0(0) * c2 * ((D_r0) ** (5 / 3)) / (2 ** (5 / 3)) * (
+				(2 * wavelength / (np.pi * D)) ** 2) * 2 * np.pi
+	Cfft = np.fft.fft2(C0_scaled)
+	S_half = np.sqrt(Cfft)
+	S_half_max = np.max(np.max(np.abs(S_half)))
+	S_half[np.abs(S_half) < 0.0001 * S_half_max] = 0
+
+	return S_half
+
+
+def I0(s):
+	# z = np.linspace(1e-6, 1e3, 1e5)
+	# f_z = (z**(-14/3))*jv(0,2*s*z)*(jv(2,z)**2)
+	# I0_s = np.trapz(f_z, z)
+
+	I0_s, _ = integrate.quad( lambda z: (z**(-14/3))*jv(0,2*s*z)*(jv(2,z)**2), 1e-4, np.inf, limit = 100000)
+	# print('I0: ',I0_s)
+	return I0_s
+
+
+def I2(s):
+	# z = np.linspace(1e-6, 1e3, 1e5)
+	# f_z = (z**(-14/3))*jv(2,2*s*z)*(jv(2,z)**2)
+	# I2_s = np.trapz(f_z, z)
+
+	I2_s, _ = integrate.quad( lambda z: (z**(-14/3))*jv(2,2*s*z)*(jv(2,z)**2), 1e-4, np.inf, limit = 100000)
+	# print('I2: ',I2_s)
+	return I2_s
+
+
+def save_In(s_max, spacing):
+	s_arr = np.arange(0,s_max,spacing)
+	I0_arr = np.float32(s_arr*0)
+	I2_arr = np.float32(s_arr*0)
+	for i in range(len(s_arr)):
+		I0_arr[i] = I0(s_arr[i])
+		I2_arr[i] = I2(s_arr[i])
+	os.makedirs('temp/In_arr', exist_ok=True)
+
+	np.save('temp/In_arr/I0_%d_%d.npy'%(s_max,spacing*1000),I0_arr)
+	np.save('temp/In_arr/I2_%d_%d.npy'%(s_max,spacing*1000),I2_arr)
+
+
+def In_m(s, spacing, In_arr):
+	idx = np.int32(np.floor(s.flatten()/spacing))
+	M,N = np.shape(s)[0], np.shape(s)[1]
+	In = np.reshape(np.take(In_arr, idx), [M,N])
+
+	return In
+
+
+def In_arr(s, spacing, In_arr):
+	idx = np.int32(np.floor(s/spacing))
+	In = np.take(In_arr, idx)
+
+	return In
+

Turbulence_Sim_v1_python/main.py

@@ -0,0 +1,25 @@
+from matplotlib import pyplot as plt
+import numpy as np
+import utilities as util
+
+
+def rgb2gray(rgb):
+    return np.dot(rgb[...,:3], [0.2989, 0.5870, 0.1140])
+
+img = rgb2gray(plt.imread('chart.png'))
+# generate S_half for tilt
+N = 225
+D = 0.2
+L = 7000
+r0 = D/2 #Denominator is the desired D/r0 ratio
+wvl = 0.5e-6
+
+
+param_obj = util.p_obj(N, D, L, r0, wvl)
+S = util.gen_PSD(param_obj)
+param_obj['S'] = S
+
+for i in range(100):
+    img_ = util.genTiltImg(img, param_obj)
+    plt.imshow(img_,cmap='gray',vmin=0,vmax=1)
+    plt.show()

Turbulence_Sim_v1_python/motion_compensate.py

@@ -0,0 +1,31 @@
+import numpy as np
+from skimage.transform import resize
+import matplotlib.pyplot as plt
+
+# Python version 
+# of original code by Prof. Stanley Chan
+
+def motion_compensate(img, Mvx, Mvy, pel):
+	m, n =  np.shape(img)[0], np.shape(img)[1]
+	img = resize(img, (np.int32(m/pel), np.int32(n/pel)), mode = 'reflect' )
+	Blocksize = np.floor(np.shape(img)[0]/np.shape(Mvx)[0])
+	m, n =  np.shape(img)[0], np.shape(img)[1]
+	M, N =  np.int32(np.ceil(m/Blocksize)*Blocksize), np.int32(np.ceil(n/Blocksize)*Blocksize)
+
+	f = img[0:M, 0:N]
+
+
+	Mvxmap = resize(Mvy, (N,M))
+	Mvymap = resize(Mvx, (N,M))
+
+
+	xgrid, ygrid = np.meshgrid(np.arange(0,N-0.99), np.arange(0,M-0.99))
+	X = np.clip(xgrid+np.round(Mvxmap/pel),0,N-1)
+	Y = np.clip(ygrid+np.round(Mvymap/pel),0,M-1)
+
+	idx = np.int32(Y.flatten()*N + X.flatten())
+	f_vec = f.flatten()
+	g = np.reshape(f_vec[idx],[N,M])
+
+	g = resize(g, (np.shape(g)[0]*pel,np.shape(g)[1]*pel))
+	return g