如何从点的X,Y列表和偏移距离中获取偏移样条线的X,Y坐标

Question

我需要制作一个偏移平行的翼型剖面曲线，但是我不知道如何使所有点与原始主要剖面曲线上的点在所需距离上等距分布。 这是我的翼型剖面示例 这是我最好的但并不理想的方法 编辑 @Patrick 解决方案，距离为0.2

- efirvida

---

@IlyaPeterov 曲线以x，y坐标存储在np.array中，我不需要非常精确的结果，只需要近似曲线即可，因为我需要在这两个轮廓之间制作有限元网格。 - efirvida

Delta 应该从每个线段的中点开始计算。让我举个例子。 - Patrick Maupin

@PatrickMaupin 我非常感激您提供的示例。 - efirvida

我没有你的点，所以无法完全做到，但我已经从正弦波中偏移了一些。结果发现你还必须处理向量的方向。 - Patrick Maupin

@efirvida 是的，但不如Patrick的好。它没有半圆形，并且总体上有点不精确。虽然这对我来说是一个很好的练习，谢谢你。 - Ilya Peterov

Answer 1

你需要特殊处理无穷/零斜率的情况，但基本方法是使用插值来计算点处的斜率，然后找到垂直斜率，然后计算该距离上的点。我已经修改了这里的示例以添加第二个图表。它可以使用你提供的数据文件，但你可能需要更改不同包络的符号计算。根据你想要包络连续的评论，我在末尾添加了一个俗气的半圆形，它非常接近为你完成此操作。实际上，在创建包络时，你可以使其更加圆润和凸出，它将工作得更好。另外，你需要重叠起点和终点，否则会有间隙。此外，它几乎肯定可以更有效率--我绝对不是numpy专家，所以这只是纯Python。

def offset(coordinates, distance):
    coordinates = iter(coordinates)
    x1, y1 = coordinates.next()
    z = distance
    points = []
    for x2, y2 in coordinates:
        # tangential slope approximation
        try:
            slope = (y2 - y1) / (x2 - x1)
            # perpendicular slope
            pslope = -1/slope  # (might be 1/slope depending on direction of travel)
        except ZeroDivisionError:
            continue
        mid_x = (x1 + x2) / 2
        mid_y = (y1 + y2) / 2

        sign = ((pslope > 0) == (x1 > x2)) * 2 - 1

        # if z is the distance to your parallel curve,
        # then your delta-x and delta-y calculations are:
        #   z**2 = x**2 + y**2
        #   y = pslope * x
        #   z**2 = x**2 + (pslope * x)**2
        #   z**2 = x**2 + pslope**2 * x**2
        #   z**2 = (1 + pslope**2) * x**2
        #   z**2 / (1 + pslope**2) = x**2
        #   z / (1 + pslope**2)**0.5 = x

        delta_x = sign * z / ((1 + pslope**2)**0.5)
        delta_y = pslope * delta_x

        points.append((mid_x + delta_x, mid_y + delta_y))
        x1, y1 = x2, y2
    return points

def add_semicircle(x_origin, y_origin, radius, num_x = 50):
    points = []
    for index in range(num_x):
        x = radius * index / num_x
        y = (radius ** 2 - x ** 2) ** 0.5
        points.append((x, -y))
    points += [(x, -y) for x, y in reversed(points)]
    return [(x + x_origin, y + y_origin) for x, y in points]

def round_data(data):
    # Add infinitesimal rounding of the envelope
    assert data[-1] == data[0]
    x0, y0 = data[0]
    x1, y1 = data[1]
    xe, ye = data[-2]

    x = x0 - (x0 - x1) * .01
    y = y0 - (y0 - y1) * .01
    yn = (x - xe) / (x0 - xe) * (y0 - ye) + ye
    data[0] = x, y
    data[-1] = x, yn
    data.extend(add_semicircle(x, (y + yn) / 2, abs((y - yn) / 2)))
    del data[-18:]

from pylab import *

with open('ah79100c.dat', 'rb') as f:
    f.next()
    data = [[float(x) for x in line.split()] for line in f if line.strip()]

t = [x[0] for x in data]
s = [x[1] for x in data]


round_data(data)

parallel = offset(data, 0.1)
t2 = [x[0] for x in parallel]
s2 = [x[1] for x in parallel]

plot(t, s, 'g', t2, s2, 'b', lw=1)

title('Wing with envelope')
grid(True)

axes().set_aspect('equal', 'datalim')

savefig("test.png")
show()

- Patrick Maupin

---

在你的代码中，x2, x1 和 y2, y1 是翼型轮廓中的两个相邻点？因此我必须从 x[0] 迭代到 x[len(x)-1]？ - efirvida

是的，但这只会给你与原曲线相同数量的点数。你甚至可能想要插值更多的点来使你的另一个曲线更平滑，因为它比你的第一个曲线更长。 - Patrick Maupin

@efirvida 另外，您可能希望从线段的中点开始计算增量，因为在那里切线更加准确，例如 (x1 + x2) / 2, (y1 + y2) / 2。 - Patrick Maupin

现在的工作效果更好了，但是最后仍然存在问题。我知道如果没有数据来计算剩余的轮廓，就像我上一条评论中所说的那样，我能否通过拟合曲线来近似计算？ - efirvida

到达一个点并不太好。在创建信封时，我会在原始数据上稍微备份一下，并可能制作一个小的半圆。你能使它更圆润和凸出，它就会工作得更好。此外，您需要重叠开头和结尾，否则将会有间隙。 - Patrick Maupin

Answer 2

如果你愿意（并且能够）安装第三方工具，我强烈推荐使用Shapely模块。以下是一个小示例，可以同时向内和向外偏移：

from StringIO import StringIO
import matplotlib.pyplot as plt
import numpy as np
import requests
import shapely.geometry as shp

# Read the points    
AFURL = 'http://m-selig.ae.illinois.edu/ads/coord_seligFmt/ah79100c.dat'
afpts = np.loadtxt(StringIO(requests.get(AFURL).content), skiprows=1)

# Create a Polygon from the nx2 array in `afpts`
afpoly = shp.Polygon(afpts)

# Create offset airfoils, both inward and outward
poffafpoly = afpoly.buffer(0.03)  # Outward offset
noffafpoly = afpoly.buffer(-0.03)  # Inward offset

# Turn polygon points into numpy arrays for plotting
afpolypts = np.array(afpoly.exterior)
poffafpolypts = np.array(poffafpoly.exterior)
noffafpolypts = np.array(noffafpoly.exterior)

# Plot points
plt.plot(*afpolypts.T, color='black')
plt.plot(*poffafpolypts.T, color='red')
plt.plot(*noffafpolypts.T, color='green')
plt.axis('equal')
plt.show()

以下是输出结果；请注意内部偏移的“蝴蝶结”（自交）已经被自动删除：

- subnivean