## Equalization

Looking back at my early strange attractors, it is obvious that the images are all very dark. The lightening technique I attempted to use on these images was level normalisation… I stretched the range of values out so that the lightest pixel was pure white, and the darkest pixel is black.

While this is often a good technique for correcting the contrast and exposure in a photograph, it was clearly not doing the job for my strange attractor images. This was because the value histogram is extremely skewed to the left. Ie; most of the pixels are in the very (very) dark range, with only a few very bright pixels. In fact when viewed on-screen, most of these dark pixels are indistinguishable from black. Fortunately these are 16 bit images, so the hidden details are there to be discovered.

That’s where Histogram equalization comes in useful.

This technique actually redistributes the value densities along the histogram so that each value range is more evenly represented. Peaks in the histogram are widened and the troughs condensed. Effectively the surplus of dark shades in my images can be spread up towards the right of the value histogram.

As you can see in these examples of the same image first normalised, and then equalized, equalization really brings these images to life. Note that the grainy texture is due to the number of samples taken and not from the equalization process.

As I wanted to be able to quickly assess a folder full of thumbnails for interesting imaged, it made sense to find a way to programmatically equalize the images in Python before they are first saved to disk.

Luckily I found a suitable algorithm on Salem’s Vision Blog. I just had to make some small alterations so it would work with 16 bit images. I also set the first bin of the histogram to zero – this is to negate the large amount of black background which would otherwise ruin the equalization effect.ote

Here is my version of this code:

def histeq(im,nbr_bins=2**16): '''histogram equalization''' #get image histogram imhist,bins = np.histogram(im.flatten(),nbr_bins,normed=True) imhist[0] = 0 cdf = imhist.cumsum() #cumulative distribution function cdf ** .5 cdf = (2**16-1) * cdf / cdf[-1] #normalize #cdf = cdf / (2**16.) #normalize #use linear interpolation of cdf to find new pixel values im2 = np.interp(im.flatten(),bins[:-1],cdf) return np.array(im2, int).reshape(im.shape)

If you’d rather used an image manipulation program to do this, the effect is available in Photoshop. I believe it is found under Levels >> Adjustments >> Equalize. If you select an area of the image first the equalization for the entire image can be based on the histogram of just the selected portion. This can be very useful for negating a black background, or a white sky etc.

In Linux you can do the same thing with ImageJ, including selecting a region to base the histogram on.

GIMP and CinePaint also have an equalize feature, but no option that I could find to base the histogram on a selected portion of the image. Also GIMP is currently geared towards 8-bit images, so no good in this case.

And for the sake of completeness, here is the latest version of my strange attractor finder. I have tidied up the code and made several optimizations for speed:

#! /usr/bin/python import png import random from math import floor import numpy as np import operator def histeq(im,nbr_bins=2**16): '''histogram equalization''' #get image histogram imhist,bins = np.histogram(im.flatten(),nbr_bins,normed=True) imhist[0] = 0 cdf = imhist.cumsum() #cumulative distribution function cdf ** .5 cdf = (2**16-1) * cdf / cdf[-1] #normalize #cdf = cdf / (2**16.) #normalize #use linear interpolation of cdf to find new pixel values im2 = np.interp(im.flatten(),bins[:-1],cdf) return np.array(im2, int).reshape(im.shape), cdf def test_render(x_coefficients, y_coefficients, size, render_max, trap_size): # starting xy_pos x, y = random.random(), random.random() # initialize array img_array = np.zeros([size, size], int) a, b, c, d, e, f = x_coefficients g, h, i, j, k, l = y_coefficients while True: # main rendering loop x, y = (a*(x**2) + b*x*y + c*(y**2) + d*x + e*y + f, g*(x**2) + h*x*y + i*(y**2) + j*x + k*y + l) # translate x, y into array indices px = int((x+trap_size)*(size/(trap_size*2))) py = int((y+trap_size)*(size/(trap_size*2))) if not(px and py): return (None, None) # catch negative numbers... the "IndexError" won't # do this resulting in "wrap-around" images try: img_array[px, py] += 1 except IndexError: return (None, None) if img_array[px, py] == render_max: return img_array, (x, y) def final_render(x_coefficients, y_coefficients, size, render_max, trap_size, xy_pos): # starting xpos x, y = xy_pos # initialize array img_array = np.zeros([size, size], int) a, b, c, d, e, f = x_coefficients g, h, i, j, k, l = y_coefficients while True: # optimised rendering loop with no error checking etc. x, y = (a*(x**2) + b*x*y + c*(y**2) + d*x + e*y + f, g*(x**2) + h*x*y + i*(y**2) + j*x + k*y + l) # translate x, y into array indices px = int((x+trap_size)*(size/(trap_size*2))) py = int((y+trap_size)*(size/(trap_size*2))) # don't test for negatives # don't catch exceptions... img_array[px, py] += 1 if img_array[px, py] == render_max: return img_array if __name__ == "__main__": import psyco psyco.full() trap_range = 2 # potential attractor is rejected if xy_pos # falls further than this distance from origin painted_area_threshold = .01 # reject attractor if less than this portion # the whole has been painted test_render_size = 200 test_render_max = 200 final_render_size = 800 final_render_max = 100000 equalize_final = True img_index = 0 imgWriter = png.Writer(final_render_size, final_render_size, greyscale=True, alpha=False, bitdepth=16) while True: # Set the tweakable values of the attractor equations: random.seed() seed = int(random.random() * 10000000) random.seed(seed) x_coefficients = [random.choice([random.random()*4-2, 0, 1])\ for x in range(6)] y_coefficients = [random.choice([random.random()*4-2, 0, 1])\ for x in range(6)] # perform small test render img_array, xy_pos = test_render(x_coefficients, y_coefficients, test_render_size, test_render_max, trap_range) #only proceed if test attractor did not fail if not xy_pos: continue #only proceed if painted area meets threshold, otherwise loop to top painted_area = sum(img_array[img_array==True].flatten()) / test_render_size ** 2. if painted_area < painted_area_threshold: continue # perform final render print "rendering: %06d_%d.png" % (img_index, seed) try: img_array = final_render(x_coefficients, y_coefficients, final_render_size, final_render_max, trap_range, xy_pos) except IndexError: continue # histogram equalization if equalize_final: img_array, cdr = histeq(img_array) # save to final render to png file f = open("a_%06d_%d.png" % (img_index, seed), "wb") imgWriter.write(f, img_array) f.close() img_index += 1 print "SAVED!"