Add files via upload

Turbulence_Sim_v1_python/Example_1_full.py

@@ -0,0 +1,47 @@
+'''
+Demo of Multi-Aperture Turbulence Simulation.
+
+Nicholas Chimitt and Stanley H. Chan "Simulating anisoplanatic turbulence
+by sampling intermodal and spatially correlated Zernike coefficients," Optical Engineering 59(8), Aug. 2020
+
+ArXiv:  https://arxiv.org/abs/2004.11210
+
+Nicholas Chimitt and Stanley Chan
+Copyright 2021
+Purdue University, West Lafayette, In, USA.
+'''
+
+from matplotlib import pyplot as plt
+import numpy as np
+import TurbSim_v1_main as util
+
+
+def rgb2gray(rgb):
+    return np.dot(rgb[...,:3], [0.2989, 0.5870, 0.1140])
+
+
+img = rgb2gray(plt.imread('chart.png'))
+N = 225 # size of the image -- assumed to be square (pixels)
+D = 0.2 # length of aperture diameter (meters)
+L = 7000 # length of propagation (meters)
+r0 = 0.05 # the Fried parameter r0. The value of D/r0 is critically important! (See associated paper)
+wvl = 0.525e-6 # the mean wavelength -- typically somewhere suitably in the middle of the spectrum will be sufficient
+obj_size = 2.06*1 # the size of the object in the object plane (meters). Can be different the Nyquist sampling, scaling
+                  # will be done automatically.
+
+param_obj = util.p_obj(N, D, L, r0, wvl, obj_size) # generating the parameter object, some other things are computed within this
+                                                   # function, see the def for details
+print(param_obj['N'] * param_obj['delta0'], param_obj['smax'], param_obj['scaling'])
+#print('hi')
+#print(param_obj['spacing'])
+#print('hi')
+S = util.gen_PSD(param_obj) # finding the PSD, see the def for details
+param_obj['S'] = S # appending the PSD to the parameter object for convenience
+
+for i in range(100):
+    img_tilt, _ = util.genTiltImg(img, param_obj) # generating the tilt-only image
+    img_blur = util.genBlurImage(param_obj, img_tilt)
+    plt.imshow(img_tilt,cmap='gray',vmin=0,vmax=1)
+    plt.show()
+    plt.imshow(img_blur, cmap='gray', vmin=0, vmax=1)
+    plt.show()

Turbulence_Sim_v1_python/Example_2_tiltonly.py

@@ -0,0 +1,39 @@
+'''
+Demo of Multi-Aperture Turbulence Simulation (Tilt-only).
+
+Nicholas Chimitt and Stanley H. Chan "Simulating anisoplanatic turbulence
+by sampling intermodal and spatially correlated Zernike coefficients," Optical Engineering 59(8), Aug. 2020
+
+ArXiv:  https://arxiv.org/abs/2004.11210
+
+Nicholas Chimitt and Stanley Chan
+Copyright 2021
+Purdue University, West Lafayette, In, USA.
+'''
+
+from matplotlib import pyplot as plt
+import numpy as np
+import TurbSim_v1_main as util
+
+
+def rgb2gray(rgb):
+    return np.dot(rgb[...,:3], [0.2989, 0.5870, 0.1140])
+
+
+img = rgb2gray(plt.imread('chart.png'))
+N = 225 # size of the image -- assumed to be square (pixels)
+D = 0.1 # length of aperture diameter (meters)
+L = 7000 # length of propagation (meters)
+r0 = D/2 # the Fried parameter r0. Most importantly, the denominator is the desired D/r0 ratio
+wvl = 0.5e-6 # the mean wavelength -- typically somewhere suitably in the middle of the spectrum will be sufficient
+
+param_obj = util.p_obj(N, D, L, r0, wvl) # generating the parameter object, some other things are computed within this
+                                         # function, see the def for details
+S = util.gen_PSD(param_obj) # finding the PSD, see the def for details
+param_obj['S'] = S # appending the PSD to the parameter object for convenience
+
+for i in range(100):
+    img_, _ = util.genTiltImg(img, param_obj) # generating the tilt-only image
+
+    plt.imshow(img_,cmap='gray',vmin=0,vmax=1)
+    plt.show()

Turbulence_Sim_v1_python/Integrals_Spatial_Corr.py

@@ -0,0 +1,94 @@
+'''
+Spatial Correlation Functions for Tilts
+
+Zhiyuan Mao and Stanley Chan
+Copyright 2021
+Purdue University, West Lafayette, In, USA.
+'''
+
+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/Motion_Compensate.py

@@ -0,0 +1,39 @@
+'''
+Warping Function for Turbulence Simulator
+
+Python version
+of original code by Stanley Chan
+
+Zhiyuan Mao and Stanley Chan
+Copyright 2021
+Purdue University, West Lafayette, In, USA.
+'''
+
+import numpy as np
+from skimage.transform import resize
+
+
+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