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]-dy*HFACTOR, mid[1]+dx*HFACTOR )
    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.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()
#!/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.