Clean-up - only fully validated code kept and all others removed.
This commit is contained in:
parent
db398fe775
commit
5b1428e91a
@ -1,233 +0,0 @@
|
|||||||
"""
|
|
||||||
@author: Adam Newton Wright, https://github.com/adamnewtonwright/GaussianBeamPropagation
|
|
||||||
"""
|
|
||||||
|
|
||||||
import numpy as np
|
|
||||||
import matplotlib.pyplot as plt
|
|
||||||
import sympy as sym
|
|
||||||
from sympy import oo
|
|
||||||
|
|
||||||
|
|
||||||
# Input Ray parameter, i.e. height and angle
|
|
||||||
def ray(y,theta):
|
|
||||||
'''
|
|
||||||
Parameters
|
|
||||||
----------
|
|
||||||
y : float or integer or sympy symbol in meters
|
|
||||||
The vertical height of a ray.
|
|
||||||
theta : float or integer in radians
|
|
||||||
The angle of divergence of the ray.
|
|
||||||
|
|
||||||
Returns
|
|
||||||
-------
|
|
||||||
mat : 2x1 matrix
|
|
||||||
[
|
|
||||||
[y],
|
|
||||||
[teta]
|
|
||||||
]
|
|
||||||
|
|
||||||
'''
|
|
||||||
|
|
||||||
mat = np.array([[y],[theta]])
|
|
||||||
return mat
|
|
||||||
|
|
||||||
# Ray Transfer Matrix for ideal lens with focal length f
|
|
||||||
def lens(f):
|
|
||||||
'''
|
|
||||||
Parameters
|
|
||||||
----------
|
|
||||||
f : float or integer or sympy symbol in meters
|
|
||||||
Thin lens focal length in meters
|
|
||||||
|
|
||||||
Returns
|
|
||||||
-------
|
|
||||||
mat : 2x2 matrix
|
|
||||||
[
|
|
||||||
[ 1, 0],
|
|
||||||
[-1/f, 1]
|
|
||||||
]
|
|
||||||
|
|
||||||
'''
|
|
||||||
|
|
||||||
mat = np.array([[1,0], [-1/f, 1]])
|
|
||||||
return mat
|
|
||||||
|
|
||||||
# Ray Transfer Matrix for propagation of distance d
|
|
||||||
def prop(d):
|
|
||||||
'''
|
|
||||||
Parameters
|
|
||||||
----------
|
|
||||||
d : float or integer or sympy symbol
|
|
||||||
Distance light is propagating along the z-axis.
|
|
||||||
|
|
||||||
Returns
|
|
||||||
-------
|
|
||||||
mat: 2x2 matrix
|
|
||||||
[
|
|
||||||
[1, d],
|
|
||||||
[0, 1]
|
|
||||||
]
|
|
||||||
|
|
||||||
'''
|
|
||||||
mat = np.array([[1,d], [0,1]])
|
|
||||||
return mat
|
|
||||||
|
|
||||||
# multiplying the matrices together. mat1 is the last matrix the light interacts with
|
|
||||||
def mult(mat1,*argv):
|
|
||||||
'''
|
|
||||||
Parameters
|
|
||||||
----------
|
|
||||||
mat1 : 2x2 ABCD matrix
|
|
||||||
Last matrix light interacts with.
|
|
||||||
*argv : 2x2 ABCD matrices
|
|
||||||
From left to right, the matrices should be entered such that the leftmost matrix interacts
|
|
||||||
with light temporally after the rightmost matrix.
|
|
||||||
|
|
||||||
Returns
|
|
||||||
-------
|
|
||||||
Mat : 2x2 matrix
|
|
||||||
The ABCd matrix describing the whole optical system.
|
|
||||||
|
|
||||||
'''
|
|
||||||
|
|
||||||
Mat = mat1
|
|
||||||
for arg in argv:
|
|
||||||
Mat = np.dot(Mat, arg)
|
|
||||||
return Mat
|
|
||||||
|
|
||||||
# Adding Gaussian beam parameters
|
|
||||||
def Zr(wo, lam):
|
|
||||||
'''
|
|
||||||
Parameters
|
|
||||||
----------
|
|
||||||
wo : float, integer, or symbol
|
|
||||||
Beam waist radius in meters.
|
|
||||||
lam : float, integer, or symbol
|
|
||||||
Wavelength of light in meters.
|
|
||||||
|
|
||||||
Returns
|
|
||||||
-------
|
|
||||||
zr : float, int, symbols
|
|
||||||
Rayleigh range for given beam waist and wavelength.
|
|
||||||
|
|
||||||
'''
|
|
||||||
|
|
||||||
zr = np.pi * wo**2 / lam
|
|
||||||
return zr
|
|
||||||
|
|
||||||
def W0(zr, lam):
|
|
||||||
'''
|
|
||||||
Parameters
|
|
||||||
----------
|
|
||||||
zr : float, integer, symbol
|
|
||||||
Rayleigh range in meters
|
|
||||||
lam : float, integer, symbol
|
|
||||||
Wavelength of light in meters
|
|
||||||
|
|
||||||
Returns
|
|
||||||
-------
|
|
||||||
w0 : float, integer, symbol
|
|
||||||
Beam waist radius in meters
|
|
||||||
|
|
||||||
'''
|
|
||||||
|
|
||||||
w0 = np.sqrt(lam * zr / np.pi)
|
|
||||||
return w0
|
|
||||||
|
|
||||||
# Remember, there should be an i in front of zr
|
|
||||||
# but this complicates the calculations, so we usually just let z = 0
|
|
||||||
# and don't explicitly deal with the i, but still do the math accordingly
|
|
||||||
#def q0_func(z,zr):
|
|
||||||
# qz = z + zr
|
|
||||||
# return qz
|
|
||||||
|
|
||||||
def q1_func(z, w0, lam, mat):
|
|
||||||
'''
|
|
||||||
Parameters
|
|
||||||
----------
|
|
||||||
z : float, int, symbol
|
|
||||||
Position of the beam waist in meters.
|
|
||||||
w0 : float, int, symbol
|
|
||||||
Radial waist size in meters (of the embedded Gaussian, i.e. W0/M).
|
|
||||||
lam : float, int, symbol
|
|
||||||
Wavelength of light in meters.
|
|
||||||
mat : float, int, symbol
|
|
||||||
The ABCD 2x2 matrix describing the optical system.
|
|
||||||
|
|
||||||
Returns
|
|
||||||
-------
|
|
||||||
z: float, int, symbol
|
|
||||||
Position of the beam waist after the optical system
|
|
||||||
zr: float, int, symbol
|
|
||||||
Rayleigh range of the beam after the optical system
|
|
||||||
'''
|
|
||||||
|
|
||||||
A = mat[0][0]
|
|
||||||
B = mat[0][1]
|
|
||||||
C = mat[1][0]
|
|
||||||
D = mat[1][1]
|
|
||||||
zr = Zr(w0, lam)
|
|
||||||
real = (A*C*(z**2 + zr**2) + z*(A*D + B*C) + B*D) / (C**2*(z**2 + zr**2) + 2*C*D*z + D**2)
|
|
||||||
imag = (zr * (A*D - B*C)) / (C**2*(z**2 + zr**2) + 2*C*D*z + D**2)
|
|
||||||
z = real
|
|
||||||
zr = imag
|
|
||||||
return z, zr
|
|
||||||
|
|
||||||
def q1_inv_func(z, w0, lam, mat):
|
|
||||||
'''
|
|
||||||
Parameters
|
|
||||||
----------
|
|
||||||
z : float, int, symbol
|
|
||||||
Position of the beam waist in meters.
|
|
||||||
w0 : float, int, symbol
|
|
||||||
Radial waist size in meters (of the embedded Gaussian, i.e. W0/M).
|
|
||||||
lam : float, int, symbol
|
|
||||||
Wavelength of light in meters.
|
|
||||||
mat : float, int, symbol
|
|
||||||
The ABCD 2x2 matrix describing the optical system.
|
|
||||||
|
|
||||||
Returns
|
|
||||||
-------
|
|
||||||
R : float, int, symbol
|
|
||||||
Radius of curvature of the wavefront in meters.
|
|
||||||
w : float, int, symbol
|
|
||||||
Radius of the beam in meters.
|
|
||||||
|
|
||||||
'''
|
|
||||||
A = mat[0][0]
|
|
||||||
B = mat[0][1]
|
|
||||||
C = mat[1,0]
|
|
||||||
D = mat[1][1]
|
|
||||||
zr = Zr(w0, lam)
|
|
||||||
real = (A*C*(z**2 + zr**2) + z*(A*D + B*C) + B*D) / (A**2*(z**2 + zr**2) + 2*A*B*z + B**2)
|
|
||||||
imag = -zr * (A*D-B*C) / (A**2 *(z**2 + zr**2) + 2*A*B*z + B**2)
|
|
||||||
R = 1/real
|
|
||||||
w = (-lam / imag / np.pi)**.5
|
|
||||||
return R, w
|
|
||||||
|
|
||||||
def plot(func, var, rang = np.arange(0,3,.01)):
|
|
||||||
'''
|
|
||||||
Parameters
|
|
||||||
----------
|
|
||||||
func : Sympy function of one variable
|
|
||||||
Sympy function defining the beam width after the last optical element.
|
|
||||||
var : sympy variable
|
|
||||||
Variable in func that will be plotted.
|
|
||||||
rang : numpy array
|
|
||||||
Array of the values along the optical axis to be plotted
|
|
||||||
|
|
||||||
Returns
|
|
||||||
-------
|
|
||||||
plot : matplotlib graph
|
|
||||||
Graph of the beam width of var
|
|
||||||
|
|
||||||
|
|
||||||
'''
|
|
||||||
func = sym.lambdify(var, func)
|
|
||||||
plt.figure()
|
|
||||||
plt.plot(rang, func(rang), color = 'b')
|
|
||||||
plt.plot(rang, -func(rang), color = 'b')
|
|
||||||
plt.grid()
|
|
||||||
plt.xlabel('Optic Axis (m)')
|
|
||||||
plt.ylabel('Beam size (m)')
|
|
||||||
plt.show()
|
|
||||||
File diff suppressed because one or more lines are too long
@ -1,70 +0,0 @@
|
|||||||
import BeamPropagation as bp # This is the script that handles the propagation
|
|
||||||
import sympy as sym # For Symbolic examples
|
|
||||||
import numpy as np # Handling of lists and for plotting
|
|
||||||
import matplotlib.pyplot as plt # Plotting
|
|
||||||
|
|
||||||
|
|
||||||
"""A Gaussian beam can be defined by it's (radial) waist wo, it's Rayleigh range zR, and the location of its waist zO"""
|
|
||||||
w0 = 5.2E-3 # 1mm beam waist
|
|
||||||
lam = 532E-9 # wavelength of 355 nm (UV)
|
|
||||||
zR = bp.Zr(w0, lam) # Rayleigh range in m
|
|
||||||
z0 = 0 # location of waist in m
|
|
||||||
|
|
||||||
"""Define first 4-f optical system using matrices"""
|
|
||||||
d1, d2, d3, f1, f2 = sym.symbols('d1 d2 d3 f1 f2')
|
|
||||||
M = bp.mult(bp.prop(d3),bp.lens(f2),bp.prop(d2), bp.lens(f1), bp.prop(d1))
|
|
||||||
|
|
||||||
"""Use script to do all the ABCD and q-parameter math, and return the waist and radius of curvature functions"""
|
|
||||||
R, w = bp.q1_inv_func(0, w0, lam, M)
|
|
||||||
|
|
||||||
"""Substitute and extract the required separation between lenses of first 4-f system"""
|
|
||||||
distance_between_dmd_first_lens = 250E-3
|
|
||||||
first_focal_length = 250.9E-3
|
|
||||||
second_focal_length = 50E-3
|
|
||||||
demag = 1/5
|
|
||||||
target_w0 = demag*w0
|
|
||||||
w = w.subs(f1, first_focal_length).subs(f2, second_focal_length).subs(d1, distance_between_dmd_first_lens).subs(d3,0)
|
|
||||||
eq = sym.solve(w - target_w0, d2)[0]
|
|
||||||
distance_between_lens_of_first_4f = list(eq.atoms())[0]
|
|
||||||
print('Distance between lenses of first 4-f system = {} mm'.format(distance_between_lens_of_first_4f * 1E3))
|
|
||||||
|
|
||||||
# Sanity check
|
|
||||||
# expansion_factor = w.subs(d2,distance_between_lens_of_first_4f).subs(d3,0)/ w0
|
|
||||||
# print('beam is w = {:.2f} x w0'.format(expansion_factor))
|
|
||||||
|
|
||||||
# """Plot beam propagation up to 3 m after the first 4-f system"""
|
|
||||||
# M = bp.mult(bp.prop(d3),bp.lens(second_focal_length),bp.prop(distance_between_lens_of_first_4f), bp.lens(first_focal_length), bp.prop(distance_between_dmd_first_lens))
|
|
||||||
# R, w = bp.q1_inv_func(0, w0, lam, M)
|
|
||||||
# bp.plot(w,d3, rang = np.arange(0,0.050,.0005))
|
|
||||||
|
|
||||||
|
|
||||||
"""Define the full optical system of two 4-f setups using matrices"""
|
|
||||||
d1, d2, d3, d4, f1, f2, f3 = sym.symbols('d1 d2 d3 d4 f1 f2 f3')
|
|
||||||
M = bp.mult(bp.prop(d4),bp.lens(f3), bp.prop(d3),bp.lens(f2),bp.prop(d2), bp.lens(f1), bp.prop(d1))
|
|
||||||
|
|
||||||
"""Use script to do all the ABCD and q-parameter math, and return the waist and radius of curvature functions"""
|
|
||||||
R, w = bp.q1_inv_func(0, w0, lam, M)
|
|
||||||
|
|
||||||
# """Find the focal length of lens required after the first 4-f system to have a collimated beam, given a certain separation between the first 4-f system and this lens"""
|
|
||||||
distance_between_4fs = 550E-3
|
|
||||||
R_coll = R.subs(d1,distance_between_dmd_first_lens).subs(d2,distance_between_lens_of_first_4f).subs(d3,distance_between_4fs).subs(d4,0).subs(f1,first_focal_length).subs(f2,second_focal_length)
|
|
||||||
f3_coll = sym.solve(1/R_coll,f3)[0]
|
|
||||||
third_focal_length = list(f3_coll.atoms())[0]
|
|
||||||
print('For a fixed separation between first 4-f and third lens of {:.3f} mm, f3 = {:.3f} mm for a collimated beam'.format(distance_between_4fs* 1E3, third_focal_length * 1E3))
|
|
||||||
|
|
||||||
# # """Plot beam propagation up to 3 m after the first 4-f system"""
|
|
||||||
# M = bp.mult(bp.prop(d4),bp.lens(third_focal_length),bp.prop(distance_between_4fs), bp.lens(second_focal_length), bp.prop(distance_between_lens_of_first_4f), bp.lens(first_focal_length), bp.prop(distance_between_dmd_first_lens))
|
|
||||||
# R, w = bp.q1_inv_func(0, w0, lam, M)
|
|
||||||
# bp.plot(w,d4, rang = np.arange(0,0.050,.0005))
|
|
||||||
|
|
||||||
third_focal_length = 501.8E-3
|
|
||||||
R_coll = R.subs(d1,distance_between_dmd_first_lens).subs(d2,distance_between_lens_of_first_4f).subs(d4,0).subs(f1,first_focal_length).subs(f2,second_focal_length).subs(f3,third_focal_length)
|
|
||||||
d3_coll = sym.solve(1/R_coll,d3)[1]
|
|
||||||
distance_between_4fs = list(d3_coll.atoms())[0]
|
|
||||||
print('For a fixed third focal length of {:.3f} mm, d3 = {:.3f} mm, for a collimated beam'.format(third_focal_length* 1E3, distance_between_4fs * 1E3))
|
|
||||||
|
|
||||||
# """Plot beam propagation up to 3 m after the first 4-f system"""
|
|
||||||
# M = bp.mult(bp.prop(d4),bp.lens(third_focal_length),bp.prop(distance_between_4fs), bp.lens(second_focal_length), bp.prop(distance_between_lens_of_first_4f), bp.lens(first_focal_length), bp.prop(distance_between_dmd_first_lens))
|
|
||||||
# R, w = bp.q1_inv_func(0, w0, lam, M)
|
|
||||||
# bp.plot(w,d4, rang = np.arange(0,0.050,.0005))
|
|
||||||
|
|
||||||
@ -1,67 +0,0 @@
|
|||||||
import numpy as np
|
|
||||||
import pyfits
|
|
||||||
import matplotlib.pyplot as plt
|
|
||||||
import skimage
|
|
||||||
from skimage.feature import blob_dog, blob_doh, blob_log, canny
|
|
||||||
from skimage.color import rgb2gray
|
|
||||||
from skimage.feature import corner_harris, corner_subpix, corner_peaks
|
|
||||||
from skimage.segmentation import slic
|
|
||||||
from skimage.filters import sobel
|
|
||||||
from scipy.signal import convolve2d
|
|
||||||
|
|
||||||
from scipy.ndimage import gaussian_filter
|
|
||||||
from skimage import measure
|
|
||||||
from scipy.optimize import curve_fit
|
|
||||||
import matplotlib.ticker as mtick
|
|
||||||
from scipy.signal import savgol_filter
|
|
||||||
|
|
||||||
import scipy
|
|
||||||
from scipy import signal
|
|
||||||
from scipy.signal import argrelextrema
|
|
||||||
|
|
||||||
|
|
||||||
import cv2
|
|
||||||
|
|
||||||
## this function will get the values along each circle. We start by defining a range of angles, computing the
|
|
||||||
#coordinates and then finding the values at the nearest pixel position.
|
|
||||||
def get_line(star,theta,radius,x_c,y_c):
|
|
||||||
#theta = np.linspace(0,2*np.pi,N_theta)
|
|
||||||
x = x_c + radius*np.cos(theta)
|
|
||||||
y = y_c + radius*np.sin(theta)
|
|
||||||
x = np.round(x)
|
|
||||||
y = np.round(y)
|
|
||||||
x = x.astype(int)
|
|
||||||
y = y.astype(int)
|
|
||||||
I = star[y,x]
|
|
||||||
|
|
||||||
return I,x,y
|
|
||||||
|
|
||||||
## a function to compute the frequecy for a certain radius
|
|
||||||
def get_radius(freq):
|
|
||||||
N_p = 36
|
|
||||||
r = N_p/(2*np.pi*freq)
|
|
||||||
return r
|
|
||||||
|
|
||||||
## a function to compute the radius for a certain frequency
|
|
||||||
def get_freq(radius):
|
|
||||||
N_p = 36
|
|
||||||
freq = N_p/(2*np.pi*radius)
|
|
||||||
return freq
|
|
||||||
|
|
||||||
def sinusoidal(theta,a,b,c):
|
|
||||||
N_p = 36
|
|
||||||
y = a + b*np.sin(N_p*theta) + c*np.cos(N_p*theta)
|
|
||||||
return y
|
|
||||||
|
|
||||||
def fit_sinusoid(I,theta,p0):
|
|
||||||
popt, pcov = curve_fit(sinusoidal,theta,I,p0)
|
|
||||||
a = popt[0]
|
|
||||||
b = popt[1]
|
|
||||||
c = popt[2]
|
|
||||||
modulation = np.sqrt(b**2 + c**2)/a
|
|
||||||
return modulation, popt
|
|
||||||
|
|
||||||
def Gaussian(x,a,x0,sigma):
|
|
||||||
y = a*(1/(np.sqrt(2*np.pi)*sigma))*np.exp(-(x-x0)**2/(2*sigma**2))
|
|
||||||
#y = np.exp(-2*(np.pi**2)*((x-x0)**2)*sigma**2)
|
|
||||||
return y
|
|
||||||
File diff suppressed because one or more lines are too long
@ -1,248 +0,0 @@
|
|||||||
import math
|
|
||||||
import numpy as np
|
|
||||||
from matplotlib import pyplot as plt
|
|
||||||
from scipy.signal import find_peaks
|
|
||||||
from scipy.optimize import curve_fit
|
|
||||||
|
|
||||||
|
|
||||||
def show_acquisition(recon_sum_norm, recon_sum, recon_sum_normalized, x0, y0):
|
|
||||||
|
|
||||||
plt.figure(figsize=(16,16), tight_layout=True)
|
|
||||||
|
|
||||||
plt.subplot(131)
|
|
||||||
plt.imshow(recon_sum_norm, cmap='gray')
|
|
||||||
plt.title('Normalization image')
|
|
||||||
|
|
||||||
plt.subplot(132)
|
|
||||||
plt.imshow(recon_sum, cmap='gray')
|
|
||||||
plt.scatter(x=y0,y=x0,s=8,color='r')
|
|
||||||
plt.title('Acquired image')
|
|
||||||
|
|
||||||
plt.subplot(133)
|
|
||||||
plt.imshow(recon_sum_normalized, cmap='gray')
|
|
||||||
plt.title('Normalized image')
|
|
||||||
|
|
||||||
|
|
||||||
def show_radii(img, x0, y0, R_MAX, R_MIN, title=None):
|
|
||||||
|
|
||||||
plt.figure(figsize=(8,8))
|
|
||||||
plt.imshow(img, cmap='gray')
|
|
||||||
plt.scatter(x=y0, y=x0, s=4, color='r')
|
|
||||||
|
|
||||||
if title is not None:
|
|
||||||
plt.title(title)
|
|
||||||
|
|
||||||
theta = np.arange(0,360,0.1)
|
|
||||||
|
|
||||||
x = np.zeros(len(theta))
|
|
||||||
y = np.zeros(len(theta))
|
|
||||||
|
|
||||||
for index,angle in enumerate(theta*np.pi/180):
|
|
||||||
x[index] = x0 + R_MAX*np.sin(angle)
|
|
||||||
y[index] = y0 + R_MAX*np.cos(angle)
|
|
||||||
|
|
||||||
plt.scatter(y,x,s=4)
|
|
||||||
|
|
||||||
for index,angle in enumerate(theta*np.pi/180):
|
|
||||||
x[index] = x0 + R_MIN*np.sin(angle)
|
|
||||||
y[index] = y0 + R_MIN*np.cos(angle)
|
|
||||||
|
|
||||||
plt.scatter(y,x,s=4)
|
|
||||||
|
|
||||||
|
|
||||||
def get_freq(radius, Np):
|
|
||||||
'''
|
|
||||||
Returns spatial frequency in cycles/pixel
|
|
||||||
'''
|
|
||||||
freq = Np/(2*np.pi*radius)
|
|
||||||
return freq
|
|
||||||
|
|
||||||
|
|
||||||
def pix_to_mm(res_radius: int, siemens_radius: int, phys_mag: float, ext_r: int,
|
|
||||||
zoom:int=1) -> float:
|
|
||||||
|
|
||||||
return res_radius * siemens_radius * phys_mag / (ext_r * zoom)
|
|
||||||
|
|
||||||
|
|
||||||
def calculate_lpmm(radius_pix: int, siemens_freq: int, siemens_radius: int,
|
|
||||||
phys_mag: float, ext_r: int, zoom:int=1) -> float:
|
|
||||||
"""Calculates resolution in linepairs per millimter (lp/mm).
|
|
||||||
|
|
||||||
Args:
|
|
||||||
radius_pix (int):
|
|
||||||
Radius in pixels at which resolution is determined.
|
|
||||||
siemens_freq (int):
|
|
||||||
Amount of black black bars in the Siemens Star.
|
|
||||||
siemens_radius (int):
|
|
||||||
Siemens Star radius in mm.
|
|
||||||
phys_mag (float):
|
|
||||||
Physical magnification due to the system's optics.
|
|
||||||
ext_r (int):
|
|
||||||
External radius in pixels.
|
|
||||||
zoom (int, optional):
|
|
||||||
Zoom applied to the acquisition. If present, the external radius
|
|
||||||
(ext_r) must be the same as in the image without zoom and this
|
|
||||||
function will calculate the correct external radius after zooming.
|
|
||||||
Defaults to 1.
|
|
||||||
|
|
||||||
Returns:
|
|
||||||
float: resolution in lp/mm.
|
|
||||||
"""
|
|
||||||
|
|
||||||
|
|
||||||
radius_mm = pix_to_mm(radius_pix, siemens_radius, phys_mag, ext_r, zoom)
|
|
||||||
theta = 2 * math.pi / siemens_freq
|
|
||||||
c = 2 * radius_mm * math.sin(theta/2)
|
|
||||||
|
|
||||||
return 1/c
|
|
||||||
|
|
||||||
|
|
||||||
def calculate_contrast(maxima, minima):
|
|
||||||
|
|
||||||
Imax = np.median(maxima)
|
|
||||||
Imax_mean = np.mean(maxima)
|
|
||||||
Imax_std = np.std(maxima)
|
|
||||||
Imax_unc = Imax_std/np.sqrt(len(maxima))
|
|
||||||
|
|
||||||
Imin = np.median(minima)
|
|
||||||
Imin_mean = np.mean(minima)
|
|
||||||
Imin_std = np.std(minima)
|
|
||||||
Imin_unc = Imin_std/np.sqrt(len(minima))
|
|
||||||
|
|
||||||
contrast = (Imax-Imin)/(Imax+Imin)
|
|
||||||
|
|
||||||
dImax2 = (2*Imax/(Imax+Imin)**2)**2
|
|
||||||
dImin2 = (2*Imin/(Imax+Imin)**2)**2
|
|
||||||
|
|
||||||
contrast_unc = np.sqrt(dImax2 * (Imax_unc**2) + dImin2 * (Imin_unc**2))
|
|
||||||
|
|
||||||
return contrast, contrast_unc, Imax, Imin
|
|
||||||
|
|
||||||
|
|
||||||
def object_resolution(res_radius: int, siemens_radius: int, siemens_freq: int,
|
|
||||||
phys_mag: float, ext_r: int, zoom:int=1) -> float:
|
|
||||||
"""Calculates resolution in mm.
|
|
||||||
|
|
||||||
Args:
|
|
||||||
res_radius (int):
|
|
||||||
Radius in pixels at which resolution is determined.
|
|
||||||
siemens_radius (int):
|
|
||||||
Siemens Star radius in mm.
|
|
||||||
siemens_freq (int):
|
|
||||||
Amount of black black bars in the Siemens Star.
|
|
||||||
phys_mag (float):
|
|
||||||
Physical magnification due to the system's optics.
|
|
||||||
ext_r (int):
|
|
||||||
External radius in pixels.
|
|
||||||
zoom (int, optional):
|
|
||||||
Zoom applied to the acquisition. If present, the external radius
|
|
||||||
(ext_r) must be the same as in the image without zoom and this
|
|
||||||
function will calculate the correct external radius after zooming.
|
|
||||||
Defaults to 1.
|
|
||||||
|
|
||||||
Returns:
|
|
||||||
float:
|
|
||||||
Real resolution at the object plane.
|
|
||||||
"""
|
|
||||||
|
|
||||||
res_r = pix_to_mm(res_radius, siemens_radius, phys_mag, ext_r, zoom)
|
|
||||||
res = 2 * np.pi * res_r / siemens_freq
|
|
||||||
|
|
||||||
return res
|
|
||||||
|
|
||||||
|
|
||||||
def find_resolution(img, x0, y0, radii, interactive=False):
|
|
||||||
|
|
||||||
d_theta = 0.0001
|
|
||||||
theta = np.arange(0,2*np.pi, d_theta)
|
|
||||||
|
|
||||||
d = int(10*np.pi/180/d_theta * 2/3)
|
|
||||||
|
|
||||||
contrast = np.zeros(len(radii))
|
|
||||||
contrast_unc = np.zeros(len(radii))
|
|
||||||
|
|
||||||
for index, R in enumerate(radii):
|
|
||||||
|
|
||||||
values = np.zeros(len(theta))
|
|
||||||
x = np.around(x0 + R*np.cos(theta)).astype('int')
|
|
||||||
y = np.around(y0 + R*np.sin(theta)).astype('int')
|
|
||||||
|
|
||||||
for i in range(len(theta)):
|
|
||||||
values[i] = img[x[i],y[i]]
|
|
||||||
|
|
||||||
# Finding maxima and minima
|
|
||||||
maxima,_ = find_peaks(values,distance=d)
|
|
||||||
minima,_ = find_peaks(-values,distance=d)
|
|
||||||
|
|
||||||
contrast[index],contrast_unc[index],Imax,Imin = calculate_contrast(
|
|
||||||
values[maxima],
|
|
||||||
values[minima])
|
|
||||||
|
|
||||||
if interactive:
|
|
||||||
plt.figure()
|
|
||||||
plt.plot(theta, values, label='profile')
|
|
||||||
plt.scatter(theta[maxima], values[maxima], label='maxima')
|
|
||||||
plt.scatter(theta[minima], values[minima], label='minima')
|
|
||||||
plt.axhline(Imax, label=f'median maximum = {Imax:.2f}')
|
|
||||||
plt.axhline(Imin, label=f'median minimum = {Imin:.2f}')
|
|
||||||
plt.xlabel('theta (rad)')
|
|
||||||
plt.ylabel('Normalized intensity')
|
|
||||||
plt.title(f'R={R} pix, contrast={contrast[index]:.3f}')
|
|
||||||
plt.legend()
|
|
||||||
plt.waitforbuttonpress()
|
|
||||||
plt.close()
|
|
||||||
|
|
||||||
ind = np.abs(contrast - 0.1).argmin()
|
|
||||||
|
|
||||||
# Forcing the contrast to be at least 0.1
|
|
||||||
if contrast[ind] < 0.1:
|
|
||||||
ind += 1
|
|
||||||
|
|
||||||
res_radius = radii[ind]
|
|
||||||
res_MTF = contrast[ind]
|
|
||||||
|
|
||||||
print(f'Found resolution at R={res_radius} pix, MTF={res_MTF}')
|
|
||||||
|
|
||||||
return res_radius, res_MTF, contrast, contrast_unc
|
|
||||||
|
|
||||||
|
|
||||||
def plot_MTF_radius(radii, contrast, contrast_unc=None):
|
|
||||||
|
|
||||||
plt.figure()
|
|
||||||
plt.plot(radii,contrast, label='MTF')
|
|
||||||
plt.axhline(0.1, label='Resolution limit') # Resolution limit at 10% of the MTF
|
|
||||||
|
|
||||||
if contrast_unc is not None:
|
|
||||||
plt.errorbar(radii,contrast,yerr=contrast_unc, label='MTF error')
|
|
||||||
|
|
||||||
plt.xlabel('Radius (pix)')
|
|
||||||
plt.ylabel('Contrast')
|
|
||||||
plt.legend()
|
|
||||||
|
|
||||||
|
|
||||||
def plot_MTF_freq(radii, contrast, contrast_unc=None):
|
|
||||||
|
|
||||||
freqs = [get_freq(R,36) for R in radii]
|
|
||||||
|
|
||||||
plt.figure()
|
|
||||||
plt.plot(freqs,contrast, label='MTF')
|
|
||||||
plt.axhline(0.1, label='Resolution limit') # Resolution limit at 10% of the MTF
|
|
||||||
|
|
||||||
if contrast_unc is not None:
|
|
||||||
plt.errorbar(freqs,contrast,yerr=contrast_unc, label='MTF error')
|
|
||||||
|
|
||||||
plt.xlabel('f (cycles/pixel)')
|
|
||||||
plt.ylabel('Contrast')
|
|
||||||
plt.title('MTF')
|
|
||||||
plt.legend()
|
|
||||||
|
|
||||||
|
|
||||||
def reciprocal_func(x, A):
|
|
||||||
return A/x
|
|
||||||
|
|
||||||
|
|
||||||
def resolution_curve_coeffs(zooms, resolutions):
|
|
||||||
|
|
||||||
popt, pcov = curve_fit(reciprocal_func, zooms, resolutions)
|
|
||||||
|
|
||||||
return popt[0]
|
|
||||||
Loading…
Reference in New Issue
Block a user