python - 确定一个点是否属于n阶Koch雪花指定的区域

我试图编写一个python脚本来执行以下计算：
输入：
（1）列表L：一些二维点的列表
（2）列表V:三角形的顶点
（3）正整数n：从该三角形创建koch雪花的顺序
输出：
列表o是l的一个子集，它包含位于区域kn上或区域kn内的来自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]

它实现了计算面积的显式公式。再说一次，这似乎和我想做的事无关。
我怎样才能解决这个问题？
提前谢谢！

最佳答案

可以用同样的方式检查点的位置，递归地创建koch雪花。步骤如下：
检查给定三角形内的点，
如果不是，则than点位于某些三角形边的负边上对于点位于负边上的每条边，递归地检查该边的“中间三角形”中的点，如果不超过递归地检查下两个可能的雪花边部分。
这种方法速度更快，因为它不创建整个多边形并检查它。
下面是使用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)