home | lab | find me | science | publications | software | toolbox | site map

Doodle Programming

Index


Printing your name in binary

2008-04-09

#include <stdlib.h>

void printBits(int x) {
   int n;
   for(n=0; n<8; n++) {
      if (0 != (x & 0x80)) { printf("1"); }
      else { printf("0"); }
      x = x<<1;
   }
}

int main() {
    char* name = "Marta";
    int i;
    for (i=0; i<5; i++) {
        printBits((int)name[i]);
        printf(" ");
    }
    return 0;
}

Which outputs the following lovely sequence:

01001101 01100001 01110010 01110100 01100001

I developed the code above dynamically despite being C, thanks to codepad.org.


Interpolating RGB colors in java, mixing hues the right way

2008-03-01

What we want: red and green makes yellow. A simple average of the RGB values would not work, except for all 255 values.

So what we do instead is to convert the RGB color to HSB (Hue, Saturation and Brightness), and then take the average of the hue. The hue axis of the HSB color space is, though, a circle. The average needs to be the closest one - hence the if/else in the code.

/** Mix colors visually: red + green = yellow, etc.*/
static public final Color mix(final Color c1, final Color c2) {
	final float[] b = Color.RGBtoHSB(c1.getRed(), c1.getGreen(), c1.getBlue(), new float[3]);
	final float[] c = Color.RGBtoHSB(c2.getRed(), c2.getGreen(), c2.getBlue(), new float[3]);
	final float[] a = new float[3];
	// Find to which side the hue values are closer, since hue space is a a circle
	// Hue values all between 0 and 1
	float h1 = b[0];
	float h2 = c[0];
	if (h1 < h2) {
		float tmp = h1;
		h1 = h2;
		h2 = tmp;
	}
	float d1 = h2 - h1;
	float d2 = 1 + h1 - h2;
	if (d1 < d2) {
		a[0] = h1 + d1 / 2;
	} else {
		a[0] = h2 + d2 / 2;
		if (a[0] > 1) a[0] -= 1;
	}
	// only Saturation and Brightness can be averaged directly:
	for (int i=1; i<3; i++) a[i] = (b[i] + c[i]) / 2;
	return Color.getHSBColor(a[0], a[1], a[2]);
}

The result is the expected for the yellowish and reddish tones of the two pipes on the left: an orangish color. And for blue and yellow, we get green.

 

The pipes shown above correspond to Drosophila brain secondary lineage bundle tracts, as registered to each other according to the axes provided by their mushroom bodies (see TrakEM2 for how to do this yourself). Blending the colors correctly for the interpolated pipe was a must.

My appreaciation to Stephan Saalfeld, who helped me hack the above, and realise the circularity of the hue axis. We both came up with the same solution, but mine had an error ... oh well.


Creating cubic mountains in Blender

2008-03-11

# Albert Cardona 2008
# Requires 2.45 or above
# Create a grid of cubes on XY, and scale their Z according to several overlapped waves

from Blender import *
from math import sin, pi, pow

editmode = Window.EditMode()
if editmode: Window.EditMode(0)

# a cube to share from
src = Mesh.New('cube')
src.verts.extend([[0,0,0], [1,0,0], [1,1,0], [0,1,0], [0,0,1], [1,0,1], [1,1,1], [0,1,1]])
src.faces.extend([[0, 1, 2, 3], [0, 1, 5, 4], [1, 2, 6, 5], [2, 3, 7, 6], [3, 0, 4, 7], [4, 5, 6, 7]])

# some color, in a single vertex
src.vertexColors = 1
src.faces[4].col[0].r = 255
src.faces[4].col[0].g = 128
src.faces[4].col[0].b = 0

scene = Scene.GetCurrent()
obsrc = scene.objects.new(src, 'src') # must be linked to be able to share from

def getZ(i, j):
	z = 0
	z += pow( sin( i / pi ) + 1, 2 )
	z += ( sin( j / pi ) + 1 )
	z += ( sin( j / pi ) + 1 )
	return z

# make a n * m grid
n = 100
m = 100
for i in range(n):
	for j in range(m):
		cube = Object.New('Mesh')
		#cube.link(src)
		cube.shareFrom(obsrc)
		cube.setLocation(i, j, 0)
		cube.setSize(1, 1, getZ(i, j))
		scene.objects.link(cube)

# remove template
scene.objects.unlink(obsrc)

if editmode: window.EditMode(1)
Redraw()

Creating the spectrogram of a song in Matlab or Octave

We take the last minute or so of an amazing Yngwie Malsteem guitar solo, and create a two-dimensional display of its sound frequencies versus time: a spectrogram.

addpath("/usr/share/octave/site/api-v22/m/octave2.9-forge/signal/")
[song,fs] = wavread("/home/albert/Music/guitar_solo.wav");
% Convert from stereo to mono, by obtaining the average of both channels
song = (song(:,1) + song(:,2))/2;
% start at the 4th minute
start = 44100 * 4 * 60
% crop song
song = song(start:end);
% set proper LUT - 'hot' only in Matlab
% In octave you have to create it yourself as a 3-column (RGB) matrix
% colormap("hot");
% display spectrogram
specgram(song, 2048, fs);

% Alternatively, call like this to avoid displaying
% the image but capturing the axes for finer plotting
% [S, f, t] = specgram(song, 2048, fs);

% one can tweak the colormap (it's a built-in var) to make it colorful

My gratitude to Ross Maddox who showed me how to do this last Summer 2007 in Telluride's neuromorphic engineering workshop (INE).

To install Octave in Ubuntu, you need two packages:

$ sudo apt-get install octave2.9 octave2.9-forge

... which you can install from synaptic as well, of course.


Converting a numerical sequence 1,2,3,...,26,27,... to a character sequence A,B,C,...,Z,AA,AB,... in Java

We can do so with a recursive function, and using the known reltionship between chars and ints:

static public final String getCharacter(int i) {
	i--;
	int k = i / 26;
	char c = (char)((i % 26) + 65); // 65 is 'A'
	if (0 == k) return new StringBuffer().append(c).toString();
	return new StringBuffer().append(getCharacter(k)).append(c).toString();
}

The above is part of the Utils class in TrakEM2 0.5n.

And yes: it could be done a lot more efficiently, computing the length of the char[] array and passing that to a private final function to fill it up, and then returning a new String made of it. If you feel like writing it, I'd appreciate if you send it to me.


Printing the numerical value of characters in Java

Sometimes one needs to see the integer correspondences between char and int values. So here is how you can print it. Just save into a file named "Print_Char_Index.java" and compile and run it:

public class Print_Char_Index {
	static public void main(String[] arg) {
		for (int i=-128; i<300; i++) {
			System.out.println(i + " : " + (char)i);
		}
	}
}

Using the above comes in handy in many situations, such as for knowing where does the uppercase alphabet starts (at 65).


Automatically placing image windows in ImageJ

There are two components: a plugin that listens for image opened events, and a macro to run the plugin when ImageJ is started up.

Edit your macros/StartupMacros.txt file and find the AutoRun macro (if the file and/or AutoRun macro don't exist, just create them.) Then write in the following:

macro "AutoRun" {
    run("ImageWindow Position");
}

Now create a file named ImageWindow_Position.java under the plugins folder, and copy the following into it:

import ij.*;
import ij.plugin.PlugIn;
import java.awt.*;

public class ImageWindow_Position implements PlugIn, ImageListener {
    public void run(String arg) {
        ImagePlus.addImageListener(this);
    }
    public void imageOpened(ImagePlus imp) {
        ImageWindow win = imp.getWindow();
        if (null == win) return;
        Dimension screen = Toolkit.getDefaultToolkit().getScreenSize();
        Rectangle box = win.getBounds();
        win.setLocation(screen.width - box.width,
                     (screen.height - box.height) / 2);
    }
    public void imageUpdated(ImagePlus imp) {}
    public void imageClosed(ImagePlus imp) {}
}

To compile the file, go to the menu "Plugins / Compile and Run".

Then restart ImageJ and your images will be placed in the right side of the screen automatically.


Counting lines of code, comments and total number of characters for a java code repository

In the example, for TrakEM2:

# Navigate T2 folder and compose one big list containing all .java code lines

import os

t2dir = '/home/albert/T2-src/'

def put(path, files, lines):
	files.append(path)
	f = open(path, 'r')
	for line in f:
		lines.append(line)
	f.close()

def read_all(source_dir, list, lines):
	for root, dirs, files in os.walk(source_dir):
		for name in files:
			if name.endswith('.java'):
				put(root + '/' + name, list, lines)
		for dir in dirs:
			read_all(dir, list, lines)

list = []
lines = []
read_all(t2dir, list, lines)

print "T2 lines: ", len(lines), "from ", len(list), " files."

# statistics
count_blank = 0
count_comment = 0
count_head_comment = 0
in_head_comment = False
in_code_comment = False
count_after_code_comment = 0
n_chars = 0

for line in lines:
	n_chars += len(line)
	line = line.strip()
	if 0 == len(line):
		count_blank += 1
	elif not in_code_comment and not in_head_comment and line.startswith('/**'):
		count_head_comment += 1
		if -1 != line.find('*/'):
			in_head_comment = True
	elif in_head_comment and -1 != line.find('*/'):
		count_head_comment += 1
		in_head_comment = False
	elif line.startswith('//'):
		count_comment += 1
	elif not in_code_comment and not in_head_comment and line.startswith('/*'):
		count_comment += 1
		if -1 != line.find('*/'):
			in_code_comment = True
	elif in_code_comment and -1 != line.find('*/'):
		count_comment += 1
		in_code_comment = False
	elif -1 != line.find('//'):
		count_after_code_comment += 1
	# else it's a line of code

print "Lines of code:", len(lines) - count_head_comment - count_comment
print "Lines of doc comments:", count_head_comment
print "Lines of code comments:", count_comment
print "Lines of code with comments appended:", count_after_code_comment
print "Total number of characters:", n_chars

The printed result looks like, as of 2008-03-02:

T2 lines:  64644 from  134  files.
Lines of code: 58109
Lines of doc comments: 1568
Lines of code comments: 4967
Lines of code with comments appended: 1761
Total number of characters: 2239517

Fetch a web page with python and parse its content to do some reverse DNS lookups

The situation: my server admins have decided to shut down DNS lookup for performance reasons. So the .addresses file generated from the file tracking php script (see next entry) contains IP numbers only from now on.

Instead of pulling my hairs in despair, I can just use python to do the trick, spawning a native OS process to do the job: either dig -x <some ip number> or even easier, using the host command.

While I was at it, I added the ability to add an optional command line argument to limit the amount of entries I wanted listed.

Note how all file reading activity is done within a try - except block, to ensure that even in the event of trouble, the file can be closed afterwards (otherwise, your operating system kernel may accumulate unclosed file handles, and eventually run out of them.)

# fix non-resolved IPs from the TrakEM2_.jar.address file

from urllib import urlopen
import traceback, os, sys

# locking version of the process spawning and output capturing
def capture(command):
	lines = []
	p = os.popen(command)
	try:
		while 1:
			line = p.readline()
			if not line: break
			lines.append(line)
	except:
		traceback.print_exc()
	p.close()
	return lines

# default maximum number of lines to output
max = 20
# capture optional command line argument
if len(sys.argv) > 1:
	max = sys.argv[1]
	if not max.isdigit():
		print max, " is not a digit"
		max = 20
	else:
		max = int(max)


url = urlopen('http://www.ini.uzh.ch/~acardona/jars/TrakEM2_.jar.addresses')

try:
	lines = url.readlines()
	lines.reverse()
	count = 1
	for line in lines:
		if count > max:
			break
		count += 1
		line = line.strip()
		# check if address is resolved
		if not line[:line.find('.')].isdigit():
			print line
			continue
		# otherwise, dig it
		b = line.split('	')
		# returns an int, not the string output
		#dig = os.system('dig -x ' + b[0])
		# dig = dig.split('\n')
		addr = None
		"""
		dig = capture('/usr/bin/dig -x ' + b[0])
		for i in range(len(dig)):
			if -1 != dig[i].find('ANSWER SECTION'):
				addr = dig[i+1]
				addr = addr[addr.rfind('	')+1:-2]
		"""
		# Even easier: (less parsing)
		host = capture('/usr/bin/host ' + b[0])[0]
		if -1 == host.rfind('not found:'):
			addr = host[host.rfind(' '):-1]

		if addr is None:
			print b[0], b[2]
		else:
			print addr, b[0], b[2]
except:
	traceback.print_exc()

url.close()

Tracking file downloads with php

The situation: a server to which I have only user access.

What I want: to count the number of downloads, and to store the IP and hostname of the visitors (if I didn't want the latter, I would just use a pure html page with a refresh of zero seconds, for google analytics to capture).

There are a million solutions. Here is mine. It is not stealthy at all -users may workaround it easily. It is not meant to be an impediment; only a resonable attempt at tracking file downloads.

<?php
	# Usage: create a link to the file to trak like:
	# http://www.your.site.com/piper.php?file=some/file.tar.gz
	# 
	# Edit accordingly:

	$base_url = 'http://www.ini.uzh.ch/~acardona/';	


	# What this file does:
	# 1 - add +1 to $file."counter"
	# 2 - append address in $file."addresses"
	# 2 - refresh to actual file

	$file = $_GET['file'];
	#echo "file is ".$file;
	$counter_log = $file.".counter";

	if (!file_exists($counter_log)) {
		# create the file
		$handle = fopen($counter_log, 'w');
		chmod($counter_log, 777);
		$c = 1;
	} else {
		# open as read+write
		$handle = fopen($counter_log, 'r+');
		$count = fread($handle, filesize($counter_log));
		if (strlen($count) > 0) $c = intval($count) + 1;
		else $c = 1;
	}
	#echo "file opened<br />\n";
	if (!is_writable($counter_log)) {
		die("Can't write to ".$counter_log);
	}
	# needed to make jars folder 777 ... bofh php configuration# $%^&! Can't put the pointer at the beggining unless fopen again?
	fclose($handle);
	$handle = fopen($counter_log, 'w');
	#echo "file ".$counter_log." is writable<br />\n";	
	if (fwrite($handle, strval($c)) === FALSE) {
		die("Couldn't write to ".$counter_log);
	}
	#echo "wrote ".$c." to file ".$counter_log."<br />\n";
	fclose($handle);
	
	# log address:
	$addr_file = $file.".addresses";
	if (!file_exists($addr_file)) {
		$handle = fopen($addr_file, 'w');
		chmod($addr_file, 777);
	} else {
		$handle = fopen($addr_file, 'a');
	}
	$val = $_SERVER["REMOTE_HOST"];
	if (!$val) $val = $_SERVER["REMOTE_ADDR"];
	$addr = $val."		".date("Ymd G:i:s")."\n";
	fwrite($handle, $addr);
	fclose($handle);

	header("Location: ".$base_url.$file);
?>

Testing the famous Quake InvSqrt function

In C:

// Compare precision of 1/sqrt and Quake's InvSqrt function
//
// Conclusion:
//   - works when compiled with:
//       gcc -lm -Wall -o inverse_sqrt.o inverse_sqrt.c
//   - works with -O1 flag
//   - fails with -O2 or -O3 flags
//
//   Check manuals with man 3 sqrt or man 3 printf ... for C functions!
//
// Compile:
//    $ gcc -Wall -O1 -lm -o invert_sqrt.o invert_sqrt.c
//
// Run:
//    $ ./invert_sqrt.o 2.35
// 
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

float InvSqrt(float x) {
        float xhalf = 0.5f * x;
        int i = *(int*)&x; // store floating-point bits in integer
        i = 0x5f3759d5 - (i >> 1); // initial guess for Newton's method
        x = *(float*)&i; // convert new bits into float
        x = x*(1.5f - xhalf*x*x); // One round of Newton's method
        return x;
}

int main(int n_args, char* argsv[]) {
        if (n_args < 2) {
                printf("Usage: ./inverse_sqrt.o <float>\n");
                return 1;
        }
        float x = atof(argsv[1]);
        printf("Number: %f\nNormal inv sqrt: %f\nQuake's InvSqrt: %f\n", x, 1/sqrt(x), InvSqrt(x));

        return 0;
}

In java:

/** Results: it's a bit faster, like 5 to 10%, but varies somehow.
 *  Considering the loss of precision, IT'S NOT WORTH IT.
 *
 * Quake function:  0.70693004
 * Standard 1/sqrt: 0.7071067811865475
 *
 * Three runs of 100000000 times each:
 *
 * invsqrt: 165 ms
 * 1/sqrt: 172 ms
 * invsqrt: 159 ms
 * 1/sqrt: 172 ms
 * invsqrt: 172 ms
 * 1/sqrt: 172 ms
 *
 * Compile:
 *   $ javac -O QuakeInvSqrt.java
 * 
 * Run:
 *   $ java QuakeInvSqrt 2.0
 */

public class QuakeInvSqrt {

        private static final float invsqrt(float x) {
                final float xhalf = 0.5f * x;
                int i = Float.floatToRawIntBits(x);
                i = 0x5f3759d5 - (i >> 1);
                x = Float.intBitsToFloat(i);
                return x*(1.5f - xhalf*x*x);
        }

        static public void main(String[] arg) {
                if (arg.length < 1) return;
                float d;
                try {
                        d = Float.parseFloat(arg[0]);
                } catch (NumberFormatException nfe) {
                        p("Argument is not a number");
                        return;
                }
                p("Quake function:  " + invsqrt(d));
                p("Standard 1/sqrt: " + 1/Math.sqrt(d));

                for (int i=0; i<3; i++)
			test(d, 100000000);
        }

        static private void test(final float d, final float max) {
                long start = System.currentTimeMillis();
                for (int i=0; i<max; i++) {
                        double any = invsqrt(d);
                }
                p("invsqrt: " + (System.currentTimeMillis() - start) + " ms");

                start = System.currentTimeMillis();
                for (int i=0; i<max; i++) {
                        double any = 1/Math.sqrt(d);
                }
                p("1/sqrt: " + (System.currentTimeMillis() - start) + " ms");
        }                                                                     
                                                                              
        static private void p(String s) { System.out.println(s); }            
} 

Generating a plot with the first 10000 prime numbers

In C, but weaved into python.

Warning, this function may take several minutes to run.

This code takes advantage of the fact that all primes larger than 3 are next to a multiple of 6: for example, 11 and 13 and just -1 and +1 of 2*6.

# in C:

#int max = 100000;
#int pri[max];
#int pri[0] = 1;
#int pri[1] = 2;
#int pri[2] = 3;
#int next = 3;
#int i=0;
#int k=6;
#int j=0;
#inline void addIfPrime(const int a) {
#       for (j=0; j<next; j++) {
#               if (k-1 % pri[j]) return;
#       }
#       pri[next] = a;
#       next++;
#}
#while (i < max) {
#       addIfPrime(k-1);
#       addIfPrime(k+1);
#}

import scipy
from scipy.weave import converters
from scipy import zeros
import pylab

max = 10000
pri = zeros(max, int)
pri[0] = 1;
pri[1] = 2;
pri[2] = 3;

# numpy arrays are read as functions: pri[j] becomes pri(j)

code = """
        /*#include <stdio.h>*/
        bool b;
        int k = 6;
        int next = 3;
        int a = k;
        int j = 1;
        while (next < max) {
                a = k-1;
                b = true;
                for (j=1; j<next; j++) {
                        if (0 == a % pri(j)) { b=false; break; }
                }
                if (b) { pri(next) = a; /*printf("found: %i\n", a);*/ next++; }
                a = k+1;
                b = true;
                for (j=1; j<next; j++) {
                        if (0 == a % pri(j)) { b=false; break; }
                }
                if (b) { pri(next) = a; /*printf("found: %i\n", a);*/ next++; }
                k+=6;
        }
        return_val = 0;
        """
val = scipy.weave.inline(code,
             ['pri', 'max'],
             type_converters=converters.blitz,
             compiler = 'gcc')

pylab.plot(pri, max, 'ro')
pylab.show()

Generating a Manhattan-style city with Blender and python

Rather crude, but illustrates some basics of python scripting in Blender.

Besides randomly generated blocks neatly arranged in a grid of streets, within the city_night.blend file there is a reflectant material assigned to the streets and a texture simulating night-lighted windows for the buildings.

Yes, it can be done a lot better. This was just a fine pass time.

# Copyright Albert Cardona 2007
# This work is under the public domain.
# USE AT YOUR OWN RISK.

# REQUIRES Blender 2.44 or later.

from Blender import *
from random import randint

editmode = Window.EditMode()    # are we in edit mode?  If so ...
if editmode: Window.EditMode(0) # leave edit mode before creating the mesh

scene = Scene.GetCurrent()

def makeBlock(h):
        verts = [[0,0,0], [1,0,0], [1,1,0], [0,1,0], [0,0,h], [1,0,h], [1,1,h], [0,1,h]]
        faces = [[0, 1, 2, 3], [0, 1, 5, 4], [1, 2, 6, 5], [2, 3, 7, 6], [3, 0, 4, 7], [4, 5, 6, 7]]
        return verts, faces

def makeTile():
        verts = [[0,0,0], [1,0,0], [1,1,0], [0,1,0]]
        faces = [[0, 1, 2, 3]]
        return verts, faces

ob_streets = None
mesh_streets = None

for x in range(-30, 30):
        for y in range(-30, 30):
                if 0 == x % 4 or 0 == y % 6:
                        # make street tile
                        verts,faces = makeTile()
                        for v in verts:
                                v[0] += x
                                v[1] += y
                        if None is ob_streets:
                                mesh = Mesh.New('Streets')
                                mesh.verts.extend(verts)
                                mesh.faces.extend(faces)
                                ob_streets = scene.objects.new(mesh, 'streets')
                                mesh_streets = mesh
                        else:
                                # takes FOREVER # ob_streets.join([ob])
                                vl = len(mesh_streets.verts)
                                mesh_streets.verts.extend(verts)
                                for f in faces:
                                        f[0] += vl
                                        f[1] += vl
                                        f[2] += vl
                                        f[3] += vl
                                mesh_streets.faces.extend(faces)
                else:
                        h = randint(0, 5)
                        if 20 == randint(0, 20):
                                h += 4
                        verts,faces = makeBlock(h)
                        mesh = Mesh.New('Block')
                        mesh.verts.extend(verts)
                        mesh.faces.extend(faces)
                        ob = scene.objects.new(mesh, 'block')
                        ob.setLocation(x, y, 0)

mesh_streets.remDoubles(0)


if editmode: Window.EditMode(1)  # optional, just being nice
# Repaint the 3D window
Redraw()



Last updated: 2012-05-08 11:10 Zurich time. Copyright Albert Cardona.