siemens_star_analysis.py

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]