A raytracer in python – part 3: samplers

Date

In the previous post, we explored a very basic way of plotting images: shooting a ray from the center of every pixel, and plot the color of the object we hit. The result is a rather flat, very jagged image

image

Border jagging arises from the fact that we are sampling with a discrete grid (our ViewPlane) an object that is smooth due to its functional expression (our spheres). The sampling is by nature approximated, and when mapped to pixels it produces either the color of the object, or the color of the background.

We would like to remove the jagged behavior. To achieve this, we could increase the resolution. That would make the jaggies less appreciable thanks to a higher number of pixels. An alternative is to keep the lower resolution, and shoot many rays per each pixels instead of only one (a technique called supersampling). To compute the final color, weighting is performed according to the number of rays that impact vs. the total number of rays sent from that pixel. This technique is known as anti-aliasing.

image

The figure details visually what said above: the real description of the sphere is smooth (left hand figure). Shooting one ray per pixel, and coloring the pixel according to hit/no-hit, produces either a fully colored pixel, or a fully background pixel (center). With supersampling, we shoot a higher number of rays, and per each pixel we perform weighting of the color.

As a result, the jaggies in the sphere are replaced with a smoother, more pleasant transition

image

Choice of Samplers

The pattern used for the supersampled rays is important. The module samplers in the python raytracer implements three common strategies.

Regular sampling

It is the most intuitive and easy, and the one used above: a regular grid of rays. It is easy to implement but it can introduce unpleasant artifacts for more complex situations. Plotting the position of the rays in the pixel will produce the following layout (for a 8x8 supersample grid)

image

As we see, the layout is regular on the grid of subcells (painted yellow and white for better visualization) that define the pixel. On the vertical and horizontal distributions (plotted on the grey bars above and on the left) we also see a regular distribution, but its regularity may introduce artifacts in some cases.

Random sampling

Random sampling positions the rays at random within the pixel. This may sound appealing, but it may instead end up as suboptimal: it can produce random clumping, in particular for a small number of samples. This will unbalance the weighting leading to an incorrect evaluation. Plotting one distribution one may obtain

image

Note the uneven distribution of the points, leaving large parts not sampled and other parts oversampled. In addition, the vertical and horizontal distribution tend to be uneven.

Jittered sampling

Jittered sampling takes the best of both worlds: the regularity of the Regular sampler with a degree of randomness from the Random sampler. The idea is to select the center of each subcell and apply randomization, so that each subcell produces only one ray, but without the artifact inducing regularity proper of the Regular sampler.

image

Computational cost impact

Unfortunately, performing supersampling makes the creation of the image considerably slower. In the following table you can see the timings (in seconds) for Jittered and Regular sampling, compared against the case with no sampling. The size of the image rendered is 200x200 pixels.

Sample size/sets No sampling Regular Random Jittered
1/1 23      
4/1   68 63 160
9/1   158 152 271
16/1   229 235 276
16/2   223 247 267
16/4   240 277 260

As we can see, supersampling introduces a considerably higher computational effort. We also see that having multiple sets (a requirement to prevent the same set of subsamples to be reused for adjacent pixels, something that again would introduce artifacts) does not really change the timings. According to the current implementation, I expect this to be verified. On the other hand, I don't expect timings for Regular and Jittered to be so different, since the creation of the values is performed once and for all at startup. This is worth investigating while looking for performance improvement. In the next post I will perform profiling of the python code and check possible strategies to reduce the timings, eventually revising the current design.

Current implementation

The current implementation of python-raytrace can be found at github. In this release, I added the samplers. Samplers are derived classes of the BaseSampler class, and are hosted in the samplers module. Derived Samplers must reimplement _generate_samples. Points are stored because we want to be able to select sets at random as well as replay the same points. Samplers also reimplements the __iter__() method as a generator of (x,y) tuples, with x and y being in the interval [0.0, 1.0). Once initialized, the Sampler can therefore be iterated over with a simple

for subpixel in sampler:
    # use subpixel

The World class can now be configured with different Samplers for the antialiasing. The default Sampler is a Regular 1 subpixel Sampler, which is the original one-ray-per-pixel sampling.


A raytracer in python - part 2: rendering multiple objects

Date

A quick addition needed to the raytracer is providing freedom to add more objects to the rendering scene. In Part 1, the design was such that only one object, a sphere, could be drawn. The new code allows much more flexibility. I added a Plane object, introduced assignment of colors to the objects, divided the source into multiple files, and fixed a bug relative to rendering direction. Let's see more specifically.

The class World is the main interface to user programming of the raytracer. This small program creates and renders a scene containing four spheres, a white one in the origin, the others along the axes, each with different colors.

import raytrace
from raytrace import objects

w=raytrace.World()
w.add_object(objects.Sphere(center=(0.0,0.0,0.0),
                            radius=10.0,
                            color=(1.0,1.0,1.0)
                           )
            )
w.add_object(objects.Sphere(center=(50.0,0.0,0.0),
                            radius=10.0,
                            color=(1.0,0.0,0.0)
                           )
            )
w.add_object(objects.Sphere(center=(0.0,50.0,0.0),
                            radius=10.0,
                            color=(0.0,1.0,0.0)
                           )
            )
w.add_object(objects.Sphere(center=(0.0,0.0,50.0),
                            radius=10.0,
                            color=(0.0,0.0,1.0)
                            )
            )
w.render()

The resulting image is the following

image

Clearly, the white sphere is not visible, as it's hidden by the blue sphere. This test allowed me to discover a problem with orientation: the green sphere was on the wrong side. I therefore had to analyze a bit the orientation and the different coordinate systems into play here

  1. The measure of the screen is given in pixels. We are used to this when it comes to screen size, for example. An image which is 320x200 means that it's 320 pixels wide (horizontal resolution) and 200 pixel high (vertical resolution). Don't fall into "thinking matrix", where NxM means N rows x M columns. It is the exact opposite.
  2. The geometry uses the cartesian system, which has the origin in the center of the picture. The x axis is oriented towards the right, the y axis towards the top, and the z axis towards the observer. This is not different from a traditional cartesian layout: for the z=0 plane, the top left pixel correspond to a (-,+) coordinate, the bottom right to a (+,-) coordinate. Points closer to the observer have positive z, honoring the right hand system. The camera in the above picture is at z=+100.0
  3. The raytracer uses pixel coordinates for the viewplane. Pixel 0,0 is at the bottom left. Changing the first index moves horizontally (along the row) from left to right. Changing the second index moves vertically (along the column) from bottom to top. As a consequence, the point at the bottom right is (hres-1,0), and at the top left is (o, vres-1). Note that this is equivalent to a cartesian (x,y) system with origin on the bottom left corner. Not a surprise, since there is a direct mapping between the pixels and the rays' origins.
  4. Finally, the pixel image indexing of pygame and PIL. For them, pixel 0,0 is top left. Like the case above, incrementing the first index also moves along the row from left to right. However, incrementing the second index moves vertically from top to bottom, which is the opposite of the raytracing index. Bottom left is (0, vres-1) and the bottom right is (hres-1, vres-1).

A remapping is therefore needed from the pixel coordinate of the rendering and the pixel coordinate of the display (e.g. pygame). The transformation is trivial, of course, but it must be kept into account, otherwise pics will be flipped horizontally.

Another interesting fact is that, according to "Ray Tracing from the ground up" it's useful to perform the raytracing operation starting from the bottom and working our way up, rendering pixels from left to right. According to them this is for coding symmetry and future convenience, so we stick to it.

Finding the foremost object

In the book, finding the foremost object is made more complex by the language used, C++. In python, you can use functional programming style to obtain the same in a very concise statement. The idea is to cast a ray, then go through all objects in the world to find the intersection point (if any) between the ray and the object. If more than one object is hit, the one with the intersection point closer to the observer commands the pixel color, since it's in front.

I achieve this with the following code

def hit_bare_bones_object(self,ray):
   def f(o):
     shadeRec = o.hit(ray)
     if shadeRec:
       return (shadeRec.parameter, o)
     else:
       return None

   try:
     foremost=sorted( \
                filter(lambda x: x is not None, \
                  map(f, self.objects)
                ), key=lambda x: x[0]
              )[0][1]
   except IndexError:
     return None

   return foremost

What does this code do ? I defined a simple internal function f which accepts an object, performs the hit and returns a tuple containing the hit point position (as a parameter of the ray, so it's a single number, not a xyz coordinate) and the object.

Now, I use this function to map all the objects defined in the world. I will obtain a list with one entry per each object, either a None (not hit) or a 2-tuple containing the parameter and the hit object. I filter out the None entries, leaving only the 2-tuples and then sort according to their first element. The 2-tuple with the lowest parameter is now at index 0, and the [1] element of this tuple is the foremost object. At any time, the list may be empty (such as if you don't have any object, or no object is hit. In that case, a IndexError will be raised and that will indicate that the ray hit nothing. I may rework on this function later on, but for this second round, it suits my needs.

It's now time to move on to samplers. Given that the code is growing in size, I created a git repository you can clone from. The release of this post is available here. The code is under BSD license.


A raytracer in python - part 1: basic functionality

Date

Some time ago I visited Sydney, and I made a tragic mistake: I entered the University bookshop. Why a mistake, you say? I am book maniac. As soon as I enter a book shop (live or on web) I end up spending up to a thousands euro every time. This time, it was not the credit card, but the limit of 20 kg on my luggage to put a limit on what I could buy. Needless to say I got really good stuff, in particular this 761 pages of pure awesomeness: "Ray tracing from the ground up". This book teaches you how to write a raytracer, step by step from a first basic skeleton to incredibly complex optic effects. It also provides code, available for download under GPL. I haven't downloaded it, but according to the book it's in C++. Since I want to keep myself fresh with python, I will write my own version in this language.

The first objection I may hear is performance. Native python is much slower than C++. I agree, and I do it on purpose. My objective is in fact not only to write a raytracer for fun (and the world is full of already awesome raytracers), but also to perform some infrequent python exercises: interfacing with C, parallelization and, hopefully, some OpenCL programming. This is my hope, at least. I will try to spend some time on it and see how far I can get.

What is raytracing and how does it work ?

Raytracing is a technique to produce a photorealistic image. It works by projecting rays from the observer to the scene, and coloring pixels on a viewplane for every ray that intersects an object.

image

This mechanism resembles how vision works, although in the opposite direction. Light rays from a lamp hit objects and their reflection happens to scatter around. Some of these reflections will enter our eyes and allow us to see the world around us. Doing the opposite, tracing from the observer, is clearly more efficient as we don't care about the rays not hitting the observer (at least in its basic implementation), only those who do.

Performing raytracing (or to be more accurate, for now just its basic form raycasting) involves the following steps:

  1. define a geometric object in space, like for example a sphere
  2. define a view panel made of pixels
  3. shoot one straight line (ray) from the center of each pixel
  4. if the ray intersects the object, mark the pixel colored, otherwise mark it with the background color

That's it. Basically, it's an exercise in geometry: finding intersections between lines and 3D objects in 3D space. This first program, ray.py, does exactly that. You will need to install the Python Imaging Library, pygame and numpy. The result is intriguing:

image

Ok, I have a very loose definition of "intriguing", but it's a start.


Export vim text (with colors) to HTML

Date Tags HTML / vim

Vim is a great, great programming tool. Even after years of experience with it you still get to discover, either by change or by sharing, fantastic tips to make an impossible task incredibly easy.

It is the case with my recent problem of exporting the visual aspect of vim (as from terminal) to an HTML document, including the highlight colors. How to do it? Well, I tried pygments and it's really great, but I asked on SuperUser for alternative methods. It appears that Vim supports this natively. Just issue

:TOhtml

to Vim, and it will create and open an HTML file containing what you see in your terminal. A big thank you to user progo for this tip.


Upgraded my mac to SSD == pure bliss

Date
image

I recently bought this 240 Gigabytes of awesomeness, conveniently packed into a 2.5" SATA box. It's an Other World Computing Mercury Extreme Solid State Hard Drive. It has no moving parts, it consumes less battery, and it's fast. Damn fast. This thing is so fast it opens applications before you lift your finger from the touchpad after the click. It's also bloody expensive.

I have a MacBook 13" unibody, the first model, also known as MacBook5,1. The hard drive is conveniently seriviceable, right under the large battery lid.

To perform the migration, I needed:

  1. The SSD drive, of course
  2. A smallPhillips screwdriver
  3. A small Torx screwdriver (I managed without one)
  4. A USB-to-SATA external drive case.
  5. Carbon Copy Cloner. The software is free but consider making a donation.

I had two or three options to perform the migration. I don't have the original OSX install disk here with me (and I don't want to go through a total reinstall in any case), but I did have the OSX from my iMac, which may be useful as a safety installation if I screw up. As I quickly learned, the iMac install DVD does not work on the MacBook, but it allows, at least, to perform some basic operation, such as disk partitioning.

My first attempt was based on the assumption that I had to physically copy the partition, low level. I thought this because I wanted to transfer the system while not active (you may have troubles), and I didn't know how the boot worked, so I assumed I had to copy the partition boot record as well... something I did with Linux once. In addition, moving stuff via USB takes a long time, and trasferring 160 GB would have left me with no laptop for a while. So here was the idea:

  1. Disassemble the original disk, and put it in the external case (yes, Mac can boot from this disk once in USB, just press "C" at boot with the USB drive powered and no DVD disk)
  2. Put the OWC disk in the laptop. Partition in two: one with the original size of my old disk, the other with the remaining space
  3. Install a bare OSX on the remaining space, and start copying low-level from the external drive to the other partition. I thought to do this via "Disk Utility" restore.
  4. When the system is done copying, boot from the first partition, delete the second one (with the bare OSX) and find a way to enlarge the partition including the now-free space.

This turned out to be not possible, because the OSX DVD I have does not allow me to install. I could, however, still use this process, using Disk Utility from the DVD, which runs. I tried with no luck, as Disk Utility complains about not being able to allocate memory. In addition, I could no longer boot from USB after this event. The reason was not clear. To address this event, I restored the original disk into the laptop, and proceeded with a different strategy, which is the one I did.

The successful strategy

1. Preparing the new disk

I booted my system back with the old disk, and put the OWC SSD in the USB enclosure. It was clear I had to transfer the system while live, so I logged in with a new bare account with admin privileges, instead of my own. This allows not to have active files or mounted FileVaults, which could give trouble.

I plugged in the SSD in the enclosure, created and formatted a single, whole partition with Disk Utility, exactly in the same format as my system partition. Check in either Disk Utility or from the command line, with the command "diskutil info /dev/disk0s2". It should be Journaled HFS+, because it's the one that allows boot. DO NOT enable case-sensitive. You will have a lot of troubles, in particular with File Vault, (yes, workarounds exist, but I haven't tried them). I loathe case-preserving/case-insensitive FileSystems, but Mac works this way. Also check to use Partition Type "GUID_partition_scheme" when partitioning.

2. Copying the system

Once formatted, I started Carbon Copy Cloner, specified the source as my system disk, the target as the SSD disk, and ticked "backup everything". The documentation in Carbon Copy Cloner is perfect. Check also the "typical usage cases" shown there, you will find the one you need. Make sure Carbon Copy Cloner tells you the cloned system partition will be bootable.

3. Swapping the disks

Carbon Copy Cloner will copy your whole system to the external disk, something that may take many hours. Let it go. When done, I turned off the computer, flipped it, and opened the battery lid. There's a small Phillips screw keeping the hard drive in place via a simple plastic brace. I unscrewed this brace and took it away, then pivoted the hard drive to lift it and unplugged the SATA cable (gently).

Now, I disassembled the SSD drive from the external casing. On the sides of the old hard drive there are four small Torx screws that provide support. The plastic brace kept two of these screws pushed down, while the other two allowed for the pivoting on the side of the battery compartment. I had to remove these screws to put them in the SSD drive, but since I didn't have a Torx screwdriver, I resorted to pliers both for removing them and putting them on. I then plugged in the SATA cable in the SSD, put the disk in the compartment, screwed the brace back and tried to run the laptop with the new setup. Note: conventional wisdom says "first test, then close", but in this case the assembly/disassembly was so trivial I preferred to put everything back in. Also, the battery does not stay in its place without the lid, so you have to close it anyway, or put the battery aside for a while.

4. Booting the new disk

The boot was smooth and uneventful. Not really as fast as I expected, but notably faster. The real difference is in launching applications and logging in. It used to take 10-15 seconds to get to work. Now I have the desktop open in just 2 seconds.

5. Conclusion

I now have the old hard drive and an external USB-SATA case. I could format it and use it as an external disk, or keep it as a safety backup system. I decided for the latter.


base32 encoding in javascript

Date

I had to perform base32 encoding in javascript, and I found nothing ready for the task, so I started cramming out code, and here it is

var baseenc = baseenc || {};

baseenc.b32encode = function(s) {
    /* encodes a string s to base32 and returns the encoded string */
    var alphabet = "ABCDEFGHIJKLMNOPQRSTUVWXYZ234567";

    var parts = [];
    var quanta= Math.floor((s.length / 5));
    var leftover = s.length % 5;

    if (leftover != 0) {
       for (var i = 0; i < (5-leftover); i++) { s += '\x00'; }
       quanta += 1;
    }

    for (i = 0; i < quanta; i++) {
       parts.push(alphabet.charAt(s.charCodeAt(i*5) >> 3));
       parts.push(alphabet.charAt( ((s.charCodeAt(i*5) & 0x07) << 2)
           | (s.charCodeAt(i*5+1) >> 6)));
       parts.push(alphabet.charAt( ((s.charCodeAt(i*5+1) & 0x3F) >> 1) ));
       parts.push(alphabet.charAt( ((s.charCodeAt(i*5+1) & 0x01) << 4)
           | (s.charCodeAt(i*5+2) >> 4)));
       parts.push(alphabet.charAt( ((s.charCodeAt(i*5+2) & 0x0F) << 1)
           | (s.charCodeAt(i*5+3) >> 7)));
       parts.push(alphabet.charAt( ((s.charCodeAt(i*5+3) & 0x7F) >> 2)));
       parts.push(alphabet.charAt( ((s.charCodeAt(i*5+3) & 0x03) << 3)
           | (s.charCodeAt(i*5+4) >> 5)));
       parts.push(alphabet.charAt( ((s.charCodeAt(i*5+4) & 0x1F) )));
    }

    var replace = 0;
    if (leftover == 1) replace = 6;
    else if (leftover == 2) replace = 4;
    else if (leftover == 3) replace = 3;
    else if (leftover == 4) replace = 1;

    for (i = 0; i < replace; i++) parts.pop();
    for (i = 0; i < replace; i++) parts.push("=");

    return parts.join("");
}

I release this code under the obnoxiously named but quite appropriate WTFPL license. I left the decode, as it is commonly said, as an exercise for the reader. Seriously, I didn't need decode, so I did not invest time into it.


Exploring Mandelbrot parameter space - part 2

Date

In the previous post, we saw that graphing the Mandelbrot starting point contains fractal features as well. We want to plot these features, and at the same time increase resolution but reducing computational cost.

The first thing to note is that apparently, the plot is symmetric (but I don't have any strict proof of it) so technically only 1/4th of it can be computed, at least during the exploratory phase. The second optimization is that we are going to plot either a white pixel (if the corresponding Mandelbrot contains at least one white pixel) or a black pixel (the Mandelbrot contains no white pixels). This dramatically reduces the time needed, because when plotting each Mandelbrot, we can exit the generator as soon as a white pixel is found. White pixels are computationally expensive because the iterative procedure for that pixel must run up to the end, hence very whitey Mandelbrots, which are the majority, require a lot of these "up-to-end" iterations. This wastes a lot of time to compute something we know as soon as the first white pixel is found. On the other hand, Mandelbrots containing no white pixels have their iteration loop cut short for basically every pixel, making them faster to compute even if we have to run on each and every pixel of the Mandelbrot.

After these improvements, I ran two calculations. One using 100x100 Mandelbrots (on the left), another using smaller, 10x10 Mandelbrots (on the right)

image

Apparently, the best visual results of the fractalization are obtained with smaller Mandelbrots. I am not really surprised, we already saw a similar situation in the previous post, where I plotted according to the number of white points. I have the feeling that increasing the resolution (the plots you see above are from 1000x1000 images) we would see better fractalization in the 100x100 case as well, but I would not bet on it. In any case, the 10x10 Mandelbrots are sufficient for our purpose, and they are much faster to compute. It's time to increase the resolution up to a staggering 900 Megapixel (30.000 x 30.000), on the full set, and see what happens. What you see took 10 full days to compute. I am not posting the full image, just a small low-resolution screen grab. The image is so big that is unreadable on Firefox, hangs Gimp and puts Preview.app in a very unpleasant situation. Here is the overview

image

As said previously, the Set appears to be symmetric. The fractal features we observe on the border are unique, in particular due to the existence of compartments with common behavior. Please note that all the terms I will use along this analysis are invented by me. Nothing of what you read is formal or agreed nomenclature. I will start from the left hand side of the border, and walk up to the top.

First, we have what I call a spike, a long, narrow line running horizontally and stretching away from the central bulge. This feature is found in the Mandelbrot set as well.

image

If you click the image and zoom in, you will notice a very interesting feature which also can be found in the Mandelbrot Set: lateral whiskers. These curved features intercept the spike in very bright areas. In traditional Mandelbrot, these areas are rich in self-similarity, as you can appreciate from this magnificent Mandelbrot image from Wikipedia. I don't know if similar self-similarity can be found in this Set as well, and the current resolution does not allow to infer any statement. According to what I understand from this answer Sam Nead gave at Math StackExchange, self-similarity is connected to what I assume is a mathematical property called renormalization. Due to my lack of expertise, I cannot explore this any further, but I encourage anyone competent on fractals to add comments, fix my mistakes and provide enlightenment on this regard.

image

The next interesting feature is what I call the prickly pear zone. This zone has very sharp convergence points, and a large number of buds with similar shapes. The central zones look slightly polygonal in nature. I call it prickly pears due to its resemblance with the Opuntia (indian fig) plant, a cactus with characteristic fruits.

image

The sharp contact points recall similar features found in the Mandelbrot set, but on the other hand, the polygonal-like structure is not found. Mandelbrot, at least in its full, unzoomed representation, contains circular structures and the main cardioid.

The third interesting feature is the tree foliage.

image

which appears as a relevant change from the repetitive, sharp-edged prickly pears area. On this regard, however, nothing beats the nervous, highly featured, strongly self-similar and spiraling behavior of the fuzzy worm. This area of the plot is incredibly featured and a real pleasure to zoom in.

image

After the fuzzy worm, we don't see any more spectacular features. The border becomes first rocky, then smooth then rocky again, and it's finally closed with a hole.

image

As usual, I cannot really say anything about the mathematical content of the task I performed and the features of the resulting fractal. However, it has been an incredibly interesting and fun exploration into a fascinating discipline. I hope you enjoyed this exploration, and feel free to leave comments if you want to add information or fix mistakes.


Exploring Mandelbrot parameter space – part 1

Some time ago, I presented an interesting python code able to draw the so-called Mandelbrot set, a fractal image with intriguing properties. Recently, Benoit Mandelbrot passed away. I want to pay homage to his work by digging into more details of his eponymous image.

In the previous post, we observed that a parameter is crucial for the resulting image: the "starting point", described in my code as (a, b) or, in the second part of the post, z0. In traditional Mandelbrot, this starting point has coordinates z0 = (0.0, 0.0), and with this input, the Mandelbrot set is obtained.

After this first evaluation, I went further on by showing how the Mandelbrot Set changes when this point z0 is modified from the original value (0.0, 0.0). The animations I provided show what happens to the Set when choosing different z0 points, such as (1.0, 0.0) or (-2.0, 0.0), for example. The first animation showed plots with z0 = (a, 0.0), with a going from -3.0 to +3.0 with a step 0.1. In other words, I produced a Mandelbrot image for z0 = (-3.0, 0.0), another one for z0 = (-2.9, 0.0) and so on, until z0 = (3.0, 0.0). In the second animated picture, I did the same for points like z0 = (0.0, a) with a on the same interval and step. Finally, the last animation shows how the Mandelbrot behaves with starting points (a, a), again with the same interval.

This exploration of the bond between a possible value for z0 and the resulting Mandelbrot has interesting properties. Let's try to see visually what's going on. Imagine there's a plane of possible values of z0

image

For each point on this plane, there's an associated Mandelbrot Set image. The traditional Mandelbrot is associated to the origin of this plane. With the experiments I did by changing the starting point, I moved along the x axis in the first animated picture, along the y axis in the second animated picture. For presentation purposes, I could choose to stick the actual Mandelbrot images on this plane, imagining that the lower left corner of each image indicates the z0point that image comes from

image

As I already stated, the lower left corner "points to" the z0 point. You see, along the x axis, images for the cases (-2, 0), (-1, 0), (0,0), (1,0), (2,0). Along the y axis, you see (0, -2), (0, -1), (0,0), (0,1), (0,2). These are only a few of the points you could actually choose. Any point of the plane has an associated Mandelbrot.

My plan for this post, is to explore the plane, that is, to see how the Mandelbrot image changes as I select points like (0.34, 0.12), among others. In order to achieve this, I need to do exactly what I did above, but with smaller images with respect to the size of the "graph paper" image. The results are quite interesting

image

The image above shows the z0 plane. It is built with the same concept given for the very simple case above: many Mandelbrot images, whose corner refers to the actual point on the plane. To demonstrate you this is the case, if I zoom in what I see is this

image

A collection of very tiny, very low resolution Mandelbrots, each of them associated to a specific point of the z0 plane. I had to keep the resolution of each individual Mandelbrot to very, very small. What I did takes a lot of computer time.

The interesting fact about this plot is that it contains fractal features as well! Looking carefully on the left, and playing with the zoom, I obtained this image

image

I suggest you observe it from afar. The ramifications typical of fractals will be evident. Every small white dot in the above image is a "mini Mandelbrot" image, whose white points allow us to see something, but not very clearly. Can we improve the situation? Let's see.

I developed a new program. Instead of placing tiny Mandelbrot images like post-it notes on a board, for each point of the z0plane I generate the corresponding Mandelbrot, then I count the number of white points it has, and I color the pixel of the z0plane of a different shade, depending on the amount of "whiteness" showed by the corresponding Mandelbrot. The result ?

image

The above image was created from the number of totally white pixels (scaled from 0 to 255 against the maximum) in a series of 100x100 Mandelbrots, one for each pixel. Not as fancy as I hoped for. I tried to play with color balance to no avail. I think there are two problems. First: the small resolution (1000x1000 is not much). Second, the poor averaging strategy, leading to poor contrast. This plot required 48 hours of computation.

Trying to get better plots requires a lot of computational time to increase the resolution, but unfortunately I only have a laptop. This kind of problem can run parallel very efficiently, but even if I put my second core to work, I won't get very far. The next step is therefore to get a better plot at a reduced cost, eventually accepting some compromises. Stay tuned.


Fortran 90 pitfall: initialization of vars at declaration

Date

I am dusting my Fortran 90 skills. One big gotcha that always leaves me baffled is the following. Suppose you write the following program

program test
  implicit none

  call testsub()
  call testsub()
end program

subroutine testsub()
  implicit none
  integer :: var

  var = 0
  print *, var
  var = 5
  print *, var
end subroutine

If you expect the output to be

0
5
0
5

you are right. This is indeed the output you get.

Now consider the following slight different testsub routine

subroutine testsub()
  implicit none
  integer :: var = 0

  print *, var
  var = 5
  print *, var
end subroutine

See the difference ? I just coalesced the first assignment of var to the declaration line. You wouldn't expect a big difference right ? Sorry to bring the news, but that's a completely different story. The output you will obtain is

0
5
5
5

What happened? The fact arises from a very subtle point of the Fortran standard: local vars with initialization at declaration are automatically SAVE, so they preserve their content between subsequent calls (in C terms, they are local static). In other words, this statement

integer :: var = 0

is totally equivalent to

integer, save :: var = 0

This is totally counter-intuitive and a huge pitfall if you don't know it. What's the rationale behind this practice ? I asked on StackOverflow, and user kemiisto referred me to this page where this behavior is explained in rather detail. The reason is due to historical bad practices.

Apparently, initialization during declaration was already possible in Fortran 77. Usage of this variable without redefinition was allowed behavior, commonly done when you initialize and assign a parameter, for example. On the other hand, redefinition within the routine body was a disallowed practice because, according to the standard, the variable technically became undefined upon reentry.

Before standardization of Fortran 90, the actual internal handling of initialization during declaration was performed with a lot of freedom by compilers, as there was nothing specifying for that in the standard. There were two possible strategies to perform it: static initialization (doing the assignment only once), and reinitialization at every new subroutine call (doing it every time). These two solutions are equivalent, if no reassignment is done inside the subroutine, i.e. the expected practice according to the standard. Compilers were free to choose which strategy to use, but in practice, most compilers used the "initialize once and consider it static" strategy, probably because it's more efficient (you assign only once), so even if the variable was technically undefined if reassignment occurred, in practice it behaved like a static variable.

While the Fortran 90 standard was defined, a lot of code was produced abusing this behavior. As time passed, forbidding it in the new release was not feasible, because it would have introduced a lot of trouble with existing code. This was probably one of those moment in history where programmers would have learned that when something is declared undefined in the standard, it's your fault if you abuse it, and you eventually pay the consequences. The Fortran committee instead condoned and ratified this practice, and now it is part of the standard. Regardless of its status, please follow my advice and stay away from it.


The Mandelbrot set, in python

Date

This code is so fascinating

from PIL import Image

max_iteration = 1000
x_center = -1.0
y_center =  0.0
size = 300

im = Image.new("RGB", (size,size))
for i in xrange(size):
    for j in xrange(size):
        x,y = ( x_center + 4.0*float(i-size/2)/size,
                  y_center + 4.0*float(j-size/2)/size
                )

        a,b = (0.0, 0.0)
        iteration = 0

        while (a**2 + b**2 <= 4.0 and iteration < max_iteration):
            a,b = a**2 - b**2 + x, 2*a*b + y
            iteration += 1
        if iteration == max_iteration:
            color_value = 255
        else:
            color_value = iteration*10 % 255
        im.putpixel( (i,j), (color_value, color_value, color_value))

im.save("mandelbrot.png", "PNG")

If you execute the python code above, the result is this beautiful and mysterious picture:

image

What you see is a Fractal, a highly fractured geometric picture, and more specifically one of the most famous representatives: the Mandelbrot Set. I am not familiar at all with the mathematics behind fractals, so I won't talk about them, or I will simply end up saying something wrong.

Such a highly featured and attractive picture is obtained from a really simple algorithm. It can be explained easily: take a grid of points on the plane, like for example a simple sheet of graph paper. Choose a scale for the axes. For every point (x,y) of the plane now do the following:

  1. start with an empty vector (a,b) containing (0.0, 0.0)
  2. create a new (a,b) vector as (a2 - b2 + x, 2ab + y)
  3. continue with the procedure given in 2. If the length of the vector (a,b) becomes larger than 2.0, plot the (x,y) point black or a shade of gray, depending how far you got in the iterations. Instead, if the length of the vector stays below 2.0 after a given maximum amount of iterations, plot the (x,y) point white.

This, as said, decides the color of one single point (hence, a pixel of the image). Repeat for all the pixels and what you obtain is the Mandelbrot set.

You can also think in terms of complex numbers, plotting on the complex plane. Instead of (a,b) we now have a single complex number z*0*, and instead of (x,y) we now have a complex number c. The repeated operation now can be described as

  1. z*0* = 0.0
  2. c = current complex number expressed by the current pixel point on the complex plane
  3. repeatedly compute*z**n*+1 = z*n*2+ c
  4. if the absolute value of*z* goes beyond 2.0, plot black or gray, otherwise plot white, as given before.

What if you change the starting point? What if we choose z*0* = 1.0 instead of 0.0, or z*0* = 1.0 i ? I tinkered with the program a bit. Here are the results. This image shows how the Mandelbrot Set changes when z*0* goes from -3 to 3 (step 0.1). The set "blows away" towards the left, when going either towards -3 or +3.

mandelbrot-change-x

If you change z*0* along the imaginary axis, going from -3.0 i to 3.0 i, what you obtain is instead:

mandelbrot-change-imaginary

It blows away towards the right. So I said, what happens if I increase both the real part and the imaginary part together? will they "balance out" ? Not really. It burns away

mandelbrot-xy

I have absolutely no idea of what all this means from the strict mathematical point of view, except the fact that I am exploring the space of the values of z*0* and how the Mandelbrot changes accordingly. I thought it was just cool to try it out.