Determine if a point belongs to the region specified by a Koch snowflake of order n

Question

I am trying to write a python script that performs the following computation:

Input: (1) List L: a list of some 2-d points (2) List V: vertices of a triangle (3) Positive integer n: the order of Koch snowflake to be created from that triangle

Output: List O, a subset of L, that contains points from L that lies on or inside the region Kn, which is the region defined by the snowflake of order n.

---

My attempt: First, I thought I'd start by implementing a standard algorithm to draw a snowflake of a given order (and side length). Here's the code I wrote:

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

but since it doesn't say anything about the region of the snowflake, I could not move any further. Next, I thought the area of the snowflake might hold some insights, so I wrote this function:

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]

It implements an explicit formula to calculate the area. Again, it didn't seem to be related to what I was trying to do.

How can I approach this problem? Thanks in advance!

Answer 1

For efficiency, when you compare the point to a side, use the following rules:

• if you are in the blue region, the point is outside,

• if you are in the orange region, the point is inside,

• otherwise you will need to test recursively, but make sure to select the green triangle where the point is, so that you recurse on only one sub-side.

This may look like a minor difference, but it allows huge savings. Indeed, at the n-th generation, the flake has 3 x 4^n sides (i.e. 3145728 at the tenth generation); if you recurse to a single sub-side, you'll be doing 12 tests only !

The version by @cdlane is the worst, as it will perform an exhaustive test each time. The version by @ante is in between as it sometimes stops early, but can still perform an exponential number of tests.

---

An easy way to implement is to assume that the side to be checked against is always (0,0)-(1,0). Then testing to which triangle the test point belongs is a simple matter as the coordinates of the vertices are fixed and known. This can be done in four comparisons against a straight line.

When you need to recurse to a sub-side, you will transform that sub-side by moving it to the origin, scaling by 3 and rotating by 60° (if necessary); apply the same transforms to the test point.

Answer 2

It is possible to check point location in a same way how Koch snowflake is created, recursively. Steps are:

• Check is point inside given triangle,
• If not, than point is on negative side of some of triangle edges. For each edge where point is on negative side, check recursively is point in 'middle triangle' of that side, if not than check recursively next two possible snowflake edge parts.

This approach is faster since it doesn't create whole polygon and checks against it.

Here is implementation using numpy for points:

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 3

Find a Python module that has a routine for performing the "point in polygon" inclusion test; use turtle's begin_poly(), end_poly(), and get_poly() to capture the vertices that your code generates and then apply the winding number test:

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