图1(a)和(b)咨询

[xirihao]

你好，我认为您的文章非常有意义，而且希望您可以做出更好的工作。
我有个问题，请问图1的(a)和(b)怎么做的呀？希望可以分享以下代码

[shangshang0912]

您的邮件已经收到。我会及时做出回复。

[xirihao]

非常希望您的分享，万分感谢。

[xwmaxwma]

我也想咨询图1的(a)和(b)如何画的。希望收到您的回复

[myu1999]

非常希望您的分享，万分感谢。

您的邮件已经收到。我会及时做出回复。

非常希望您的分享，万分感谢。

请问你收到作者的邮件回复了吗，方便分享一下怎么画的吗？

[shangshang0912]

您的邮件已经收到。我会及时做出回复。

[dalision]

I have tried some, but the curve is much smoother than the one in the paper. Any further discussions?

import torch
import numpy as np
import matplotlib.pyplot as plt
from PIL import Image
from matplotlib.cm import get_cmap
def load_image_grayscale(image_path, size=None):
    """读取图像并转为灰度Tensor."""
    img = Image.open(image_path).convert('L')  # 转换为灰度图
    if size:
        img = img.resize(size)
    img_tensor = torch.from_numpy(np.array(img)).float()
    return img_tensor

def compute_fft(img_tensor):
    """对输入灰度图做二维FFT并中心化."""
    # 扩展维度: (H,W) -> (1,1,H,W)
    img_tensor = img_tensor.unsqueeze(0).unsqueeze(0)
    fft_result = torch.fft.fft2(img_tensor)
    fft_shifted = torch.fft.fftshift(fft_result)
    # squeeze 回到 (H, W)
    return fft_shifted.squeeze()

def compute_frequency_grid(width, height):
    """
    生成 2D 的角频率网格 (单位: 弧度/采样).
    np.fft.fftfreq(N) 返回范围大概 [-0.5, 0.5), 我们乘以 2*pi 得到 [-π, +π).
    """
    freq_x = np.fft.fftfreq(width) * 2 * np.pi   # [-pi, pi)
    freq_y = np.fft.fftfreq(height) * 2 * np.pi  # [-pi, pi)
    freq_x_shifted = np.fft.fftshift(freq_x)     # 让数组中心为 0 频率
    freq_y_shifted = np.fft.fftshift(freq_y)
    X, Y = np.meshgrid(freq_x_shifted, freq_y_shifted)
    R = np.sqrt(X**2 + Y**2)  # 径向坐标
    return R

def radial_average(log_amp_spectrum, R, bin_size=0.01):
    """
    对 2D 对数振幅谱做径向平均，返回频率刻度和平均值。
    bin_size = 0.01 -> 在 [0, pi] 范围内每隔0.01做一个分箱.
    """
    data_2d = log_amp_spectrum.numpy()
    
    # 显示 2D 对数振幅谱
    # plt.figure(figsize=(8, 6))
    # plt.imshow(data_2d,
    #            cmap='viridis',
    #            extent=[-np.pi, np.pi, -np.pi, np.pi])  # 这里extent配合了2*pi*fftfreq的做法
    # plt.colorbar(label="Log Amplitude")
    # plt.title("2D FFT Log Amplitude Spectrum")
    # plt.xlabel("Frequency (radians/sample)")
    # plt.ylabel("Frequency (radians/sample)")
    # plt.show()
    
    # 设置频率分桶 (0 到 π)
    freq_bins = np.arange(0, np.pi + bin_size, bin_size)
    radial_means = np.zeros(len(freq_bins) - 1)

    # 逐个桶计算平均
    for i in range(len(freq_bins) - 1):
        low, high = freq_bins[i], freq_bins[i+1]
        mask = (R >= low) & (R < high)
        if np.any(mask):
            radial_means[i] = data_2d[mask].mean()
        else:
            radial_means[i] = 0

    # （可选）减去第一个桶的平均值，相当于去掉 DC 分量的影响
    radial_means = radial_means - radial_means[0]
    
    # 取分桶中心点作为最终频率刻度
    freq_bins_center = 0.5 * (freq_bins[:-1] + freq_bins[1:])
    return freq_bins_center, radial_means

def plot_log_amplitude_vs_frequency(freq, radial_mean):
    """绘制“频率 vs. 径向平均对数振幅”曲线."""
    plt.figure(figsize=(10, 6))
    plt.plot(freq/np.pi, radial_mean, linestyle='-')
    plt.title('Log Amplitude Spectrum vs Frequency')
    plt.xlabel('Frequency (radians/sample)')  # 这里与上面保持一致
    plt.ylabel('Log Amplitude')
    plt.grid(True)
    plt.show()

def calculate_relative_log_frequency(image_path):
    img_tensor = load_image_grayscale(image_path, size=(512, 512))
    
    # 2. 执行二维 FFT
    fft_shift = compute_fft(img_tensor)
    
    # 3. 计算对数振幅谱
    log_amplitude_spectrum = torch.log(torch.abs(fft_shift) + 1)
    
    # 4. 计算频率网格 R
    h, w = log_amplitude_spectrum.shape
    R = compute_frequency_grid(w, h)
    
    # 5. 计算径向平均
    freq_bins_center, radial_mean_val = radial_average(log_amplitude_spectrum, R, bin_size=0.025*np.pi)
    return freq_bins_center/np.pi, radial_mean_val


if __name__ == "__main__":
    cmap = get_cmap('Blues') 
    image_filenames = ["1.PNG", "2.PNG", "3.PNG"]
    new_labels = ["1", "2", "3"]
    for i, filename in enumerate(image_filenames):
        relative_log_frequency, frequency_values = calculate_relative_log_frequency(filename)
        k = i / len(image_filenames)
        color = cmap(k+0.2)
        label = new_labels[i]
        plt.plot(relative_log_frequency,frequency_values, color=color, label=label, linewidth=2.5, alpha=0.8)


    plt.xlabel('Frequency',fontsize=16)
    plt.ylabel('Relative Log amplitude',fontsize=16)
    # plt.xlim(0.0, 1.0)
    plt.grid(False)
    plt.legend()

    x_ticks = np.linspace(0, 1, 6)
    x_tick_labels = [f'{x:.1f}π' for x in x_ticks]
    plt.xticks(x_ticks, x_tick_labels, fontsize=12)
    
    plt.minorticks_off()
    plt.show()