Koch Curve algorithm in Python without using Turtle/Logo logic

Question

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

Answer 1

#!/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.