[PYTHON] Drawing with Bezier curve and Fourier transform

This is the first post. I intend to write it so that even high school students can understand it. Please forgive that the implementation is difficult to read.

The characters in the picture are Suzu Kagura (green) and Yamatoiori (light blue) from .LIVE. The work on this page is a derivative work. Intellectual property rights related to characters belong to Upland Co., Ltd.

Work example

This is the one. test.gif

Rough flow

Make a collection of ordered dots that form a picture when tied ↓ Fast Fourier transform of x and y coordinates respectively ↓ Can be drawn by superposition of circular motion



If you already have a collection of points that will make a picture when connected, skip to the Discrete Fourier Transform.

Draw a picture in PowerPoint


In PowerPoint "Curve", write a picture __ one stroke __. Writing with one stroke is quite difficult.

Rewrite .pptx to .zip and unzip

If you rewrite .pptx to .zip and decompress it, the following xml will appear. I referred to I wish the normal distribution could be drawn in PowerPoint.

Unzipped file/ppt/slides/slide1.xml

    <a:path w="6982104" h="6267205">
            <a:pt x="2424619" y="3964459"/>
            <a:pt x="2459747" y="3967701"/>
            <a:pt x="2494875" y="3970944"/>
            <a:pt x="2534866" y="3945004"/>
            <a:pt x="2574858" y="3919064"/>
            <a:pt x="2615930" y="3853132"/>
            <a:pt x="2664568" y="3808817"/>
            <a:pt x="2713206" y="3764502"/>
            <a:pt x="2757521" y="3720187"/>
            <a:pt x="2826695" y="3679115"/>

Cubic Bézier curves are specified by four points. See "Cubic Bézier curves" in English wikipedia for more information. The meaning of the contents of xml


    <a:pt The first control point and the fourth control point on the previous curve/>
    <a:pt 2nd control point/>
    <a:pt 3rd control point/>
    <a:pt 4th control point and 1st control point on the curve after/>

It seems that it looks like this.

Bezier curve implementation

Control points as $ P_ {0}, P_ {1}, P_ {2}, P_ {3} $ Bezier curve B(t) = (1-t)^{3}P_{0}+3(1-t)^{2}tP_{1}+3(1-t)t^{2}P_{2}+t^{3}P_{3},\quad 0\leq t\leq 1 Is expressed as. I implemented it as it is.


import numpy as np

class BezFunc:
    #p is a list of 4 control points
    def __init__(self, p):
        self.p0 = np.array(p[0])
        self.p1 = np.array(p[1])
        self.p2 = np.array(p[2])
        self.p3 = np.array(p[3])
    # 0<=t<=1
    def __call__(self, t):
        return ((1-t)**3)*self.p0+3*((1-t)**2)*t*self.p1+3*(1-t)*(t**2)*self.p2+(t**3)*self.p3

Extract control points from xml and draw a curve

Use BeautifulSoup to retrieve control points from xml.


from bs4 import BeautifulSoup
import numpy as np

class Points:
    def __init__(self, xml_path):
        with open(xml_path, 'r') as f:
            soup = BeautifulSoup(f, 'xml')

        # [[Control point 1,Control point 2,Control point 3,Control point 4], [Control point 1,Control point 2, …
        bez_points_list = []

        temp_pt = soup.find('a:pathLst').find('a:moveTo').find('a:pt')
        temp_point = (int(temp_pt['x']), int(temp_pt['y']))

        for temp_cubicBezTo in soup.find('a:pathLst').find_all('a:cubicBezTo'):
            temp_list = [temp_point]
            for temp_pt in temp_cubicBezTo.find_all('a:pt'):
                temp_point = (int(temp_pt['x']), int(temp_pt['y']))
        self.bez_points_list = bez_points_list

    # n :Division number
    def calc_points(self, n):
        points = []

        for bez_points in self.bez_points_list:
            f = BezFunc(bez_points)
            for i in range(n):
                t = i/n

        self.points = np.array(points)
    def get_points(self):
        return self.points

How to use


xml_path = 'presentation/ppt/slides/slide1.xml'
p = Points(xml_path)
p.calc_points(n=5) #5 divisions(suitable)
points = p.get_points()
array([[2424619.   , 3964459.   ], #Coordinates of the first point
       [2445734.704, 3966170.848], #Coordinates of the second point
       [2467083.832, 3966482.104], #Coordinates of the third point
       [3557916.696, 5271117.176],
       [3503553.224, 5255947.584],
       [3548753.368, 5184643.528]])

points is an array ofshape = (number of points, 2)numpy. If you connect these in order from the top, you can draw a picture.


import matplotlib.pyplot as plt

plt.figure(figsize=(8, 5), dpi=200)
plt.plot(points[:,0], points[:,1])


Discrete Fourier transform

Discrete Fourier transform and inverse transform

F(t) = \sum_{x=0}^{N-1}f(x)e^{-i\frac{2\pi t x}{N}}\\

f(x) = \frac{1}{N}\sum_{t=0}^{N-1}F(t)e^{i\frac{2\pi t x}{N}}

Then, when the discrete Fourier transform is performed and the inverse transform is performed, the original function is restored. Think of $ F (t) and \ quad f (x) $ as an array of N elements, respectively. this is,

Discrete Fourier transform

F(t) = \frac{1}{N}\sum_{x=0}^{N-1}f(x)e^{-i\frac{2\pi t x}{N}}\\

f(x) = \sum_{t=0}^{N-1}F(t)e^{i\frac{2\pi t x}{N}}

But the same is true.


F(t) = A(t)( \cos \phi (t) + i \sin \phi (t))

Euler's formula

e^{i\theta}=\cos \theta + i \sin \theta

With, the inverse transformation

f(x)&= \sum_{t=0}^{N-1}F(t)e^{i\frac{2\pi t x}{N}}\\
&= \sum_{t=0}^{N-1}A(t)( \cos \phi (t) + i \sin \phi (t))(\cos \frac{2\pi t x}{N} + i \sin \frac{2\pi t x}{N})\\
&= \sum_{t=0}^{N-1}A(t)(\cos (\phi (t) +\frac{2\pi t x}{N}) + i \sin (\phi (t) + \frac{2\pi t x}{N}))\\
&= \sum_{t=0}^{N-1}A(t)\cos (\phi (t) +\frac{2\pi t x}{N}) + i \sum_{t=0}^{N-1}A(t) \sin (\phi (t) + \frac{2\pi t x}{N})\\


And if $ f (x) $ is a real number, the second term will be zero, so

f(x)= \sum_{t=0}^{N-1}A(t)\cos (\phi (t) +\frac{2\pi t x}{N})

You can see that the inverse transformation can be represented by the projection of the superposition of circular motions onto a straight line.

Fast Fourier transform

Discrete Fourier transform formula

F(t) = \frac{1}{N}\sum_{x=0}^{N-1}f(x)e^{-i\frac{2\pi t x}{N}}

You can implement it as is, but using the Fast Fourier Transform (FFT) will result in a very short execution time. (Experience story)

numpy fft is

F(t) = \sum_{x=0}^{N-1}f(x)e^{-i\frac{2\pi t x}{N}}

Is designed to be calculated, so divide by N.


import numpy as np

F = np.fft.fft(One-dimensional array) / len(One-dimensional array)

This will give you a one-dimensional array of complex numbers.

F(t) = A(t)( \cos \phi (t) + i \sin \phi (t))

Because it was

A(t) = | F(t) |, \quad \phi(t) = \arg F(t)

You can calculate the amplitude and initial phase of the circular motion with.


import math

def complex_to_wave(F):
    results = []
    for i, temp_complex in enumerate(F):
        amp = abs(temp_complex) # A(t)Equivalent to
        freq = i #Equivalent to t
        phase = math.atan2(temp_complex.imag, temp_complex.real) # phi(t)Equivalent to
        results.append((amp, freq, phase))
    return sorted(results, reverse = True)

Diameter A (t) = amp Frequency t/N = freq/N Initial phase phi (t) = phase If you add up the circular motions of

f(x)= \sum_{t=0}^{N-1}A(t)\cos (\phi (t) +\frac{2\pi t x}{N})

You can reproduce the coordinate values.

Draw a picture with Pillow

With pillow, you can draw shapes with python. Please check the detailed usage.

That's what we're doing in the implementation below.

Subtract the average value of each from the one-dimensional array of x and y coordinates, and adjust the scale appropriately. ↓ Fast Fourier transform of each array is performed to obtain a one-dimensional array of complex numbers. ↓ From the complex number, calculate the amplitude, initial phase, and frequency of the overlapping circular motions. ↓ Draw in order from the circle with the largest amplitude.


import math
import numpy as np
from PIL import Image, ImageDraw
from tqdm import tqdm_notebook as tqdm

def draw(draw_info, save_name):

    original_points = np.array(draw_info['points'])

    epicycles_x = complex_to_wave(np.fft.fft(original_points[:,0])/len(original_points[:,0]))
    epicycles_y = complex_to_wave(np.fft.fft(original_points[:,1])/len(original_points[:,1]))

    width = 640
    height = 360
    center = (width // 2, height // 2)
    frames = len(epicycles_x)
    frame_interval = 5
    center_of_horizontal_circle = (400, 70)
    center_of_vertical_circle = (70, 200)

    images = []
    drawing_points = original_points + np.array([center_of_horizontal_circle[0],

    for frame in tqdm(range(0, frames, frame_interval)):
        im = Image.new('RGB', (width, height), draw_info['background_color'])
        draw = ImageDraw.Draw(im)
        top_of_horizontal_circle = draw_epicycles(draw, epicycles_x, center_of_horizontal_circle,
                                                     frame, frames, vertical = False)
        top_of_vertical_circle = draw_epicycles(draw, epicycles_y, center_of_vertical_circle,
                                                   frame, frames, vertical = True)

        drawing_point = (top_of_horizontal_circle[0], top_of_vertical_circle[1])
        draw.line((top_of_horizontal_circle, drawing_point), fill = (0, 0, 0), width = 2)
        draw.line((top_of_vertical_circle, drawing_point), fill = (0, 0, 0), width = 2)
        draw_lines(draw, drawing_points[:frame], draw_info['colors'], draw_info['width'])

                   save_all = True, append_images = images[1:],
                   optimizer = False, duration = 17, loop = 0)

def draw_lines(draw, drawing_points, drawing_colors, drawing_width):
    drawing_width = int(drawing_width)
    drawing_points = tuple([tuple(i) for i in drawing_points])
    length = len(drawing_points)
    for i in range(length - 1):
        draw.line((drawing_points[i], drawing_points[i + 1]),
                  fill = drawing_colors[i + 1], width = drawing_width)
def draw_epicycles(draw, epicycles, center_of_circle, frame, frames, vertical):
    if vertical:
        init_theta = math.pi / 2
        init_theta = 0
    for amp, freq, phase in epicycles:

        theta = 2*math.pi*freq*frame/frames + phase + init_theta
        center_of_next_circle = (center_of_circle[0] + amp * math.cos(theta),
                                 center_of_circle[1] + amp * math.sin(theta))

        draw.ellipse((center_of_circle[0] - amp, center_of_circle[1] - amp,
                      center_of_circle[0] + amp, center_of_circle[1] + amp),
                     outline = (127, 127, 127))
        draw.line((center_of_circle, center_of_next_circle),
                  fill = (127, 127, 127), width = 2)
        center_of_circle = center_of_next_circle
    return center_of_circle


draw_info = {'points' : (points - np.mean(points, axis=0))/13000,
             'width' : 3,
             'colors' : [(0, 190, 255) for _ in range(len(points))],
             'background_color' : (255, 255, 255)}
draw(draw_info, 'save_name')



GIF tends to have a large file size, so please convert it to mp4.


You can also do this.


Recommended Posts

Drawing with Bezier curve and Fourier transform
Image processing from scratch with python (5) Fourier transform
I wanted to play with the Bezier curve
Wavelet transform of images with PyWavelets and OpenCV
FFT (Fast Fourier Transform): Formulas and Implementation Examples for Implementation
[# 2] Make Minecraft with Python. ~ Model drawing and player implementation ~
Calculate Pose and Transform differences in Python with ROS
[Python] One-liner FFT (Fast Fourier Transform) and other madness
Real-time drawing with matplotlib
Drawing with Python Tinker
Qiskit: Quantum Fourier Transform
With and without WSGI