5050L-MTF/show.py

import numpy as np
import matplotlib.pyplot as plt
from scipy.ndimage import gaussian_filter
from scipy.fft import fft, fftshift, fftfreq
from ROI import interface_roi
from line import interface_corr


def generate_test_edge(width=512, height=512, edge_blur=2.0):
    """生成测试用的边缘图像，用于演示MTF计算"""
    # 创建理想的锐利边缘
    edge_image = np.ones((height, width))
    edge_position = width // 2
    edge_image[:, :edge_position] = 0  # 左侧为暗区，右侧为亮区

    # 添加高斯模糊模拟实际成像系统的模糊
    blurred_edge = gaussian_filter(edge_image, sigma=edge_blur)

    return blurred_edge


def extract_esf(edge_image, roi_size=100):
    """从边缘图像中提取边缘扩散函数(ESF)"""
    height, width = edge_image.shape

    # 找到边缘位置（亮度变化最大的列）
    edge_profile = np.mean(edge_image, axis=0)  # 沿垂直方向平均得到水平方向的亮度分布
    edge_gradient = np.abs(np.gradient(edge_profile))
    edge_position = np.argmax(edge_gradient)

    # 提取感兴趣区域(ROI)，围绕边缘位置
    start = max(0, edge_position - roi_size // 2)
    end = min(width, edge_position + roi_size // 2)

    # 提取边缘区域并沿垂直方向平均，得到ESF
    esf_region = edge_image[:, start:end]
    esf = np.mean(esf_region, axis=0)  # 沿垂直方向平均

    # 归一化ESF到[0, 1]范围
    esf = (esf - np.min(esf)) / (np.max(esf) - np.min(esf))

    return esf, start, end


def calculate_mtf(esf, pixel_size=1.0):
    """
    从边缘扩散函数(ESF)计算调制传递函数(MTF)

    参数:
        esf: 边缘扩散函数
        pixel_size: 像素尺寸(mm)，用于计算实际空间频率

    返回:
        frequencies: 空间频率 (线对/毫米)
        mtf: 调制传递函数值
    """
    # 对ESF求导得到线扩散函数(LSF)
    lsf = np.gradient(esf)

    # 对LSF进行傅里叶变换得到SFR
    n = len(lsf)
    lsf_fft = fft(lsf)
    sfr = fftshift(lsf_fft)

    # 计算幅度谱并归一化得到MTF
    mtf = np.abs(sfr)
    mtf = mtf / mtf[int(n / 2)]  # 直流分量归一化到1

    # 计算空间频率 (线对/毫米)
    frequencies = fftfreq(n, d=pixel_size)
    frequencies = fftshift(frequencies)
    frequencies = frequencies[frequencies >= 0]  # 只保留正频率

    # 只保留正频率部分的MTF
    mtf = mtf[int(n / 2):]

    return frequencies, mtf, lsf


def trans(frequencies,pixel_size,x_unit='LP/mm',H_mm=1.0):
    unit = x_unit
    F = 50.0  # 焦距，单位毫米
    if unit == 'LP/mm':
        x = frequencies
        xlabel = 'Spatial Frequency (LP/mm)'
    elif unit == 'L/mm':
        x = frequencies * 2.0  # 1 LP = 2 L
        xlabel = 'Spatial Frequency (L/mm)'
    elif unit == 'Cycles/mm':
        x = frequencies  # Cycles/mm 等同于 LP/mm
        xlabel = 'Spatial Frequency (Cycles/mm)'
    elif unit == 'Cycles/pixel':
        x = frequencies * pixel_size  # = (LP/mm) * (mm/pixel)
        xlabel = 'Spatial Frequency (Cycles/pixel)'
    elif unit == 'LP/PH':
        x = frequencies * H_mm  # = (LP/mm) * H(mm)
        xlabel = 'Spatial Frequency (LP/PH)'
    elif unit == 'LW/PH':
        x = frequencies * (2.0 * H_mm)  # = (LP/mm) * 2 * H(mm)
        xlabel = 'Spatial Frequency (LW/PH)'
    elif unit == 'C/Mrad':
        x = (F *frequencies)/1000  # = (LP/mm) * H(mm) / 1000
        xlabel = 'Spatial Frequency (Cycles/Mrad)'

    else:
        x = frequencies
        xlabel = f'Spatial Frequency ({unit})'

    return x,xlabel


def plot_results(edge_image, esf, frequencies, mtf, lsf,image_high=480,x_range=None, x_unit='LP/mm',
                 pixel_size=1.0):
    """可视化结果：边缘图像、ESF和MTF曲线"""
    fig, axes = plt.subplots(2, 2, figsize=(18, 5))
    axes = axes.flatten()  # 展平成一维数组，方便索引

    # 显示边缘图像
    axes[0].imshow(edge_image, cmap='gray')
    axes[0].set_title('Edge Image')
    axes[0].axis('off')

    # 显示边缘扩散函数(ESF)
    axes[1].plot(esf)
    axes[1].set_title('Edge Spread Function (ESF)')
    axes[1].set_xlabel('Position (pixels)')
    axes[1].set_ylabel('Normalized Intensity')
    axes[1].grid(True)

    # 显示MTF曲线（支持多种空间频率单位）
    # ph_pixels = edge_image.shape[0]
    h_mm = image_high * pixel_size  # 画幅高度（mm）
    unit = x_unit

    x,xlabel = trans(frequencies,pixel_size,unit,h_mm)

    axes[2].plot(x, mtf)
    axes[2].set_title('Modulation Transfer Function (MTF)')
    axes[2].set_xlabel(xlabel)
    axes[2].set_ylabel('MTF Value')
    axes[2].set_ylim(0, 1.1)
    # 根据单位设置横坐标显示范围
    if x_range is not None:
        axes[2].set_xlim(0,x_range)
    else:
        axes[2].set_xlim(x.min(), x.max())

    axes[2].grid(True)


    axes[2].grid(True)

    # 打印一些关键频率的MTF值
    key_mtf_values = [0.5, 0.2, 0.1]
    print("关键MTF对应的空间频率:")

    for t in key_mtf_values:
        f = find_freq_for_mtf(frequencies, mtf, t)
        if f is not None:
            f,xlabel = trans(f,pixel_size,unit,h_mm)
            print(f"MTF={t:.3f}: {f:.3f} {xlabel}")
        else:
            print(f"MTF={t:.3f}: 未找到对应频率")



    # 显示线扩散函数(LSF)曲线
    axes[3].plot(lsf)
    axes[3].set_title('Line Spread Function (LSF)')
    axes[3].set_xlabel('Position (pixels)')
    axes[3].set_ylabel('LSF Value')
    axes[3].grid(True)

    plt.tight_layout()
    plt.tight_layout(pad=2.0)
    plt.subplots_adjust(top=0.9, bottom=0.1, left=0.08, right=0.95, wspace=0.3, hspace=0.5)
    plt.show()

def find_freq_for_mtf(freqs, mtf_vals, target):
    freqs = np.asarray(freqs)
    mtf_vals = np.asarray(mtf_vals)
    idxs = np.where((mtf_vals[:-1] >= target) & (mtf_vals[1:] <= target))[0]
    if idxs.size == 0:
        return None
    i = idxs[0]
    m1, m2 = mtf_vals[i], mtf_vals[i + 1]
    f1, f2 = freqs[i], freqs[i + 1]
    if m2 == m1:
        return f1
    w = (target - m1) / (m2 - m1)
    return f1 + w * (f2 - f1)

def main(image_path=None):
    # 生成测试边缘图像
    edge_image = None
    # image_path = r'output/result_hough.bmp'  # 将此路径改为你的图像文件路径（Windows 风格路径）

    # 读取图像（支持灰度或 RGB，使用 matplotlib.imread）
    img = plt.imread(image_path)

    # 如果是带 alpha 的 RGBA，去掉 alpha
    if img.ndim == 3 and img.shape[2] == 4:
       img = img[..., :3]

    # 如果是彩色图像，转换为灰度
    if img.ndim == 3:
       edge_image = np.mean(img, axis=2)
    else:
       edge_image = img.copy()

    # 若像素值为 0-255，归一化到 0-1
    if edge_image.max() > 1.0:
       edge_image = edge_image.astype(np.float32) / 255.0

    # 可选：确保图像足够大或裁剪/缩放以适应 ESF 提取（根据需要调整）
    # 例如裁剪到 512x512： edge_image = edge_image[:512, :512]

    # 提取边缘扩散函数
    # esf, _, _ = extract_esf(edge_image, roi_size=200)
    x,y = edge_image.shape
    if x > y:
         min_len = y
    else:
         min_len = x

    esf, _, _ = extract_esf(edge_image, roi_size=min_len)

    # # 计算MTF，假设像素尺寸为0.012毫米
    # frequencies, mtf, lsf = calculate_mtf(esf, pixel_size=0.012)
    #
    # # 显示结果
    # plot_results(edge_image, esf, frequencies, mtf, lsf)
    # 计算MTF，假设像素尺寸为0.012毫米
    pixel_size = 0.012
    frequencies, mtf, lsf = calculate_mtf(esf, pixel_size=pixel_size)

    # 显示结果
    x_unit = 'C/Mrad'  # 可选: 'C/Mrad','LP/mm','L/mm','Cycles/mm','Cycles/pixel','LP/PH','LW/PH'
    plot_results(edge_image, esf, frequencies, mtf, lsf, x_unit=x_unit, pixel_size=pixel_size)




if __name__ == "__main__":
    roi_img = interface_roi()
    input_path = interface_corr(roi_img)
    main(input_path)
    main(r"output/result_hough.bmp")