Koch Curve algorithm in Python without using Turtle/Logo logic

Inspired by the motivations expressed in this question.

I created some code to generate the point list (and also display a pretty picture) using just analytic geometry (no Logo/Turtle, no trigonometry). Also, instead of recursively inserting new points to a growing list, I calculated the final length, and then assigned values to the already initialized positions.

I would like to know if there is some obvious improvement to make it faster, cleaner, more pythonic and/or more canonical.

```
#!/bin/env python
# coding: utf-8

def kochenize(a,b):
HFACTOR = (3**0.5)/6
dx = b[0] - a[0]
dy = b[1] - a[1]
mid = ( (a[0]+b[0])/2, (a[1]+b[1])/2 )
p1 = ( a[0]+dx/3, a[1]+dy/3 )
p3 = ( b[0]-dx/3, b[1]-dy/3 )
p2 = ( mid[0]-dyHFACTOR, mid[1]+dxHFACTOR )
return p1, p2, p3

def koch(steps, width):
arraysize = 4**steps + 1
points = [(0.0,0.0)]*arraysize
points[0] = (-width/2., 0.0)
points[-1] = (width/2., 0.0)
stepwidth = arraysize - 1
for n in xrange(steps):
segment = (arraysize-1)/stepwidth
for s in xrange(segment):
st = s*stepwidth
a = (points[st][0], points[st][1])
b = (points[st+stepwidth][0], points[st+stepwidth][1])
index1 = st + (stepwidth)/4
index2 = st + (stepwidth)/2
index3 = st + ((stepwidth)/4)*3
result = kochenize(a,b)
points[index1], points[index2], points[index3] = result
stepwidth /= 4
return points

if __name__ == '__main__':
TOTALWIDTH = 1000.
points = koch(7, TOTALWIDTH)

# If you want a pretty picture (and have Cairo and PIL modules installed)
if True:
import cairo, Image

width = int(TOTALWIDTH)
height = int(TOTALWIDTH*0.32)
surface = cairo.ImageSurface(cairo.FORMAT_ARGB32, width, height)
cr = cairo.Context(surface)
cr.set_source_rgb(1,1,1)
cr.rectangle(0, 0, width, height)
cr.

#!/bin/env python
# coding: utf-8

def kochenize(a,b):
    HFACTOR = (3**0.5)/6

Stylistically, this would better as a global constant although there'll be some runtime cost to that.

dx = b[0] - a[0]
    dy = b[1] - a[1]
    mid = ( (a[0]+b[0])/2, (a[1]+b[1])/2 )
    p1 = ( a[0]+dx/3, a[1]+dy/3 )
    p3 = ( b[0]-dx/3, b[1]-dy/3 )
    p2 = ( mid[0]-dy*HFACTOR, mid[1]+dx*HFACTOR )

The [0] and [1] all over the place obscures your algorithm. Checkout collections.namedtuple

return p1, p2, p3

def koch(steps, width):
    arraysize = 4**steps + 1
    points = [(0.0,0.0)]*arraysize
    points[0] = (-width/2., 0.0)
    points[-1] = (width/2., 0.0)
    stepwidth = arraysize - 1
    for n in xrange(steps):
        segment = (arraysize-1)/stepwidth

This loop is slightly confusing to read. To understand the logic of what stepwidth, n and steps are all doing I have to look for it in several places. I'd probably write a generator which produces the decreasing values of stepwidth to encapsulate that

for s in xrange(segment):
            st = s*stepwidth
            a = (points[st][0], points[st][1])
            b = (points[st+stepwidth][0], points[st+stepwidth][1])
            index1 = st + (stepwidth)/4
            index2 = st + (stepwidth)/2
            index3 = st + ((stepwidth)/4)*3

They way you've coded that, the symmetry between the index is non-obvious. The code would be clearer if you hand't reduce the fractions. Also, index1, index2, index3 are begging to be in a list index

result = kochenize(a,b)
            points[index1], points[index2], points[index3] = result            
        stepwidth /= 4
    return points

if __name__ == '__main__':

I recommend putting the contents of this inside a function called main which you call from here

TOTALWIDTH = 1000.
    points = koch(7, TOTALWIDTH)

    # If you want a pretty picture (and have Cairo and PIL modules installed)
    if True:

If you want this to be optional check command line arguments or check whether the necessary modules are installed. Putting if True as anything other then a quick hack is frowned upon.

import cairo, Image

        width = int(TOTALWIDTH)
        height = int(TOTALWIDTH*0.32)
        surface = cairo.ImageSurface(cairo.FORMAT_ARGB32, width, height)
        cr = cairo.Context(surface)

I dislike cr, it doesn't obviously stand for anything.

cr.set_source_rgb(1,1,1)
        cr.rectangle(0, 0, width, height)
        cr.fill()
        cr.translate(width*0.5, height*0.95)
        cr.scale(1, -1)

        cr.set_source_rgb(0,0,0)
        cr.set_line_width(0.5)
        cr.move_to(*points[0])
        for n in range(len(points)):
            cr.line_to(*points[n])
        cr.stroke()

        im = Image.frombuffer("RGBA", (width, height), surface.get_data(), "raw", "BGRA", 0,1)
        im.show()

If you want more speed, you should look into numpy. It will allow you to do vector operations which will be faster then what you can do in straight python.

StackExchange Code Review Q#4832, answer score: 2