pythonmath

Fourier Series Implementation cannot approximate batman shape


I tried to implement a formula, from which a coefficients of Fourier Series could be calculated. (I used 3B1B's video about it: Video) and writing code for that, my first test subject was singular contour of batman logo, I first take a binary picture of batman logo and use marching squares algorithm to find contour of it. after that i rescale values and get this results: POINTS FROM THE MARCHING SQUARES ALGORITHM

And Here is Code for creating this points: (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)

Now we go to the main part, where we implement the fourier series algorithm. I divide the time interal (t) into the amount of points provided and i make assumtion that all of that points relative to t have same distances between eachother. I use approximation of integral as the sum of the points as provided into the formula.FORMULA APPROXIMATION

And Here is code implementing it (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

(NOTE: For this code t1 is basically deltaT, because it equals to 1/len(complex_points) And Now, in the next slide i animate whole process, where i also wrote additional code snippet for creating a gif. If my implementation were correct it shouldn't have anu difficulty creating a batman shape, but we can observe really weird phenomenons throught the gif.

this is code snippet for this part

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'.")

Now this is the result GIF

Observation #1 when coefficients number is 0, 1, 2 or 3 it doesnt draw anything.

Observation #2

As coefficients number raises, we get the wobbly circular shape, where the lower part of the image is slightly more identical tot he original imagine, but messes up on its wings

Observation #3

As we get closer to the len(complex_numbers), the situacion changes and we get this weird shapes, different from circular

Observation #4

When we surpass the len(complex_number), it draws a random gibberish

Observation #5

When the number of the divisions inside the t value in animation.py code is altered we get completely different images.

EDIT 1

here is actual .txt data provided for further testing.

https://pastebin.com/Q51pT09E

After all of this information given, can you guys help me out whats wrong with my code


Solution

  • In the definition of the Fourier series, you can see that n goes from negative infinity to positive infinity. The issue in your code is that you forgot to compute the coefficients associated with negative values of n.

    Here is a simple example that shows how to compute the coefficients (from -50 to 50) associated with an ellipse, and build a curve from them:

    import numpy as np
    import matplotlib.pyplot as plt
    
    def get_ellipse():
        t = np.linspace(0, 1, 100)
        X = 2 * np.cos(2 * np.pi * t)
        Y = np.sin(2 * np.pi * t)
        return (X, Y)
    
    def calculate_Fourier(X, Y, N):
        complex_points = [complex(x, y) for x, y in zip(X, Y)]
        t = np.linspace(0, 1, len(complex_points), endpoint=False)
        coefficients = np.zeros(2 * N + 1, dtype=complex)
        for i in range(len(coefficients)):
            n = i - N
            coefficients[i] = np.sum(complex_points * np.exp(-2j * np.pi * n * t) * t[1])
        return coefficients
    
    def build_curve(coefficients, num_points):
        N = (len(coefficients) - 1) / 2
        t = np.linspace(0, 1, num_points)
        curve = np.zeros(num_points, dtype=complex)
        for i in range(len(coefficients)):
            c = coefficients[i]
            n = i - N
            curve += c * np.exp(2j * np.pi * n * t)
        return curve
    
    X, Y = get_ellipse()
    coefficients = calculate_Fourier(X, Y, 50)
    curve = build_curve(coefficients, 50)
    
    plt.plot(curve.real, curve.imag, color="blue")
    plt.show()
    

    Result:

    enter image description here

    Remark: if the number of coefficients is too high, you will get numerical instability.