248 lines
7.2 KiB
Python
248 lines
7.2 KiB
Python
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] |