H. Chase Stevens Logo
blog · schedule · résumé · public key · contact

Most Common Tags

  1. programming
  2. python
  3. code
  4. philosophy
  5. evolution
  6. game design
  7. probability
  8. video games
  9. genetic algorithms
  10. government



Posts with tag "code":

Generating Contrast and Saliency Maps on the GPU (Using Python and OpenCL)

posted on: Saturday, April 13, 2013 (8:21 am) by Chase Stevens

In the field of Computational Neuroscience, saliency maps are a means of graphically representing the areas of any visual scene presenting the most "bottom-up" saliency to a human observer (i.e. those most likely to draw the viewer's attention). Although the generation of these maps is not particularly difficult on a conceptual level, doing so is quite computationally expensive if using a serial approach. Below, I provide code for quickly generating the component contrast maps needed to build a saliency map by parallelizing the task on the GPU, as adapted from MATLAB code provided by Vicente Ordonez of SUNY. To run this, you'll need pyopencl v0.92, numpy, and PIL.

Christian Bale Saliency map of Christian Bale
Original image Saliency map

from PIL import Image
import itertools
import numpy as np
import pyopencl as cl
import pyopencl.array as cl_array
from pyopencl.elementwise import ElementwiseKernel
import math

## initialization of GPU context, queue, kernel
ctx = cl.create_some_context()
queue = cl.CommandQueue(ctx)
kernel_args = "int width, int len_min_width, float *image, float *contrast"
kernel_code = '''
contrast[i/3] = contrast[i/3] +
     (i < width) ||
     (i > len_min_width) || 
     (i % width == 0) ||
     ((i+1) % width == 0) ||
     ((i+2) % width == 0) ||
     ((i+3) % width == 0) ||
     ((i-1) % width == 0) ||
     ((i-2) % width == 0)
     pown((image[i] - image[i-width]), 2) +
     pown((image[i] - image[i-width-3]), 2) +
     pown((image[i] - image[i-width+3]), 2) +
     pown((image[i] - image[i-3]), 2) +
     pown((image[i] - image[i+3]), 2) +
     pown((image[i] - image[i+width]), 2) +
     pown((image[i] - image[i+width-3]), 2) +
     pown((image[i] - image[i+width+3]), 2)
 ) : 0''' # for each rgb value in image, if not part of a pixel on the image
          # border, sum the squares of the differences between it and the
          # corresponding color values of each surrounding pixel, then add
          # this to the corresponding output pixel's value
contrast_kernel = ElementwiseKernel(ctx,

flatten = lambda l: list(itertools.chain(*l))

def contrast_map(image):
    width = image.size[0]
    height = image.size[1]
    ## creating numpy arrays
    image = np.array(flatten(list(image.getdata()))).astype(np.float32) #image array
    contrast_map = np.zeros((height*width)).astype(np.float32) #blank write array
    ## send arrays to the gpu:
    image_gpu = cl_array.to_device(ctx,queue,image) 
    contrast_map_gpu = cl_array.to_device(ctx,queue,contrast_map) 
    contrast_kernel(width*3,(image.size - width - 1),image_gpu,contrast_map_gpu) #executing kernel
    contrast_map = contrast_map_gpu.get().astype(np.float32) #retrieving contrast map from gpu
    contrast_map = np.nan_to_num(contrast_map) #conversion of NaN values to zero
    ## normalization:
    contrast_map += max(contrast_map.min(),0)
    contrast_map /= contrast_map.max()
    contrast_map *= 255
    return contrast_map.astype(np.uint8)

def saliency_map(image):
    width = image.size[0]
    height = image.size[1]
    resizes = int(math.floor(math.log(min(width,height),2)) - 3) #calculating necessary number of
                                                                 #images for gaussian pyramid
    resized_images = [image.resize((width/factor,height/factor),Image.BICUBIC) 
                      for factor in [2**x for x in range(resizes)]] #resizing images
    contrast_resized_images = map(contrast_map, resized_images) #generating contrast maps

    ## resizing contrast maps back to original image size:
    images = list()
    for i in range(len(contrast_resized_images)):
        image = Image.new("L",resized_images[i].size)

    ## combining images:    
    result = np.zeros_like(images[0]).astype(np.float32)
    for image in images:
        result += image.astype(np.float32)
    ## normalization:
    result += max(result.min(),0)
    result /= result.max()
    result *= 255
    output = Image.new("L",(width,height))
    return output

Tags: code, cognitive science, gpgpu, image processing, numpy, opencl, programming, python

Visualizing k-means clustering (using d3.js)

posted on: Friday, March 22, 2013 (8:22 am) by Chase Stevens

Having learned a little about d3.js at a hackathon I attended, I decided to put it to good use. Below is an interactive implementation of the k-means clustering algorithm, as visualized using d3.js. Click to place points. Centroids are (initially) set on randomly chosen points, and the number of clusters is set by a supposed rule of thumb. You can choose to step through the algorithm one centroid-readjustment/simultaneous point assignment at a time, or view the process in real time.


//config values:
var colors = ['coral','green','steelblue','seagreen','yellow','cyan','salmon','skyblue','pink']
var radius = 8;
var trans_secs = 2;
var update_ms = 150;
var svg_size = {width:600,height:400};

var mouse_x;
var mouse_y;
var new_points = false;
var points = new Array();
var centers = new Array();
var algo_started = false;
var svg = d3.select('#kmeans').append('svg')
canvas = document.getElementById('kmeans');
var stream = false;
var stop_update = false;

//@ http://jsfromhell.com/array/shuffle [v1.0]
function shuffle(o){ //v1.0
    for(var j, x, i = o.length; i; j = parseInt(Math.random() * i), x = o[--i], o[i] = o[j], o[j] = x);
    return o;

function safeShuffle(o){
	var x = new Array();
	for(var i = 0; i < o.length; i++){ 
		x[i] = o[i]; 
	return shuffle(x);

function addPoint(){
	//add new point to canvas (d3)
	if (algo_started) {return true;}
	var rel_x = mouse_x - canvas.offsetLeft;
	var rel_y = mouse_y - canvas.offsetTop;
	if (points.length >= 1 && (points[points.length-1].x == rel_x)
		&& (points[points.length-1].y == rel_y)){ //making sure double-clicks don't get counted as two points
		return true;
	new_points = true;

//tracking mouse:
document.onmousemove = getMousePosition;
function getMousePosition(mp) {
	mouse_x = mp.pageX;
	mouse_y = mp.pageY;
	return true;

function capitalize(string)
    return string.charAt(0).toUpperCase() + string.slice(1);

function createCenters () {
	//creating centroids
	document.getElementById('create').style.display = 'none';
	document.getElementById('move').style.display = 'block';
	algo_started = true;
	rand_points = safeShuffle(points);
	for (var i = 0; i < Math.round(Math.sqrt(rand_points.length / 2)); i++){
		centers.push({x:rand_points[i].x, y:rand_points[i].y, color:colors[i%colors.length], old_x:null, old_y:null});
	.attr('cx', function(d){return d.x;})
	.attr('cy', function(d){return d.y;})
	.attr('r', radius)
	.attr('style', function(d){return 'fill:' + d.color;})
   	.text(function(d){return capitalize(d.color) + ' cluster center'});

function moveCenters () {
	//recalculate centroid positions given their assigned points
	for (var i = 0; i < centers.length; i++){
		var center = centers[i];
		center.old_x = center.x;
		center.old_y = center.y;
		var new_x = 0;
		var new_y = 0;
		var point_count = 0;
		for (var j = 0; j < points.length; j++){
			if (points[j].nearest_center == center){
				new_x += points[j].x;
				new_y += points[j].y;
		center.x = Math.round(new_x / point_count);
		center.y = Math.round(new_y / point_count);
	.duration(trans_secs * 1000)
	.attr('cx', function(d){return d.x;})
	.attr('cy', function(d){return d.y;});

	//check to see if centroids have failed to move any more (i.e. algorithm complete)
	var all_finalized = true;
	var svg_centers = svg.selectAll('.center')[0];
	for (var i = 0; i < centers.length; i++){
		var center = centers[i];
		var svg_center = {x: Math.round(svg_centers[i].cx.baseVal.value),
			y: Math.round(svg_centers[i].cy.baseVal.value)};
		all_finalized &= (center.x == center.old_x 
			&& center.y == center.old_y 
			&& Math.round(center.x) == svg_center.x 
			&& Math.round(center.y) == svg_center.y
	if (all_finalized){
		stream = false;
		document.getElementById('movecenters').disabled = true;
		document.getElementById('togglestream').disabled = true;
		document.getElementById('final').innerHTML = "Finalized!";
		stop_update = true;

function eucDist(a,b){
	return Math.sqrt(Math.pow((a.x - b.x),2) + Math.pow((a.y - b.y),2));

function updatePoints(){
	//either creates new svg objects based on clicked points,
	//or re-colors points based on current centroid locations
	if (new_points && !algo_started){
		.attr('cx', function(d){return d.x;})
		.attr('cy', function(d){return d.y;})
		.attr('r', radius)
		new_points = false;
	if (algo_started){
		for (var i = [0]; i < points.length; i++){
			var point = points[i];
			var center_candidate = null;
			var center_candidate_dist = 99999;
			var svg_centers = svg.selectAll('.center')[0]
			for (var j = 0; j < centers.length; j++){
				var svg_center = svg_centers[j];
				center_dist = eucDist(point,{x:svg_center.cx.baseVal.value,y:svg_center.cy.baseVal.value});
				if (center_dist < center_candidate_dist){
					center_candidate_dist = center_dist;
					center_candidate = centers[j];
			point.nearest_center = center_candidate;
			point.color = 'light' + point.nearest_center.color;
		.attr('style', function(d){return 'fill:' + d.color;});

function updateDisplay(){
	document.getElementById('numpoints').innerHTML = points.length;
	document.getElementById('numclusters').innerHTML = Math.round(Math.sqrt(points.length / 2));

function toggleStream(){
	stream ^= true;

function update(){
	if (stop_update){return true;}
	if (stream) { moveCenters(); }

Tags: code, d3, javascript, machine learning, programming, web

Reconstructing Images Using Damaged Copies (in Python)

posted on: Thursday, September 13, 2012 (4:01 pm) by Chase Stevens

I've recently been working on implementing various methods of image reconstruction in python. The idea is to, given several imperfect (that is to say, noisy, incomplete, photoshopped, or otherwise damaged) copies of some original image, attempt to arrive at something close to (if not precisely) the original image by combining said copies. Through these means, an approximation of the original image can be generated should the original be lost, itself damaged, irretrievable, or otherwise unavailable or useless. In writing the various functions for doing this, I implemented techniques used in signal averaging and variants thereof. I also implemented a “modal image" function which, for each pixel, uses the “modal" RGB values across all image copies or, failing that, performs a simple mean of values.

Examples and Analysis

Original Christian Bale photograph

For the following examples, I modified the above image of actor Christian Bale. Ironically enough, in testing for this article, I overwrote the original image and had to employ the use of Google's reverse image-search in order to find it.

Damaged images
Damaged Christian Bale #1 Damaged Christian Bale #2 Damaged Christian Bale #3
Delta: 105.122025 Delta: 88.026425 Delta: 68.5942
Function results (listed with difference from original image as given by get_delta, lower is better):
<function average_images_add at 0x00000000030E8F28> : 155.35693875
<function average_images_sub at 0x00000000030E95F8> : 72.4844575
<function average_images at 0x0000000002EF0208> : 43.92254625
<function average_noisefilter_all_sub at 0x0000000002EF0278> : 51.1805645833
<function average_noisefilter_all_delta at 0x0000000002EF02E8> : 36.9071316667
<function modal_image at 0x0000000002EF0358> : 42.53322
Christian Bale Reconstructed Using modal_image Christian Bale Reconstructed Using average_noisefilter_all_delta
Function: modal_image
Delta: 42.53322
Function: average_noisefilter_all_delta
Delta: 36.9071316667
As is readily visible, in this example the naïve “voting”-type approach used by modal_image is deceived in any case where forms of damage across multiple images “agree” with each other. Given a larger number of copies or a less artificial form of damage, this would likely cease to be an issue; in theory, modal_image could even go so far as to reconstruct the original image perfectly. average_noisefilter_all_delta produced the best results on average, although, to its detriment, its output relies on the order of the list of image copies passed to it. In addition, while it manages to be closer to the original image, the subtraction of image differences it employs creates a slightly “jarring” visual effect (as seen above). The inherent issue in reconstructing damaged images is one of knowledge. To humans, it seems obvious that the original image of Christian Bale didn't have streaks of white across it. However, this is a mere assumption based on our vast pool of experience in viewing photographs. Who's to say that the original photo didn't have white scribbles on it? The computer is hampered, from our perspective, by its inability to make these assumptions when reconstructing the original image, so although it produces results better than a mere superposition of the copies, they rarely are better than what could be accomplished by human hands.

Noisy images
Noisy Christian Bale #1 Noisy Christian Bale #2 Noisy Christian Bale #3
Delta: 92.715475 Delta: 93.352175 Delta: 93.08885
Function results (listed with difference from original image as given by get_delta, lower is better):
<function average_images_add at 0x00000000030E8F28> : 103.01672125
<function average_images_sub at 0x00000000030E95F8> : 81.929801875
<function average_images at 0x0000000002EF0208> : 82.971985
<function average_noisefilter_all_sub at 0x0000000002EF0278> : 69.6356495833
<function average_noisefilter_all_delta at 0x0000000002EF02E8> : 75.23682
<function modal_image at 0x0000000002EF0358> : 65.2880416667
Christian Bale De-Noised Using modal_image Christian Bale De-Noised Using average_noisefilter_all_sub
Function: modal_image
Delta: 65.2880416667
Function: average_noisefilter_all_sub
Delta: 69.6356495833
In this case, modal_image produced the best results, although its output is still quite noisy. Again, having more copies would significantly improve the end product. The output also appears to be slightly brighter than would be expected. average_noisefilter_all_sub comes in at a close second, although (to the human eye) its results appear to be quite a bit noisier than modal_image.

I tried to use these image reconstruction methods on compressed images with jpeg artifacts as well, although the results were much poorer. Output deltas ended up being worse than the least compressed copy of a given set, although better than an average copy. modal_image and average_images_add seemed to perform best. The overall problem was again one of knowledge: could less compressed images be prioritized or weighted given their closeness to the source, results would likely be much better. However, the computer has no means by which to determine what is closest to the original, and thus fails to leverage this. Some kind of programmatic detection of jpeg artifacts could be of great help in improving this situation.

As a final note before my code, I just recently launched Broverbs. Please do check it out if you've got the time!

from __future__ import division
from PIL import Image as PILImage
from Tkinter import *
from tkFileDialog import askopenfilename as openfile
from itertools import product
from random import shuffle

class Image():
    def __init__(self,filename=None,dimensions=None):
        if filename!=None:
            self.pic = PILImage.open(filename)
            self.imgData = self.pic.load()
            self.size = self.pic.size
            if dimensions!=None:
                self.size = dimensions
                self.size = [0,0]
            self.pic = None
            self.imgData = {}
            for x,y in product(range(self.size[0]),range(self.size[1])):
    def setpixel(self,x,y,r,g,b):
        self.imgData[x,y] = tuple(map(max,zip((r,g,b),(0,0,0))))
    def getpixellist(self):
        return product(range(self.size[0]),range(self.size[1]))
    def show(self):
        tempImage = PILImage.new('RGB',self.size)
        for x,y in self.getpixellist():

def get_delta(image1,image2):
    if image1.size != image2.size: raise ValueError("Image sizes do not match")
    delta = 0
    for x,y in image1.getpixellist():
        delta += sum(abs_pixel_diff(image1,image2,x,y))
    return delta / (image1.size[0] * image1.size[1])

mode = (lambda l: max([(l.count(i), i) for i in set(l)])[1] if len(set(l)) < len(l) else sum(l)/float(len(l)))

pixel_operation = lambda f: lambda image1,image2,x,y: map(f,zip(image1.imgData[x,y],image2.imgData[x,y]))
abs_pixel_diff = pixel_operation(lambda (x,y): abs(x-y))
pixel_diff = pixel_operation(lambda (x,y): x-y)
pixel_avg = pixel_operation(lambda (x,y): (x+y)/2)
pixel_sum = pixel_operation(lambda (x,y): x+y)

def image_operation(operation,image1,image2):
    if image1.size != image2.size: raise ValueError("Image sizes do not match")
    result = Image(dimensions=image1.size)
    for x,y in result.getpixellist():
        r,g,b = operation(image1,image2,x,y)
    return result

get_delta_image = lambda image1,image2: image_operation(abs_pixel_diff,image1,image2)
subtract_image = lambda minuend,subtrahend: image_operation(pixel_diff,minuend,subtrahend)
average_image = lambda image1,image2: image_operation(pixel_avg,image1,image2)
add_image = lambda image1,image2: image_operation(pixel_sum,image1,image2)

def get_average_image(image_list):
    output = Image(dimensions=set1[0].size)
    for x,y in output.getpixellist():
        r,g,b = reduce(lambda a,b: (a[0]+b[0],a[1]+b[1],a[2]+b[2]),[image.imgData[x,y] for image in image_list])
    return output

def generalized_average(image_list, combination_function, function1, function2):
    set1 = image_list[0::2]
    image1 = get_average_image(set1)
    set2 = image_list[1::2]
    image2 = get_average_image(set2)
    return combination_function(function1(image1,image2),function2(image1,image2))

def average_images_add(image_list): 
    return generalized_average(image_list,add_image,average_image,subtract_image)

def average_images_sub(image_list):
    return generalized_average(image_list,subtract_image,average_image,subtract_image)

def average_images(image_list): 
    return generalized_average(image_list,subtract_image,average_image,get_delta_image)

def generalized_noisefilter(image_list,function):
    set1 = image_list[0::2]
    image1 = get_average_image(set1)
    set2 = image_list[1::2]
    image2 = get_average_image(set2)
    noise = function(image1,image2)
    denoise_set = [subtract_image(image,noise) for image in image_list]
    return get_average_image(denoise_set)

def average_noisefilter_all_sub(image_list):
    return generalized_noisefilter(image_list,subtract_image)

def average_noisefilter_all_delta(image_list):
    return generalized_noisefilter(image_list,get_delta_image)

def modal_image(image_list):
    result = Image(dimensions=image_list[0].size)
    for x,y in result.getpixellist():
        r = mode([image.imgData[x,y][0] for image in image_list])
        g = mode([image.imgData[x,y][1] for image in image_list])
        b = mode([image.imgData[x,y][2] for image in image_list])
    return result

def test_func(f,images,original,trials=25):
    results = list()
    for x in range(trials):
        print "Trial #"+str(x+1)
    return sum(results)/len(results)

function_list = [average_images_add,average_images_sub,average_images,average_noisefilter_all_sub,average_noisefilter_all_delta,modal_image]

def test_functions(image_list,original,functions=function_list,trials=10):
    out = ''
    for f in functions:
        out += (str(f) + ' ')
        out += (': ')
        out += (str(test_func(f,image_list,original,trials)))
        out += ('n')
    print out

Tags: code, image processing, information theory, programming, python, voting

More Stupid Higher-Order Function Tricks

posted on: Thursday, March 15, 2012 (11:23 am) by Chase Stevens

A while ago, a friend of mine who was unfamiliar with python asked me how to go about parsing a CSV into a list of lists. Moreover, the data included integers which needed to be parsed as such (and not as strings). My solution was as follows:

(lambda f: [map((lambda x: x.isdigit() and int(x) or x),(line.split(','))) for line in f.readlines()]) #takes a file
Just yesterday, another friend of mine had accidentally merged and subsequently sorted two excel spreadsheets, and overwrote one with the result. She asked if I could take the merged file and the remaining original file and re-construct the file she had overwritten. Incorporating the CSV parsing function I had written earlier, I wrote this:
from tkFileDialog import asksaveasfile, askopenfile
full = (lambda f: [map((lambda x: x.isdigit() and int(x) or x),(line.split(','))) for line in f.readlines()])(askopenfile())
exclude = (lambda f: [map((lambda x: x.isdigit() and int(x) or x),(line.split(','))) for line in f.readlines()])(askopenfile())
(lambda f: map(f.write,(map(','.join,filter((lambda x: x not in exclude),full)))))(asksaveasfile())

Tags: code, csv, filter, lambda, map, programming, python, tkinter

One-Line Mode Function

posted on: Tuesday, March 13, 2012 (12:59 pm) by Chase Stevens

More functional programming fun.

(lambda l: max([(l.count(i), i) for i in set(l)])[1])

Tags: lambda, math, programming, python, statistics, code