确定一个点是否属于 n 阶科赫雪花指定的区域

Question

我正在尝试编写一个执行以下计算的 python 脚本：

输入： (1) List L：一些二维点的列表 (2) List V：三角形的顶点 (3) 正整数n：从那个三角形创建的科赫雪花的顺序

输出： 列表 O，L 的子集，包含 L 中位于区域 Kn 上或内部的点，该区域由 n 阶雪花定义。

---

我的尝试： 首先，我想我会先实现一个标准算法来绘制给定顺序（和边长）的雪花。这是我写的代码：

import turtle
from test import test

world= turtle.Screen()
t= turtle.Turtle()

def koch(t, order, size):
    if order == 0:
        t.forward(size)
    else:
        for angle in [60, -120, 60, 0]:
           koch(t, order-1, size/3)
           t.left(angle)

def koch_fractal(t, order, size, main_polygon_sides= 3):
    for i in range(main_polygon_sides):
        koch(t, order, size)
        t.right(360/main_polygon_sides)

koch_fractal(t, 2, 100)
world.mainloop()

但由于它没有说明雪花区域，我无法继续前进。接下来，我认为雪花的区域可能会有一些见解，所以我写了这个函数：

from math import sqrt
koch_cache={}
def koch_fractal_area(n, side):
    original_area = (sqrt(3)/4) * side**2 #Area of the original triangle 
    koch_cache[0] = original_area
    for i in range(n+1):
        if i not in koch_cache:
         koch_cache[i] = koch_cache[i-1] + (3*4**(i-1))*(sqrt(3)/4) * (side/(3**i))**2
    return koch_cache[n]

它实现了一个明确的公式来计算面积。同样，这似乎与我想要做的事情无关。

我该如何解决这个问题？ 提前致谢！

Answer 1

可以以与创建 Koch 雪花相同的方式递归地检查点位置。步骤是：

• 检查是给定三角形内的点，
• 如果不是，则该点位于某些三角形边的负侧。对于点位于负侧的每个边缘，递归检查该侧的“中间三角形”中的点，如果不是递归检查下两个可能的雪花边缘部分。

这种方法更快，因为它不会创建整个多边形并对其进行检查。

这里是使用 numpy 作为积分的实现：

import numpy

def on_negative_side(p, v1, v2):
    d = v2 - v1
    return numpy.dot(numpy.array([-d[1], d[0]]), p - v1) < 0

def in_side(p, v1, v2, n):
    if n <= 0:
        return False
    d = v2 - v1
    l = numpy.linalg.norm(d)
    s = numpy.dot(d / l, p - v1)
    if s < 0 or s > l:  # No need for a check if point is outside edge 'boundaries'
        return False
    # Yves's check
    nd = numpy.array([-d[1], d[0]])
    m_v = nd * numpy.sqrt(3) / 6
    if numpy.dot(nd / l, v1 - p) > numpy.linalg.norm(m_v):
        return False
    # Create next points
    p1 = v1 + d/3
    p2 = v1 + d/2 - m_v
    p3 = v1 + 2*d/3
    # Check with two inner edges
    if on_negative_side(p, p1, p2):
        return in_side(p, v1, p1, n-1) or in_side(p, p1, p2, n-1)
    if on_negative_side(p, p2, p3):
        return in_side(p, p2, p3, n-1) or in_side(p, p3, v2, n-1)
    return True

def _in_koch(p, V, n):
    V_next = numpy.concatenate((V[1:], V[:1]))
    return all(not on_negative_side(p, v1, v2) or in_side(p, v1, v2, n)
        for v1, v2 in zip(V, V_next))

def in_koch(L, V, n):
    # Triangle points (V) are positive oriented
    return [p for p in L if _in_koch(p, V, n)]

L = numpy.array([(16, -16), (90, 90), (40, -40), (40, -95), (50, 10), (40, 15)])
V = numpy.array([(0, 0), (50, -50*numpy.sqrt(3)), (100, 0)])
for n in xrange(3):
    print n, in_koch(L, V, n)
print in_koch(L, V, 100)

Answer 2

为了提高效率，当您将点与边进行比较时，请使用以下规则：

• 如果你在蓝色区域，点在外面，

• 如果你在橙色区域，点在里面，

• 否则您将需要递归测试，但请确保选择该点所在的绿色三角形，以便仅在一个子侧进行递归。

这可能看起来很小，但可以节省大量资金。事实上，在n-th 代，薄片有3 x 4^n 边（即第十代的3145728）；如果你递归到一个子方，你将只做12 测试！

@cdlane 的版本是最差的，因为它每次都会执行详尽的测试。 @ante 的版本介于两者之间，因为它有时会提前停止，但仍然可以执行指数级的测试。

---

一种简单的实现方法是假设要检查的一侧始终是(0,0)-(1,0)。然后测试测试点属于哪个三角形是一件简单的事情，因为顶点的坐标是固定的并且是已知的。这可以通过与直线进行四次比较来完成。

当您需要递归到子边时，您将通过将子边移动到原点、缩放 3 和旋转 60°（如果需要）来变换子边；对测试点应用相同的变换。

Answer 3

查找具有a routine for performing the "point in polygon" inclusion test 的Python 模块；使用turtle 的begin_poly()、end_poly() 和get_poly() 捕获代码生成的顶点，然后应用缠绕数测试：

from turtle import Turtle, Screen
from point_in_polygon import wn_PnPoly

points = [(16, -16), (90, 90), (40, -40), (40, -95)]

screen = Screen()
yertle = Turtle()
yertle.speed("fastest")

def koch(turtle, order, size):
    if order == 0:
        turtle.forward(size)
    else:
        for angle in [60, -120, 60, 0]:
            koch(turtle, order - 1, size / 3)
            turtle.left(angle)

def koch_fractal(turtle, order, size, main_polygon_sides=3):
    for _ in range(main_polygon_sides):
        koch(turtle, order, size)
        turtle.right(360 / main_polygon_sides)

yertle.begin_poly()
koch_fractal(yertle, 2, 100)
yertle.end_poly()

polygon = yertle.get_poly()

yertle.penup()

inside_points = []

for n, point in enumerate(points):
    yertle.goto(point)
    yertle.write(str(n), align="center")

    winding_number = wn_PnPoly(point, polygon)

    if winding_number:
        print(n, "is inside snowflake")
        inside_points.append(point)
    else:
        print(n, "is outside snowflake")

print(inside_points)

yertle.hideturtle()

screen.exitonclick()

% python3 test.py
0 is inside snowflake
1 is outside snowflake
2 is inside snowflake
3 is outside snowflake
[(16, -16), (40, -40)]