我尝试实现一个计算傅里叶级数系数的公式（参考了3Blue1Brown的相关视频），并编写代码测试。我的第一个测试对象是蝙蝠侠标志的单轮廓——先提取蝙蝠侠标志的二值图像，用Marching Squares算法提取轮廓，缩放后得到了如下点集：

提取轮廓点的代码（Contour_Classifier.py） #

import numpy as np
import matplotlib.pyplot as plt
from skimage import measure, draw

def read_binary_image(file_path):
    # Open the file and read line by line
    with open(file_path, 'r') as file:
        lines = file.readlines()

    height, width = len(lines), len(lines[0])
    print(height, width)
    # Process lines into a 2D numpy array
    image_data = []

    for i in range(height + 2):
        arr = []
        for j in range(width + 2):
            arr.append(0)
        image_data.append(arr)

    for i in range(2, height + 1):
        for j in range(2, width + 1):
            if(lines[i - 2][j - 2] != '1'):
                image_data[i][j] = 0
            else:
                image_data[i][j] = 1

    # Convert list to numpy array for easier manipulation
    image_array = np.array(image_data)

    return image_array

def display_image(image_array):
    # Display the binary image using matplotlib
    plt.imshow(image_array, cmap="gray")
    plt.axis('off')  # Hide axes
    plt.show()

# Example usage
file_path = 'KOREKT\images\sbetmeni.txt'  # Replace with the path to your file
image_array = read_binary_image(file_path)
#display_image(image_array)

#----------------------------------------------------------------------------------------------------------
#-------------------------------------------Finding Contours-----------------------------------------------
#----------------------------------------------------------------------------------------------------------

contours = measure.find_contours(image_array, level=0.5, positive_orientation='high')

fixed_contours = []
for contour in contours:
    fixed_contour = np.column_stack((contour[:, 1], contour[:, 0]))  # Swap (row, column) to (column, row)
    fixed_contour[:, 1] = image_array.shape[0] - fixed_contour[:, 1]  # Invert the y-axis
    # Normalize coordinates between [0, 1]
    fixed_contour[:, 0] /= image_array.shape[1]  # Normalize x (width)
    fixed_contour[:, 1] /= image_array.shape[0]  # Normalize y (height)

    fixed_contour[:, 0] *= 250  # Normalize x (width)
    fixed_contour[:, 1] *= 250  # Normalize y (height)

    fixed_contours.append(fixed_contour)
contours = fixed_contours

print(fixed_contours[0])

def visualize_colored_contours(contours, title="Colored Contours"):
    # Create a plot
    plt.figure(figsize=(8, 8))

    for i, contour in enumerate(contours):
        # Extract X and Y coordinates
        x, y = zip(*contour)
        # Plot the points with a unique color
        plt.plot(x, y, marker='o', label=f'Contour {i+1}')

    plt.title(title)
    plt.xlabel("X")
    plt.ylabel("Y")
    plt.legend()
    plt.grid(True)
    plt.axis("equal")
    plt.show()

# Visualize the normalized contours
visualize_colored_contours(contours)

傅里叶系数计算实现 #

接下来是核心的傅里叶级数算法实现。我将时间区间t按点的数量划分，假设每个点对应的t间隔相等，用点的求和来近似积分，公式近似如下：

对应的实现代码（Fourier_Coefficients.py）：

import numpy as np

def calculate_Fourier(points, num_coefficients):
    complex_points = []
    for point in points:
        complex_points.append(point[0] + 1j * point[1])


    t = np.linspace(0, 1, len(complex_points), endpoint=False)

    c_k = np.zeros(num_coefficients, dtype=np.complex128)

    for i in range(num_coefficients):
        c_k[i] = np.sum(complex_points * np.exp(-2j * np.pi * i * t) * t[1])

    return c_k

注：代码中的t[1]其实就是Δt，等于1/len(complex_points)

动画生成代码 #

我还编写了代码将傅里叶逼近过程生成为GIF，理论上如果实现正确，应该能还原蝙蝠侠形状，但实际生成的GIF出现了很多奇怪的现象。动画代码如下：

import numpy as np
import matplotlib.pyplot as plt
import imageio
from Fourier_Coefficients import calculate_Fourier
from Countour_Classifier import contours



# List to store file names for GIF creation
png_files = []

# Generate plots iteratively
for i in range(len(contours[0])):


    contour_coefficients = []

    for contour in contours:
        contour_coefficients.append(calculate_Fourier(contour, i))

    # Fourier coefficients (complex numbers) and frequencies
    coefficients = contour_coefficients[0]  # First contour
    frequencies = np.arange(len(coefficients))

    # Time parameters
    t = np.linspace(0, 1, len(coefficients))  # One period
    curve = np.zeros(len(t), dtype=complex)

    # Use the first (i + 1) coefficients
    for j in range(len(coefficients)):
        c, f = coefficients[j], frequencies[j]
        curve += c * np.exp(1j * 2 * np.pi * f * t)

    # Plotting
    plt.figure(figsize=(8, 8))
    plt.plot(curve.real, curve.imag, label="Trajectory", color="blue")
    plt.scatter(0, 0, color="black", label="Origin")
    plt.axis("equal")
    plt.title(f"Fourier Series with {i + 1} Coefficients")
    plt.xlabel("Real Part (X)")
    plt.ylabel("Imaginary Part (Y)")
    plt.legend()
    plt.text(-0.5, -0.5, f"Using {i + 1} coefficients", fontsize=12, color="red")

    # Save the figure as a PNG file
    filename = f"fourier_{i + 1}_coefficients.png"
    plt.savefig(filename)
    plt.close()

    # Append the file name to the list
    png_files.append(filename)

# Create a GIF from the PNG files
gif_filename = "fourier_series.gif"
with imageio.get_writer(gif_filename, mode='I', duration=0.5) as writer:
    for filename in png_files:
        image = imageio.imread(filename)
        writer.append_data(image)

print("Plots saved as PNG files and GIF created as 'fourier_series.gif'.")

生成结果与问题观察 #

实际生成的GIF出现了诸多异常，我总结了以下观察：

• 观察1：当系数数量为0、1、2或3时，没有任何绘制内容。
• 观察2：随着系数数量增加，得到一个晃动的圆形，图像下半部分和原图有一点相似，但翅膀部分完全不对。
• 观察3：当系数数量接近len(complex_numbers)时，形状发生变化，不再是圆形，变成奇怪的样子。
• 观察4：当系数数量超过len(complex_number)时，绘制出完全随机的乱码形状。
• 观察5：修改动画代码中t的划分数量时，得到的图像完全不同。

补充：我已提供测试用的txt数据供进一步排查。

以上就是我的全部信息，有没有大佬能帮我找出代码里的问题？

备注：内容来源于stack exchange，提问作者NikoMolecule