[PYTHON] Generate that shape of the bottom of a PET bottle

I thought about programming that pentagonal one at the bottom of a PET bottle of carbonated drink, which everyone has always seen.

From the conclusion, about 70 points were made by self-evaluation. If you are worried about the result, scroll to the bottom first to see the completed image, and then decide whether to read further.

The shape of the bottom of a PET bottle of carbonated beverage

This is a pentagonal one. pet01.jpg

This shape is commonly referred to as ** Petaloid **. It means "petal-oid". It is so called because it looks like a petal when viewed from below.

If you have a carbonated PET bottle at home, take a look at the bottom.

Why is this shape?

PET bottles come in many shapes, but basically only carbonated drinks have this bottom shape. Since the pressure from the inside of carbonated drinks is strong, the bottom will swell if it is shaped like a tea bottle. It is said that petaloid is used to maintain the strength because it causes problems such as rolling and bursting as it is.

I would like to introduce an experiment in which baking soda was actually put into a PET bottle of tea to test its strength.

There is a deep secret in the shape of PET bottles-Diamond Online

Draw a petaloid

Target

The goal of this article is to programmatically generate petaloid shapes. We aim for quality as it is (don't expect too much because my judgment is more sloppy than the general public).

Search for materials by google

I wondered if it could be easily expressed with a single formula, but I did not get any results. Although pressure test data will come out ...

Maybe I'm designing with CAD.

Observe

If you can't rely on the internet, you have to do it yourself. Fortunately, I have as many real things as I have. First, use a cutter to cut out only the bottom of the PET bottle.

pet02.jpg

In addition, try halving vertically from the center. From the 12 o'clock direction to the 6 o'clock direction in the image above.

pet03.jpg

What i found

--Five of the same shape are lined up on the circumference --Valleys and mountains alternate -The valley part is a smooth curve from the center --The material is especially thick in the center

First of all, as you can see at a glance, five of the same shape are lined up on the circumference, so if you can reproduce one of them, you can see that you can achieve the goal just by copying and arranging five.

In terms of shape, there are alternating peaks and valleys, and the valleys are approximately 6 to 7 mm wide and appear to have lines running parallel from the center to the outer circumference. The shape of the mountain is quite complicated, but it looks like a fan when viewed from directly above (or directly below) and a trapezoid when viewed from the side. To be precise, the slope is made with a smooth curve between the peak and the valley.

Looking at the cross section, the cross section of the valley part is a smart curve that looks like a quarter of an elliptical curve. On the other hand, the cross section of the mountain part is a little complicated.

analyze.jpgcurve.jpg

Also, this is not very relevant to the article, but the bottom part was quite thick, about 3mm. At first I was wondering if I could cut it with scissors, but I had to use a thread saw because the blade wasn't toothy.

Set a policy

From this information, consider how to reproduce the shape.

First, as a basic policy, let's think about things based on the center of the bottle. Specifically, as shown in the figure below, we will try to reproduce the solid by calculating the cross-sectional shape in the direction of the azimuth angle θ each time and accumulating it with the center as the reference.

concept.jpg

Next, regarding the reproduction of the shape of the cross section, it is necessary to generate the following three different types of curves.

--Oval arc in the valley --Complex curves in the mountains --Complex curve of sloped part

Perhaps it could be summarized in one formula, but I can't think of it in my head, so this time I decided to go in the direction of using a spline curve (spline in case of trouble).

Reproduction of curves by splines

As mentioned above, let's consider how to generate three different curves based on one variable called azimuth $ θ $.

First, consider expressing the curve of the valley part with a spline function. As you can see, this curve is roughly in the shape of an arc on an ellipse, so it's not difficult to represent it with a spline.

For the sake of simplicity, we will use it as a perfect circle instead of an ellipse. If you want to return it to an ellipse, you can scale it vertically.

fig01.png

Draw a circle as shown above. Arrange $ N $ control points at equal intervals on the arc in the 4th quadrant (lower right) and draw a cubic spline curve connecting them. This curve is made to pass through all control points. It doesn't exactly match the arc, but it's a curve that is indistinguishable to the human eye. In addition, $ N = 9 $.

fig02.png

Next, I will make a curve for the mountain part. Here, I will process the curve of the valley and imitate the curve of the mountain as roughly observed.

Specifically, the control point $ (x_i, y_i) $ on the arc is extended by multiplying $ k_i $ from the center $ (0,0) $ of the circle ($ 1 \ leq i \ leq N $). Make a rough outline as shown in the figure below.

x'_i = x_i * k_i \\
y'_i = y_i * k_i

fig03_.png

While staring at the captured image, fine-tune the value of $ k_i $ so that it looks almost the same. The final parameter is

\begin{align}
K &= [k_1, k_2, k_3, ... , k_9]  \\
&= [1.00, 1.05, 1.18, 1.46, 1.36, 1.18, 1.07, 1.01, 1.00]
\end{align}

It was made. The end points (k1 and k9) must be 1.0 or the result will be strange.

Increasing the number of control points $ N $ to reproduce a more accurate curve may improve the quality, but it is quite difficult to adjust. Considering after it was completed, it might have been better to add one or two more points in order to improve the accuracy of the bottom part (fourth point from the left).

Spline interpolation of this scheme gives: I adjusted $ k_i $ while looking at the real PET bottle from the side, so it looks a lot like that.

fig04.png

Finally, make a curve of the slope between the peaks and valleys. This is done by blending the peak and valley curves with a variable called the deformation ratio $ D $ (what we are doing is just linear interpolation. % B7% 9A% E5% BD% A2% E8% A3% 9C% E9% 96% 93)).

\begin{align}
x''_i &= x_i * D + x'_i * (1-D) \\
&= x_i * (D + k_i - Dk_i) \\
y''_i &= y_i * D + y'_i * (1-D) \\
&= y_i * (D+k_i - Dk_i)
\end{align}
\quad (0 \leq D \leq 1)

fig05.png

The image above shows the curve of the slope from 0.0 to 1.0 in 0.2 steps. D = 0 is the reddish purple line and D = 1 is the yellowish green line. You can create a curve that smoothly fills the gap between peaks and valleys, such as contour lines.

Relationship between azimuth θ and deformation ratio D

When deciding the policy, I decided to make a curve of the cross section for each of multiple azimuth angles $ θ $. Also, since we introduced a variable called deformation ratio $ D $ to make a curve, we need to find the value of $ D $ based on the value of $ θ $.

I will draw a graph of θ (0 to 360 degrees) on the horizontal axis and D (0 to 1) on the vertical axis to clarify this relationship. First, for the top of the mountain, set $ D = 1 $. The range of the angle considered to be the top is 22 degrees from the visual measurement (roughly).

graphA02.png

Also, set $ D = 0 $ for the bottom of the valley. This is about 10 degrees. Since it was difficult to understand because it did not appear in the graph, the lower limit of the Y axis is set to -0.2.

graphA01.png

The white part of the gap between orange and blue is the slope between the mountain and the valley. This part will be filled between 0 and 1 with a curve, but how should we fill it? It's cool to say that it's mathematically like this, but I'm not that smart. I'm going to use the easing function here to fill it in nicely.

Reference: Easing Function Cheat Sheet

After trying various things, I chose ʻeaseInOutQuad` because it looked good.

\begin{align}
y &= easeInOutQuad(x)  \\
&= 
\begin{cases}
2x^2  &(0 \leq x < \frac{1}{2}) \\
-2x^2 + 4x - 1 &(\frac{1}{2} \leq x \leq 1) 
\end{cases}
\end{align}

The final graph looks like this: Note that the range on the vertical axis has been changed to 0 to 1.

graphA03.png

3D model generation

Now that we have the curves in all directions, we will use them to create a 3D model.

First, the azimuth angles θ are carved at equal intervals (for example, every 1 °) to generate a curve. For the curve, calculate the coordinates discretized at M points. Focus on the points on the two curves and select the three points to form a triangular polygon.

polygon.jpg

Complete!

You can draw directly with OpenGL etc. as it is, but this time I will export the file in DXF format and display it with 3D modeling software.

Since I wrote the code in Python this time, I used the ezdxf library that can handle DXF format in Python. The following image shows the output file displayed by software called Modo.

render1.jpg render2.jpg

Here is the result of assigning materials and rendering [^ 1].

render5.JPG

This is the result displayed in Blender. The AutoCAD DXF Import plug-in must be enabled.

render4.jpg

It looks like it looks good from the side, but it doesn't look very similar from the bottom. Maybe it's because the unevenness isn't enough, or because it's not rounded. Hmm ...

Source code

I don't think most people will read it, so I'll fold it up.

View Python source code

Requires NumPy, SciPy, ezdxf libraries.

petaloid.py


from typing import List, Tuple
import math
import numpy as np
from scipy.interpolate import splprep, splev
import ezdxf


#Definition of type hints
Point2D = Tuple[float, float]
Point3D = Tuple[float, float, float]


def mix(a: float, b: float, x: float) -> float:
    """Linear internal division x=When 0 a, x=When 1 b"""
    return b * x + a * (1.0 - x)


def easeInOutQuad(x: float) -> float:
    """Easing function"""
    if x < 0.5:
        return 2.0 * x * x
    else:
        return -2 * x * x + 4 * x - 1


def calcDeformRate(t: float) -> float:
    """
Find the curve deformation ratio D corresponding to the azimuth angle t
    :param t:Azimuth(The unit is degree).. When the remainder divided by 72 is 0, it corresponds to the valley, and when it is 36 °, it corresponds to the convex part.
    :return:Deformation ratio D(0~1)
    """
    t = math.radians(math.fmod(t, 72.0))
    M = math.radians(72.0)  #Angle range for one protrusion
    k = math.radians(5.0)  #Angle range of the valley
    tL = math.radians(25.0) - k  #One side boundary of the mountain
    tR = (M - k) - tL  #The other side of the mountain

    if t < k or t > M - k:
        return 0.0  #Valley part
    elif t < k + tL:
        return easeInOutQuad((t-k) / tL)  #Inclined
    elif t < tR:
        return 1.0  #Mountain part
    else:
        return easeInOutQuad(((M-k) - t) / tL)  #Inclined(Reverse)


def curve(angle: float, div_num: int) -> List[Point2D]:
    """
Calculate the outer curve according to the specified azimuth
    :param angle:Azimuth(0~360°)
    :param div_num:The number of points that make up the curve. The more it is, the smoother it is
    :return:A sequence of points representing a curve
    """

    # degree -> radian
    angle = math.fmod(angle, 360.0 / 5)  #Since it is a pentagon, 360 degrees is divided into 5 equal parts.

    #Find the deformation ratio D corresponding to the azimuth(D in the valley=0, D in the mountains=1)
    D = calcDeformRate(angle)

    #Find the deformation ratio of each control point
    #In the valley part, all control points are on the left(1.0)And it becomes an arc
    #In the mountain part, all control points draw a specified curve with the ratio on the right
    #Other parts(Mountain slope)Now draws an interpolated curve according to the value of D
    r_rate = [
        mix(1.00, 1.00, D),
        mix(1.00, 1.05, D),
        mix(1.00, 1.18, D),
        mix(1.00, 1.46, D),
        mix(1.00, 1.36, D),
        mix(1.00, 1.18, D),
        mix(1.00, 1.07, D),
        mix(1.00, 1.01, D),
        mix(1.00, 1.00, D)]
    #Find the position of the control point
    xs = []
    ys = []
    n = len(r_rate)
    for i in range(n):
        deg = i * 90.0 / (n - 1)
        rad = math.radians(deg)
        x = r_rate[i] * math.sin(rad)
        y = r_rate[i] * math.cos(rad)
        xs.append(x)
        ys.append(y)
    #Draw a cubic spline curve to pass through all control points
    spline = splprep([xs, ys], s=0, k=3)[0]
    detail = np.linspace(0, 1, num=div_num, endpoint=True)
    ix, iy = splev(detail, spline)
    points = list(zip(ix, iy))
    return points


def main():
    #Make a small protrusion on the center bottom
    num_bottom_point = 8  #Score of microprojections at the center bottom
    bottom_width = 0.2  #Center bottom size
    bottom_height = 0.05  #Central bottom height
    bottom_points = []
    for i in range(num_bottom_point):
        n = i / num_bottom_point
        angle = n * math.pi
        x = n * bottom_width
        #Y with cos=[1.0,1.0+height]Make a curve
        y = 1.0 + (math.cos(angle) + 1) / 2 * bottom_height
        bottom_points.append((x, y))

    #Create a 3D vertex list representing petaloid shapes
    num_curve: int = 180  #Number of azimuth divisions(360 every time,180 is every 2 times)
    num_point: int = 50  #Approximate number of curves in one direction(Including the center point)
    num_point_curve = num_point - num_bottom_point  # (Excluding the center point)
    aspect_adjust: float = 0.75  #Aspect ratio adjustment value
    vertices: List[Point3D] = []
    for i in range(num_curve):
        #Find the curve for each direction
        theta_deg = 360.0 / num_curve * i
        theta_rad = math.radians(theta_deg)
        points1 = curve(theta_deg, num_point_curve)

        #Connect with the small protrusion in the center
        points2: List[Point2D] = [*bottom_points]
        for p in points1:
            #Horizontal curve(X)Crush a little in the direction and make the vacant part a flat surface(What a delicate process ...)
            x = bottom_width + p[0] * (1.0 - bottom_width)
            y = p[1]
            points2.append((x, y))

        #Convert to 3D coordinates
        c = math.cos(theta_rad)
        s = math.sin(theta_rad)
        for p in points2:
            x2 = p[0]
            y2 = p[1]
            x3 = x2 * c
            y3 = y2 * aspect_adjust  #Adjust aspect ratio
            z3 = x2 * s
            vertices.append((x3, y3, z3))

    #Divide into polygons
    indices: List[Tuple[int, int, int]] = []
    for angle in range(num_curve):
        a0 = angle
        a1 = (angle + 1) % num_curve
        for i in range(num_point - 1):
            i00 = a0 * num_point + i
            i01 = a0 * num_point + i + 1
            i10 = a1 * num_point + i
            i11 = a1 * num_point + i + 1
            if i != 0:  #Only the center is a polygon, so exclude it so that it does not overlap
                indices.append((i00, i10, i11))
            indices.append((i11, i01, i00))

    #Output with dxf
    #I divided it into vertex and index, but since dxf's 3D face has only vertices,
    #Synthesize a list of vertices
    triangles: List[Tuple[Point3D, Point3D, Point3D]] = []
    for i in indices:
        v1 = vertices[i[0]]
        v2 = vertices[i[1]]
        v3 = vertices[i[2]]
        triangles.append((v3, v2, v1))
    doc = ezdxf.new('R2010')  # MESH requires DXF R2000 or later
    msp = doc.modelspace()
    for t in triangles:
        msp.add_3dface(t)
    doc.saveas("petaloid.dxf")


if __name__ == '__main__':
    main()

The remaining challenges

I'm not motivated anymore. If you have a strange person who wants to challenge a more complete petaloid, please refer to it.

--Reproduction of a more realistic constriction ――When you look at the real thing, both ends of the valley part look more constricted. In particular, the valley on the center side has a U-shaped constriction, but this time I gave up thinking that the calculation would be too complicated to complete. I don't know what to do specifically, but it would be nice to have an element that changes not only with the azimuth but also with the distance from the center or the zenith angle. --Spline-independent shape creation ――This time, I relied on splines (a lot of curved lines), but I wonder if it can be expressed with one simpler function. I think that the shape of the mountain is like a quartic function, but [Curve fitting](https://ja.wikipedia.org/wiki/%E6%9B%B2%E7%B7%9A%E3 % 81% 82% E3% 81% A6% E3% 81% AF% E3% 82% 81) has not been tried. --Design differences between manufacturers ――As far as I checked, it seems that there is no standard for PET bottles. The caliber of the cap seems to have a de facto standard and is standardized, but the shape of the bottle itself seems to depend on the manufacturer. The PET bottle I used this time is from Wilkinson, so it may be different from the one you have. However, the bottles in my house all seem to have the same shape regardless of the beverage manufacturer (maybe they are from the same container manufacturer, is it an oligopolistic market ...).

At the end

Q. Why only the bottom? What about the torso? A. I'm exhausted \ _ (: 3 "∠) \ _

[^ 1]: After importing DXF, Mesh CleanUp is executed to combine polygons. In addition, a subdivision surface is used to smooth the shape.

Recommended Posts

Generate that shape of the bottom of a PET bottle
A story that reduces the effort of operation / maintenance
[Python] A program that counts the number of valleys
Make a BOT that shortens the URL of Discord
# Function that returns the character code of a string
A story that analyzed the delivery of Nico Nama.
[Python] A program that compares the positions of kangaroos.
Create a shape on the trajectory of an object
A tool that automatically turns the gacha of a social game
A Python script that compares the contents of two directories
A memo that reproduces the slide show (gadget) of Windows 7 on Windows 10.
When incrementing the value of a key that does not exist
pandas Fetch the name of a column that contains a specific character
A formula that simply calculates the age from the date of birth
A story that struggled to handle the Python package of PocketSphinx
Distinguishing the agari shape of mahjong
A function that measures the processing time of a method in python
The story of creating a site that lists the release dates of books
I made a slack bot that notifies me of the temperature
Generate a list of consecutive characters
[python] A note that started to understand the behavior of matplotlib.pyplot
The story of writing a program
The story of making a module that skips mail with python
[Python] A program that rotates the contents of the list to the left
A story that visualizes the present of Qiita with Qiita API + Elasticsearch + Kibana
[Python] A program that calculates the number of chocolate segments that meet the conditions
I made a calendar that automatically updates the distribution schedule of Vtuber
[Python] A program that calculates the number of socks to be paired
Generate a list packed with the number of days in the current month.
Semi-automatically generate a description of the package to be registered on PyPI
I wrote a corpus reader that reads the results of MeCab analysis
The story of developing a web application that automatically generates catchphrases [MeCab]
Hit a method of a class instance with the Python Bottle Web API
How to create a wrapper that preserves the signature of the function to wrap
A note about the functions of the Linux standard library that handles time
The story of making a package that speeds up the operation of Juman (Juman ++) & KNP
From a book that makes the programmer's way of thinking interesting (Python)
[Python] Note: A self-made function that finds the area of the normal distribution
This and that of the inclusion notation.
A class that hits the DMM API
Measure the relevance strength of a crosstab
A quick overview of the Linux kernel
A memo explaining the axis specification of axis
Get the filename of a directory (glob)
The story of blackjack A processing (python)
A code that corrects the yoon / sokuon (sokuon)
Notice the completion of a time-consuming command
[Python] A program that rounds the score
Arrange the numbers in a spiral shape
A collection of Numpy, Pandas Tips that are often used in the field
The story of IPv6 address that I want to keep at a minimum
Let's do clustering that gives a nice bird's-eye view of the text dataset
The story of making a web application that records extensive reading with Django
The story of Django creating a library that might be a little more useful
[Discode Bot] I created a bot that tells me the race value of Pokemon
Implementation of a model that predicts the exchange rate (dollar-yen rate) by machine learning
The story of making a Line Bot that tells us the schedule of competitive programming
How to make a Raspberry Pi that speaks the tweets of the specified user
I made a github action that notifies Slack of the visual regression test
A model that identifies the guitar with fast.ai
How to calculate the volatility of a brand