[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.

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.

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

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.

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.

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.

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 $.

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

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.

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)

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

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.

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.

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.

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.

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

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

• Actually, the shape of the central part was a little distorted due to the influence of the spline function, so I used a cosine curve to make it look natural, but it does not affect the overall appearance, so I will omit the explanation. There is.

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.

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.