## Thursday, 6 August 2009

### Sims bakery: Steak and kidney pie

Another interesting pie that I bought with the lambs fry and bacon one, was a steak and kidney pie. On the outside it looks like any other potentially satisfying pie. But the oozing innards unfortunately resemble its sibling. There was too little meat substance and too much gravy. Unlike the gravy in the lambs fry pie, the steak and kidney gravy was not floury and pasty, and had some flavour. But the texture the predominance of gravy brought to the pie just didn't make it a good one. Sims also do an interesting sounding Venison pie, but at a whopping price of \$6.00NZ with what experience dictates will be a distinct lack of meatiness, I'm not springing for it.

### Mollweide projections

Distracted!

With working code to generate mollweide projections of my planetary surface and display it in a pyglet window, I got a little distracted. I ditched the window with its default width and height and generated images of a width and height with correct aspect ratio. Given there are 360 degrees of longitude and 180 degrees of latitude, the width needed to be twice the height. With this working, I randomly generated seeds for the noise and images of the corresponding projection for each seed. I kind of wanted to find a nice looking surface with many separate islands, and not just one or two planet wrapping islands.

A simple loop generates the images. I don't really need to use the Wichmann Hill generator that comes with Python, but I've been using it get to get small extractable seeds so that I can reseed the generator and get reproducible sequences. Irrelevantly, I am not illustrating this in the following code snippet.

`    rg = random.WichmannHill()    for i in range(2000):        v = int(rg.random() * sys.maxint)        terrain.set_noise_seed(v)        generate_mollweide_image("mollweide/water-%011d.jpg" % v, 1000, 500)`
Keep in mind that the projection attaches to the sides and not the top and bottom. So an island that features on the edge of the left or right hand side will wrap to the other side.

75893359 (like the majority of other seeds) produces a surface with one overly large island wrapping the lower half of the planet. 1320502026 produces a nicer planet with one larger island, several of medium size and a scattering of smaller ones. While 1898084921 produces one huge island, it is unique and interesting in having two huge inland bodies of water. Lakes, mountains and a range of interesting islands of different sizes are present in 2029406243. And there I am distracted into looking for interesting planetary surfaces again. At least it resulted in me optimising the code from taking 31 seconds to generate the height data of a given seed (640x480) to taking 1.8 seconds to generate the height data for a larger area (1000x500).
`    def generate_mollweide_image(fileName, width, height):        global heights, height_min, height_max        print "generating heights"        t = time.time()        xscale = width        yscale = int(height * 1.3)  # Cover more lines and leave no black ones.        heights = terrain.mollweide_projection(8, xscale, yscale, width, height)        height_min = min(heights, key=lambda x: sys.maxint if x is None else x)        height_max = max(heights, key=lambda x: -sys.maxint if x is None else x)        dt = time.time() - t        print "generated heights", dt, "seconds", (height_min, height_max)        print "generating colours"        t = time.time()        imageData = generate_image_data(width, height)        dt = time.time() - t        print "generated colours", dt, "seconds"        print "saving image"        mapImage = pyglet.image.create(width, height)        mapImage.set_data('RGB', width * 3, imageData.tostring())        mapImage.save(fileName)`
All of the height generation now happens in a C extension (`terrain.mollweide_projection`), which contributed to half of that optimisation.

Here's the original Python height generation routine:
`def generate_mollweide_heights_xy():    global heights, height_min, height_max    heights = [ None ] * (w.height * w.width)    xscale = w.width    yscale = int(w.height * 1.3)  # Cover more lines and leave no black ones.    longitude0 = 0.0    for yi in xrange(yscale):        v = float(yi) / yscale        # Stay within the domain of -pi/2 to pi/2        latitude = (v - 0.5) * math.pi        dtheta = 0.149        theta = latitude        pi_sin_lat = math.pi * math.sin(latitude)        i = 0        while math.fabs(dtheta) > 5e-8:            if i == 30:                break            dtheta = -(theta + math.sin(theta) - pi_sin_lat) / (1.0 + math.cos(theta))            theta += dtheta            i += 1        theta = theta / 2.0        for xi in xrange(xscale):            u = float(xi) / xscale            xj = (u - 0.5) * math.cos(theta)            yj = math.sin(theta)            xc = int((xj + 0.5) * w.width)            yc = int(((yj + 1.0) / 2.0) * w.height)            longitude = u * math.pi * 2.0 - longitude0            x = math.cos(longitude) * math.cos(latitude)            y = math.sin(longitude) * math.cos(latitude)            z = math.sin(latitude)            heights[yc * w.width + xc] = terrain.fBm(8, x, y, z)    filtered_heights = [ x for x in heights if x is not None ]     height_min = min(filtered_heights)    height_max = max(filtered_heights)`
It might be useful for someone who would also otherwise try and plug in the formulae present at Wikipedia or Mathworld.

I couldn't get the inverse transformation to produce anything that looked vaguely correct. And the forward transformation was much the same until I ditched a few elements and figured out the domain of the variables.

The corresponding C extension function:
`static PyObject *terrain_mollweide(PyObject *self, PyObject *args) { PyObject *list, *ob; int octaveCount; float u, v, du, dv; float longitude, latitude; float dtheta, theta, sin_lat, cos_lat, cos_theta; float x, y, z; float xf, yf, xiterations, yiterations; int i, width, height, arr_size, arr_idx; float *arr; if (!PyArg_Parse(args, "(iffii)", &octaveCount, &xiterations, &yiterations, &width, &height))  return NULL; arr_size = width * height; arr = malloc(arr_size * sizeof(float)); if (arr == NULL) {  PyErr_SetString(PyExc_RuntimeError, "failed to allocate working memory");  return NULL; } for (i = 0; i < arr_size; i++)  arr[i] = -1.0f; du = 1.0f / xiterations; dv = 1.0f / yiterations; v = 0.0f; while (v < 1.0f) {  latitude = (v - 0.5f) * F_PI;  theta = latitude;  sin_lat = sinf(latitude);  cos_lat = cosf(latitude);  i = 0;  dtheta = 0.149f;  while (fabs(dtheta) > 4.9e-08f) {   if (i == 30)    break;   dtheta = -(theta + sinf(theta) - F_PI * sin_lat) / (1.0f + cosf(theta));   theta += dtheta;   i += 1;  }  theta = theta / 2.0f;  yf = ((sinf(theta) + 1.0f) / 2.0f);  cos_theta = cosf(theta);  u = 0.0f;  while (u < 1.0f) {   xf = (u - 0.5f) * cos_theta + 0.5f;   longitude = u * F_PI * 2.0f;   x = cosf(longitude) * cos_lat;   y = sinf(longitude) * cos_lat;   z = sin_lat;   arr_idx = width * (int)(yf * height) + (int)(xf * width);   if (arr_idx < 0 && arr_idx >= arr_size) {    free(arr);    PyErr_SetString(PyExc_RuntimeError, "internal out of range error");    return NULL;   }   arr[arr_idx] = fBm(octaveCount, x, y, z);   u += du;  } /* while(u < 1.0f) */  v += dv; } /* while(v < 1.0f) */ list = PyList_New(width * height); if (list == NULL) {  free(arr);  return NULL; } for (i = 0; i < arr_size; i++) {  if (arr[i] == -1.0f) {   Py_INCREF(Py_None);   ob = Py_None;  } else {   ob = PyFloat_FromDouble(arr[i]);   if (ob == NULL) {    Py_DECREF(list);    free(arr);    return NULL;   }  }  PyList_SET_ITEM(list, i, ob); } free(arr); return list;}`
The image generation and colouring is now the slowest part of the process. Even worse is that the array and string shenanigans that pyglet insists on encounters `MemoryError` exceptions with larger sized projections. But really, neither of these things matter as the map projection is serving its purpose.

## Wednesday, 5 August 2009

### Map projections with Pyglet

In 2001, I wrote code to generate fractal-based landscapes and wanted to visualise what it produced using map projections. The projection I was most familiar with was the mercator projection, and while it showed that the generation was working, it wasn't ideal for visualisation due to the distortion present within it. So I played around with other projections and eventually settled on an orthographic projection. The end result was an animated gif of around a megabyte in size that showed a generated planet rotating. Unfortunately, at some point in time, the code was lost.

With a bit of time to spare to create some similar code, this time I decided on using the mollweide projection. Pygame has always been a little too complicated, sure there is documentation, but more often than not when I try and use a feature there will be an non-obvious little detail that prevents it from working correctly if at all. So this time, I decided to use Pyglet.

I've used Pyglet in the past, and had a simple script that used OpenGL to display a grid and move around it on the XZ plane. Unpromisingly, this script no longer worked and the grid wouldn't display. It turns out that Pyglet has changed the way their window events work. My resize event was being called:

`@window.eventdef on_resize(width, height):    # Override the default on_resize handler to create a 3D projection    glViewport(0, 0, width, height)    glMatrixMode(GL_PROJECTION)    glLoadIdentity()    gluPerspective(60.0, width / float(height), 0.1, 1000.0)    glMatrixMode(GL_MODELVIEW)`
But it turns out that now if you do not return `pyglet.event.EVENT_HANDLED` from this event, Pyglet decides to assert some default behaviour. Some helpful bloke from the Pyglet mailing list pointed me in the right direction with this.

During the period where I was unaware that this non-backwards compatible change had occurred I tried two different examples from the Pyglet documentation. Both errored due to bugs in the current release of Pyglet (1.1.3). It seems that Pyglet isn't in a stable state at the moment, and given that 1.1.3 is a maintenance release or something, I am not sure it ever is. I tracked down bug reports for the problems and noted they were considered to be fixed in the trunk, both at different points in time months before 1.1.3 was released. It would be nice to be able to see a list of all the known bugs with the current release (especially since any I encounter already seem to be fixed) so that I as a user do not have to waste my time unexpectedly experiencing them, tracking down the problem and trying to report them. The best approach with regard to Pyglet when I experience a bug seems to be ignoring the bug reporting system and just looking to see if a fix has been made in the trunk and manually patching it into my own code.

I mentioned something above about PyGame. I'd also like to say that sure there is documentation for Pyglet, but more often than not when I try and use a feature there will be an non-obvious little detail that prevents it from working correctly if at all.

In any case, eventually I reached the point where I had written code to do a mollweide projection and had it displaying in the PyGame window, looking the following way. I am also saving out the images I generate, in order to take not of different approaches and mistakes made during the process of development. By itself Pyglet can only save out PNG images, so this was what it saved out to a PNG file. Now one picture is a little bright and the other is a little dark. The PNG images saved by Pyglet have darker colours. I initially tried to compensate by continually brightening the map colours. The brighter image, which is identical to what is displayed by Pyglet, was saved out as a JPG image after I had installed PIL. If PIL is present, than Pyglet can handle a wider range of image formats.