Running statistics using RunningStats (how to run on Python using SWIG and Distutils)

There is a more accurate method for calculating focal statistics than usually implemented in numerical packages. This is documented here and the C++ code is download from here

The SWIG setup file is simple:

	/* File: RunningStats.i */
	%module RunningStats
	%{
	#define SWIG_FILE_WITH_INIT
	#include "RunningStats.h"
	%}
	%include "RunningStats.h"

view raw RunningStats.i hosted with ❤ by GitHub

which you run using:

swig -c++ -python RunningStats.i

view raw swig_runningstats.sh hosted with ❤ by GitHub

The Distutils setup script is:

	"""
	setup.py file for RunningStats
	"""
	from distutils.core import setup, Extension
	RunningStats_module = Extension('_RunningStats',
	sources=['RunningStats_wrap.cxx', 'RunningStats.cpp'],
	)
	setup (name = 'RunningStats',
	version = '0.1',
	author = "John D. Cook / Daniel Buscombe",
	description = """Compute running stats on an array""",
	ext_modules = [RunningStats_module],
	py_modules = ["RunningStats"],
	)

view raw setup.py hosted with ❤ by GitHub

compile using "python setup.py install"

Then in python, a simple implementation of calculating the global stats from an array would be something like:

	import numpy as np
	# generate some data
	x = np.random.rand(1000,1)
	# get the class
	import RunningStats
	# get an instance of the class
	rs = RunningStats.RunningStats()
	# global stats
	for k in x:
	rs.Push(k[0])
	print rs.Variance()
	print rs.Mean()
	print rs.Kurtosis()
	print rs.Skewness()
	print rs.StandardDeviation()
	rs.Clear()

view raw simple_stats.py hosted with ❤ by GitHub

 To show the convergence of locals means on the global mean, you might use the RunningStats class in the following way:

	import numpy as np
	import matplotlib.pylab as plt
	# generate some data
	x = np.random.rand(1000,1)
	# get the class
	import RunningStats
	# get an instance of the class
	rs = RunningStats.RunningStats()
	# convergence on global mean
	convav = []
	for k in x:
	rs.Push(k[0])
	convav.append(rs.Mean())
	plt.plot(x); plt.plot(convav,'r'); plt.show()

view raw conv_stats.py hosted with ❤ by GitHub

 Finally, to compute running averages on the array you might use something like the following, where the function to window the data efficiently was downloaded from here

	import numpy as np
	import matplotlib.pylab as plt
	# generate some data
	x = np.random.rand(1000,1)
	# get the class
	import RunningStats
	# get an instance of the class
	rs = RunningStats.RunningStats()
	# define a window size for the running average
	window_size = 10
	from itertools import product
	# use sliding window_nd from here: http://www.johnvinyard.com/blog/?p=268
	xs = sliding_window_nd(np.squeeze(x),window_size,1)
	runav = []
	for j in xrange(len(xs)):
	for k in xs[j]:
	rs.Push(k)
	runav.append(rs.Mean())
	rs.Clear()
	i = np.linspace(1,len(x),len(runav))
	plt.plot(x); plt.plot(i,runav,'r');
	# compute using a larger window size
	window_size = 50
	xs = sliding_window_nd(np.squeeze(x),window_size,1)
	runav = []
	for j in xrange(len(xs)):
	for k in xs[j]:
	rs.Push(k)
	runav.append(rs.Mean())
	rs.Clear()
	i = np.linspace(1,len(x),len(runav))
	plt.plot(i,runav,'k');
	plt.show()

view raw running_stats.py hosted with ❤ by GitHub

Obviously, as well as the mean, the above could be used for variance, standard deviation, kurtosis and skewness.

Pymix fix: short note on compiling PyMix on Ubuntu

A fix for users of Python 2.6 or 2.7 when installing PyMix:

	# 3 steps
	# 1) in setup.py
	# change
	from distutils.core import setup, Extension,DistutilsExecError
	#to:
	from distutils.core import setup, Extension
	from distutils.errors import DistutilsExecError
	# 2) in setup.py
	# change
	numpypath = prefix + '/lib/python' +pyvs + '/site-packages/numpy/core/include/numpy' # path to arrayobject.h
	# to
	# (locate arrayobject.h) - for me this path was
	numpypath = '/usr/share/pyshared/numpy/core/include/numpy'
	# Finally
	# 3) in /usr/local/lib/python2.7/dist-packages/AminoAcidPropertyPrior.py
	# change
	as = alpha.pop(0) # (line 170)
	# to
	dummy = alpha.pop(0) # 'as' is a reserved keyword

view raw pymix_fix.py hosted with ❤ by GitHub

How to pass arguments to python programs ...

... the right way (in my opinion). This is the basic template I always use. Not complicated, but the general formula might be useful to someone.

It allows you to pass arguments to a program with the use of flags, so they can be provided in any order. It has a 'help' flag (-h), will raise an error if no arguments are passed.

In this example one of the inputs (-a) is supposed to be a file. If the file is not provided, it uses Tkinter to open a simple 'get file' window dialog.

It then parses the other arguments out. Here I have assumed one argument is a string, one is a float, and one is an integer. If any arguments are not provided, it creates variables based on default variables. It prints everything it's doing to the screen.

	# these are the libraries needed
	import sys, getopt
	from Tkinter import Tk
	from tkFileDialog import askopenfilename
	import numpy as np
	# get list of input arguments and pre-allocate arrays
	argv = sys.argv[1:]
	inputfile = ''; arg2 = ''
	arg3 = ''; arg4 = ''
	# parse inputs to variables
	try:
	opts, args = getopt.getopt(argv,"ha:b:c:d:")
	except getopt.GetoptError:
	print 'my_script.py -a <inputfile> -b <arg2> -c <arg3> -d <arg4>'
	sys.exit(2)
	for opt, arg in opts:
	if opt == '-h': # h is 'help'
	print 'my_script.py -a <inputfile> -b <arg2> -c <arg3> -d <arg4>'
	sys.exit()
	elif opt in ("-a"):
	inputfile = arg
	elif opt in ("-b"):
	arg2 = arg
	elif opt in ("-c"):
	arg3 = arg
	elif opt in ("-d"):
	arg4 = arg
	# prompt user to supply file if no input file given
	if not inputfile:
	print 'An input file is required!!!!!!'
	# we don't want a full GUI, so keep the root window from appearing
	Tk().withdraw()
	# show an "Open" dialog box and return the path to the selected file
	inputfile = askopenfilename(filetypes=[("JPG files","*.jpg")])
	# print given arguments to screen and convert data type where necessary
	# let's assume arg2 is meant to be an integer, arg3 is a float, and arg4 is a string
	# let's further assume arg2 and arg3 are meant to be numpy arrays
	if inputfile:
	print 'Input file is %s' % (inputfile)
	if arg2:
	arg2 = np.asarray(arg2,int)
	print 'Argument 2 is %s' % (str(arg2))
	if arg3:
	arg3 = np.asarray(arg3,float)
	print 'Argument 3 is %s' % (str(arg2))
	if arg4:
	print 'Argument 4 is %s' % (arg4)
	# if arguments not supplied, use defaults
	if not arg2:
	arg2 = 100
	print '[Default] Argument 2 is %s' % (str(arg2))
	if not arg3:
	arg3 = 0.1
	print '[Default] Argument 3 is %s' % (str(arg3))
	if not arg4:
	arg4 = 'hello'
	print '[Default] Argument 4 is %s' % (arg4)

view raw starting_right.py hosted with ❤ by GitHub

Kivy python script for capturing displaying webcam

Kivy is very cool. It allows you to make graphical user interfaces for computers, tablets and smart phones ... in Python. Which is great for people like me. And you, I'm sure. And it's open source, of course.

Today was my first foray into programming using Kivy. Here was my first 'app'. It's not very sophisticated, but I find Kivy's documentation rather impenetrable and frustratingly lacking in examples. Therefore I hope more people like me post their code and projects on blogs.

Right, it detects and displays your webcam on the screen, draws a little button in the bottom. When you press that button it takes a screengrab and writes it to an image file. Voila! Maximise the window before you take a screen shot for best effects.

	import kivy
	kivy.require('1.7.2')
	from kivy.app import App
	from kivy.uix.floatlayout import FloatLayout
	from kivy.uix.camera import Camera
	from kivy.uix.button import Button
	from kivy.core.window import Window
	class CamApp(App):
	# Function to take a screenshot
	def screengrab(self,*largs):
	outname = self.fileprefix+'_%(counter)04d.png'
	Window.screenshot(name=outname)
	def build(self):
	# create a floating layout as base
	camlayout = FloatLayout(size=(600, 600))
	cam = Camera() #Get the camera
	cam=Camera(resolution=(1024,1024), size=(300,300))
	cam.play=True #Start the camera
	camlayout.add_widget(cam)
	button=Button(text='Take Picture',size_hint=(0.12,0.12))
	button.bind(on_press=self.screengrab)
	camlayout.add_widget(button) #Add button to Camera Layout
	self.fileprefix = 'snap'
	return camlayout
	if __name__ == '__main__':
	CamApp().run()

view raw simple_kivy_cam.py hosted with ❤ by GitHub

Simple Python script to write a list of coordinates to a .geojson file

This has probably been done better elsewhere, but here's a really simple way to write a geojson file programmatically to display a set of coordinates, without having to bother learning the syntax of yet another python library. Geojson is a cool little format, open source and really easy.

Below, arrays x and y are decimal longitude and latitude respectively. The code would have to be modified to include a description for each coordinate, which could be included easily on line 15 if in the form of an iterable list.

	f = open('points.geojson', 'w')
	f.write('\n')
	f.write('{\n')
	f.write(' "type": "FeatureCollection",\n')
	f.write(' "features": [\n')
	for k in range(len(x)):
	f.write('\n')
	f.write(' {\n')
	f.write(' "type": "Feature",\n')
	f.write(' "geometry": {\n')
	f.write(' "type": "Point",\n')
	f.write(' "coordinates": ['+str(y[k])+', '+str(x[k])+']\n')
	f.write(' },\n')
	f.write(' "properties": {\n')
	f.write(' "'+str(k)+'": "description here"\n')
	f.write(' }\n')
	if k==len(x)-1:
	f.write(' }\n')
	else:
	f.write(' },\n')
	f.write('\n')
	f.write('\n')
	f.write(' ]\n')
	f.write('}\n')
	f.close()

view raw xy2geojson.py hosted with ❤ by GitHub

Doing Spectral Analysis on Images with GMT

Mostly an excuse to play with, and show off the capabilities of, the rather fabulous Ipython notebook! I have created an Ipython notebook to demonstrate how to use GMT to do spectral analysis on images.

This is demonstrated by pulling wavelengths from images of rippled sandy sediments.

GMT is command-line mapping software for Unix/Linux systems. It's not really designed for spectral analysis in mind, especially not on images, but it is actually a really fast and efficient way of doing it! 

This also requires the convert facility provided by Image Magick. Another bodacious tool I couldn't function without.

My notebook can be found here.

Fantastic images of ripples provided by http://skeptic.smugmug.com

Taking the Grunt out of Geotags

Here's a matlab script to automatically compile all geo-tagged images from a directory and all its subdirectories into a Google Earth kml file showing positions and thumbnails

	close all;clear all;clc
	addpath('/path/where/googleearth/toolbox/is')
	% get main directory (Projects directory)
	pathFolder=uigetdir();
	% get list of subfolders
	d = dir(pathFolder);
	isub = [d(:).isdir];
	nameFolds = {d(isub).name}';
	nameFolds(ismember(nameFolds,{'.','..'})) = [];
	% prompt user to select which subfolders to process
	[s,v] = listdlg('PromptString','Select folders:',...
	'SelectionMode','multiple',...
	'ListString',nameFolds);
	% use only those directories
	directories=char(nameFolds(s));
	for md=1:size(directories,1)
	% get list of subfolders
	d = dir(strtrim([pathFolder,filesep,directories(md,:)]));
	isub = [d(:).isdir];
	nameFolds = {d(isub).name}';
	nameFolds(ismember(nameFolds,{'.','..'})) = [];
	subdirectories{md}=char(nameFolds);
	end
	mattime=[]; timestring=[]; lat=[]; lon=[]; filenames={};
	direc=[];
	counter=0;
	for dd=1:size(directories,1)
	for sd=1:size(subdirectories{dd},1)
	if isempty(strmatch(strtrim(subdirectories{dd}(sd,:)),'dng')) &&...
	isempty(strmatch(strtrim(subdirectories{dd}(sd,:)),'crw'))
	files1=ReadImDir(strtrim([pathFolder,filesep,strtrim(directories(dd,:)),...
	filesep,strtrim(subdirectories{dd}(sd,:))]),'JPG');
	files2=ReadImDir(strtrim([pathFolder,filesep,strtrim(directories(dd,:)),...
	filesep,strtrim(subdirectories{dd}(sd,:))]),'jpg');
	files=char(files1,files2); %,files3
	for k=1:size(files,1)
	in=strtrim([pathFolder,filesep,strtrim(directories(dd,:)),filesep,strtrim(subdirectories{dd}(sd,:)),...
	filesep,strtrim(files(k,:))]);
	if isempty(regexp(in,'ppm.jpg', 'once'))
	counter=counter+1;
	filenames{counter}=in;
	[stat,res]=system(['exiftool -h "',in,'"']);
	latstr=res(regexp(res,'>GPS Latitude<')+22:regexp(res,'>GPS Latitude<')+41);
	if ~isempty(latstr)
	deg=str2double(strtrim(latstr(1:2)));
	mins=str2double(strtrim(latstr(8:9)));
	secs=str2double(strtrim(latstr(16:end)));
	lat=[lat; deg+ mins/60 + secs/3600];
	lonstr=res(regexp(res,'>GPS Longitude<')+22:regexp(res,'>GPS Longitude<')+42);
	deg=str2double(strtrim(lonstr(2:4)));
	mins=str2double(strtrim(lonstr(10:11)));
	secs=str2double(strtrim(lonstr(17:end)));
	lon=[lon; deg+ mins/60 + secs/3600];
	timestr=res(regexp(res,'>GPS Date/Time<')+23:regexp(res,'>GPS Date/Time<')+41);
	y=str2double(timestr(1:4)); m=str2double(timestr(6:7)); d=str2double(timestr(9:10));
	H=str2double(timestr(12:13)); M=str2double(timestr(15:16)); S=str2double(timestr(18:19));
	if S==0
	S=S+1;
	end
	mattime=[mattime;datenum([y,m,d,H,M,S])];
	timestring=[timestring;datestr(datenum([y,m,d,H,M,S]))];
	direc=[direc;dd];
	else
	lat=[lat;NaN]; lon=[lon;NaN];
	timestr=res(regexp(res,'>Create Date<')+21:regexp(res,'>Create Date<')+39);
	y=str2double(timestr(1:4)); m=str2double(timestr(6:7)); d=str2double(timestr(9:10));
	H=str2double(timestr(12:13)); M=str2double(timestr(15:16)); S=str2double(timestr(18:19));
	if y==2012
	y=y+1;
	end
	if S==0
	S=S+1;
	end
	mattime=[mattime;datenum([y,m,d,H,M,S])];
	timestring=[timestring;datestr(datenum([y,m,d,H,M,S]))];
	direc=[direc;dd];
	end
	end % if not ppm.jpg
	end %k
	end % if not dng or crw
	end %sd
	end %dd
	thumbfiles=cell(size(filenames));
	for i=1:length(lat)
	thumbfiles{i}=[filenames{i},'_thumb.jpg'];
	end
	kmlStr=cell(1,length(lat));
	for i=1:length(lat)
	system(['exiftool -b -ThumbnailImage ',filenames{i},' > ',thumbfiles{i}])
	kmlStr{i} = ge_point_new(-lon(i),lat(i),0,...
	'description',['<a href="',filenames{i},'">link to file</a>'] ,...
	'iconURL',thumbfiles{i},...
	'name',filenames{i}(regexp(filenames{i},'IMG'):end));
	end
	tmp=[];
	for i=1:length(lat)
	tmp=[tmp,eval(['kmlStr{',num2str(i),'}'])];
	end
	ge_output('points.kml',...
	ge_folder('points',tmp));
	system(['google-earth "',pwd,filesep,'points.kml" &'])
	save('photo_res','mattime','lat','lon','direc','timestring','filenames','directories')

view raw geophotos2kml.m hosted with ❤ by GitHub

Like most of my matlab scripts, it's essentially just a wrapper for a system command. It requires exiftool for the thumbnail creation (available on both *nix and for Windows). It also needs the (free) Google Earth toolbox for Matlab.

Right now it's just set up for JPG/jpg files, but could easily be modified or extended. As well as the kml file (Google Earth will launch as soon as the script is finished) it will save a matlab format file contining the lats, longs, times and filenames.

Outputs look a little like this (a selection of over 30,000 photos taken in Grand Canyon)

A cooler shell prompt

In ~/.bashrc:

export PS1='\[\033[1;34m\]\t\[\033[0m\] \[\033[1;35m\]\u\[\033[0m\]:\[\033[1;35m\]\W\[\033[0m\] \[\033[1;92m\]$(__git_ps1 "(%s)")\[\033[0m\]$ '

view raw prompt.sh hosted with ❤ by GitHub

time: username: directory

much better!

How to start a BASH script

Use chance!

	# flip a coin
	FLIP=$(($(($RANDOM%10))%2))
	# if heads, use cowsay to present your splash
	if [ $FLIP -eq 1 ]
	then
	# start with some wisdom
	fortune -s | cowsay -n | zenity --text-info --title="Your message here" --width 500 --height 500
	else
	# or start with some other msg
	figlet DAN ROCKS | zenity --text-info --title="Your message here" --width 600 --height 200
	fi

view raw cow flip.sh hosted with ❤ by GitHub

Flip a coin. If heads, allow fortune to generate a quotation and cowsay to present it:

If tails, figlet your 'splash' message:

Discuss

Destination point given distance and bearing from start point

You have starting points:

xcoordinate, a list of x-coordinates (longitudes)

ycoordinate, a list of y-coordinates (latitudes)

num_samples, the number of samples in the plane towards the destination point

bearings:

heading, a list of headings (degrees)

and distances:

range, a list of distances in two directions (metres)

This is how you find the coordinates (x, y) along the plane defining a starting point and two end points, in MATLAB.

	R = 6371 % mean radius of the Earth, in km
	% pre-allocate for results
	x=cell(1,num_points)
	y=cell(1,num_points)
	for k=1:num_points
	% COORDINATES OF STARTING POINT
	lon1=xcoordinate(k) ;
	lat1=ycoordinate(k) ;
	d=range(k)/1000; %range in km
	b1=heading(k); % bearing in degrees
	% FIND FINAL BEARING FROM INITIAL BEARING
	if b1=>180
	b2 = b1-180;
	else
	b2 = b1 + 180;
	end
	% CONVERT TO RADIANS
	b1= b1.*(pi/180);
	b2= b2.*(pi/180);
	dOverR=d/R; % angular radians
	% LATITUDE, END POINT 1
	lat2 = asin( sin(lat1*(pi/180))*cos(dOverR) + ...
	cos(lat1*(pi/180))*sin(dOverR)*cos(b1) ) / (pi/180);
	% LONGITUDE, END POINT 1
	lon2 = lon1 + atan2( sin(b1)*sin(dOverR)*cos(lat1), ...
	cos(dOverR)-sin(lat1)*sin(lat2) ) ;
	% LATITUDE, END POINT 2
	lat3 = asin( sin(lat1*(pi/180))*cos(dOverR) + ...
	cos(lat1*(pi/180))*sin(dOverR)*cos(b2) ) / (pi/180);
	% LONGITUDE, END POINT 2
	lon3 = lon1 + atan2( sin(b2)*sin(dOverR)*cos(lat1), ...
	cos(dOverR)-sin(lat1)*sin(lat2) ) ;
	% FIND COORDINATES IN THAT PLANE
	y{k}=linspace( min([lat2,lat3]),max([lat2,lat3]), num_samples(k));
	x{k}=linspace( min([lon2,lon3]),max([lon2,lon3]), num_samples(k));
	end

view raw EndPoint_bearing_Distance.m hosted with ❤ by GitHub

Routine BASH memory aid

Some common things I do in BASH

	#Do something if a file exists
	if [ -f $target_dir"/$target_file" ]
	then do something here
	fi
	#Calculate pi using bc
	pi=`echo "4*a(1)" | bc -l`
	#Convert degs to radians
	rad=`echo "$deg*($pi/180)" | bc -l`
	#Calculate tangent in degrees
	tandeg=$(echo "s($rad)/c($rad)" | bc -l)
	#List subdirectories in pwd
	subdir=$(ls -d */)
	#Clear screen
	printf "\033c"
	#Calling python to do a math operation rather than bc , e.g. a ceil:
	round_up=$(python -c "from math import ceil; print ceil($variable)")
	#Pass arguments to script from command line
	args=("$@")
	first_arg=${args[0]}
	second_arg=${args[1]}
	#Pass number lines in file to variable
	numrecords=$(wc -l $target_dir"/file.ext" | awk -F" " '{print $1}')
	#Delete every file which is not a *.ext file
	find $target_dir -maxdepth 1 -type f ! -iname "*.ext" | xargs -I '{}' rm '{}'
	#Pipe the first two columns of file to new file
	awk '{print $1 " " $2 " " (-$3) }' $xyzfile".dat" > $xyzfile"new.dat"
	#Max and min of a column (n) in file
	n=2
	lims=$(awk 'BEGIN{ max = -999999999 ; min = 9999999999 } $n > max {max = $n} $n < min {min = $n} END{ print min, max }' $file)
	#Find a replace in file (NaNs with Os)
	sed -i 's/NaN/0/g' file.ext
	#Select a directory using zenity and pass to variable
	export my_dir=`zenity --file-selection --title="Select a Directory" --directory`

view raw bash memory aid.sh hosted with ❤ by GitHub

Coordinate system conversion using cs2cs from PROJ

I hate geodesy. It's confusing, multifarious, and ever-changing (and I call myself a scientist!). But we all need to know where we, and more importantly our data, are. I have need of a single tool which painlessly does all my conversions, preferably without me having to think too much about it.

You could use an online tool, but when you're scripting (and I'm ALWAYS scripting) you want something a little less manual. Enter cs2cs from the PROJ.4 initiative. It will convert pretty much anything to anything else. Perfect.

In Ubuntu:
sudo apt-get install proj-bin

In Fedora I searched yum for proj-4

If the man page is a little too much too soon to take in, this wonderful page helps you decide which flags and parameters to use, with a rudimentary understanding of what you want

For example, if I want to convert coordinates in Arizona Central State Plane to WGS84 Latitude/Longitude, I use this:

cs2cs -f "%.6f" +proj=tmerc +lat_0=31 +lon_0=-111.9166666666667 +k=0.9999 +x_0=213360 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs infile > outfile

view raw cs2cs 1.sh hosted with ❤ by GitHub

(the State Plane Coordinate System is new to me, but my view of the whole field of geodesy is 'let's make things even more complicated, sorry, accurate!')

Let's decompose that a bit:

-f "%.6f"

view raw proj1.sh hosted with ❤ by GitHub

tells it I want decimal degrees with 6 values after the decimal point (you could leave this out if you want deg, min, sec)

+proj=tmerc

view raw proj2.sh hosted with ❤ by GitHub

says use the Transverse Mercator projection

+lat_0=31 +lon_0=-111.9166666666667 +k=0.9999 +x_0=213360 +y_0=0 +ellps=GRS80 +datum=NAD83 +to_meter=0.3048006096012192 +no_defs

view raw proj3.sh hosted with ❤ by GitHub

are all parameters related to the input coordinates in Arizona Central State Plane (the 'to_meter' is a conversion from feet, the default unit, to metres)

+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs

view raw proj4.sh hosted with ❤ by GitHub

are the parameters related to the output coordinates in Latitude/Longitude

infile is a two column list of points e.g.
2.19644131000e+005 6.11961823000e+005
2.17676234764e+005 6.11243565478e+005
2.19763457634e+005 6.11234534908e+005
2.19786524782e+005 6.11923555789e+005
2.19762476867e+005 6.11378246389e+005

outfile will be created with the results:
-111.846501    36.517758 0.000000
-111.868478    36.511296 0.000000
-111.845175    36.511202 0.000000
-111.844912    36.517412 0.000000
-111.845185    36.512498 0.000000

By default it will always give you the third dimension (height) even if you didn't ask for it, like above, so I tend to trim it by asking awk to give me just the first 2 columns separated by a tab:

awk '{print $(NF-2),"\t",$(NF-1)}' outfile > outfile2

view raw awk trim.awk hosted with ❤ by GitHub

-111.846501      36.517758
-111.868478      36.511296
-111.845175      36.511202
-111.844912      36.517412
-111.845185      36.512498

Raspberry pi launch program on boot, and from desktop

The following is an example of how to launch a program running in a terminal upon boot, in Raspbian wheezy:

	sudo mkdir ~/.config/lxsession
	sudo mkdir ~/.config/lxsession/LXDE
	sudo nano ~/.config/lxsession/LXDE/autostart

view raw rpi1.sh hosted with ❤ by GitHub

and add:

lxterminal --geometry=50x10 --title="MY TITLE" -e python /usr/bin/my_script.py

view raw rpi2.sh hosted with ❤ by GitHub

This will allow the same program to be launched from the desktop upon double-click:

lxshortcut -o ~/Desktop/my_launcher.desktop

view raw rpi3.sh hosted with ❤ by GitHub

add

lxterminal --geometry=50x10 --title="MY TITLE" -e python /usr/bin/my_script.py

view raw rpi4.sh hosted with ❤ by GitHub

Compiling points2grid on Fedora with g++

Gridding large [x,y,z] point cloud datasets. Matlab and python have both failed me with their inefficient use of memory. So I'm having to turn to other tools.

Points2Grid is a tool, written in C++, which promises to do the job in a simple intuitive way. It accepts ascii and las (lidar) formats and outputs gridded data, plus gridded data densities.

First install LIBLAS which is Lidar data translation libraries. The recommended version is 1.2.1. I installed liblas 1.6.1 because it was the oldest version on the site that had a working link. I installed this by issuing:

	cd libLAS-1.6.1/
	mkdir makefiles
	cd makefiles
	cmake -G "Unix Makefiles" ../
	make
	sudo make install

view raw p2g1.sh hosted with ❤ by GitHub

Make a note of where the libraries are installed. On my system: /usr/local/include/liblas/

Download points2grid from here. 

In Interpolation.h, adjust MEM_LIMIT variable in the  file to reflect the memory limitations of your system. From the README guide, 

"As a rule of thumb, set the MEM_LIMIT to be (available memory in bytes)/55. For jobs that generate

grids cells greater than the MEM_LIMIT, points2grid goes "out-of-core", resulting in significantly worse performance."

Find this value using:

	mem=$(cat /proc/meminfo | grep MemTot | awk -F" " '{print $2}')
	echo "$mem/55" | bc

view raw p2g2.sh hosted with ❤ by GitHub

I then had to point the code to where the liblas directories had been installed. In Interpolation.cpp, replace:

	#replace
	#include <liblas/laspoint.hpp>
	#include <liblas/lasreader.hpp>
	#with
	#include </usr/local/include/liblas/liblas.hpp>
	#include </usr/local/include/liblas/reader.hpp>
	#replace
	liblas::LASHeader const& header = reader.GetHeader();
	#with
	liblas::Header const& header = reader.GetHeader();
	#replace
	liblas::Reader reader(ifs); //liblas::LASReader reader(ifs);
	#with
	liblas::LASReader reader(ifs);
	#replace
	liblas::LASPoint const& p = reader.GetPoint();
	#with
	liblas::Point const& p = reader.GetPoint();

view raw p2gg.sh hosted with ❤ by GitHub

It should then install with a simple 'make' command. Happy gridding!

Automatic Clustering of Geo-tagged Images. Part 2: Other Feature Metrics

In the last post I took a look at a simple method for trying to automatically cluster a set of images (of Horseshoe Bend, Glen Canyon, Arizona). Some of those picture were of the river (from different perspectives) and some of other things nearby.

With a goal of trying to separate the two classes, the results were reasonably satisfactory with image histograms as a feature metric, and a little user input. However, it'd be nice to get a more automated/objective way to separate the images into two discrete clusters. Perhaps the feature detection method requires a little more scrutiny?

In this post I look at 5 methods using the functionality of the excellent Mahotas toolbox for python downloaded here:

1) Image moments (I compute the first 5)
2) Haralick texture descriptors which are based on the co-occurrence matrix of the image
3) Zernike Moments
4) Threshold adjacency statistics (TAS)
5) Parameter-free threshold adjacency statistics (PFTAS)

The code is identical to the last post except for the feature extraction loop. What I'm looking for is as many images of the river at either low distances (extreme left in the dendrogram) or very high distances (extreme right in the dendrograms). These clusters have been highlighted below. There are 18 images of the river bend.

In order of success (low to high):

1) Zernike (using a radius of 50 pixels and 8 degrees, but other combinations tried):

	features = np.zeros([len(imlist),25])
	for i,f in enumerate(imlist):
	im = np.array(Image.open(f).convert('L'))
	h=mahotas.features.zernike_moments(im, 50, 8)
	features[i] = h.flatten()

view raw zernicke.py hosted with ❤ by GitHub

 This does a poor job, with no clusters at either end of the dendrogram.

2) Haralick:

	features = np.zeros([len(imlist),13])
	for i,f in enumerate(imlist):
	im = np.array(Image.open(f))
	h = mahotas.texture.haralick(im)[0]
	features[i] = h.flatten()

view raw haralick.py hosted with ❤ by GitHub

There is 1 cluster of 4 images at the start

3) Moments:

	features = np.zeros([len(imlist),5])
	for i,f in enumerate(imlist):
	im = np.array(Image.open(f).convert('L'))
	Z=np.array([],np.float)
	for k in range(5):
	Z=np.append(Z,mahotas.moments(im, k, k))
	features[i]=Z

view raw class moments.py hosted with ❤ by GitHub

There is 1 cluster of 6 images at the start, and 1 cluster of 2 at the end

4) TAS:

	features = np.zeros([len(imlist),162])
	for i,f in enumerate(imlist):
	im = np.array(Image.open(f))
	h=mahotas.features.tas(im)
	features[i] = h.flatten()

view raw tas class.py hosted with ❤ by GitHub

There is 1 cluster of 7 images at the start, and 1 cluster of 5 near the end, and both clusters are on the same stem

5) PFTAS:

	features = np.zeros([len(imlist),162])
	for i,f in enumerate(imlist):
	im = np.array(Image.open(f))
	h=mahotas.features.pftas(im)
	features[i] = h.flatten()

view raw pftas class.py hosted with ❤ by GitHub

The cluster at the end contains 16 out of 18 images of the river bend at the end. Plus there are no user-defined parameters. My kind of algorithm. The winner by far!

So there you have it, progress so far! No methods tested so far is perfect. There may or may not be a 'part 3' depending on my progress, or lack of!

Automatic Clustering of Geo-tagged Images. Part 1: using multi-dimensional histograms

There's a lot of geo-tagged images on the web. Sometimes the image coordinate is slightly wrong, or the scene isn't quite what you expect. It's therefore useful to have a way to automatically download images from a given place (specified by a coordinate) from the web, and automatically classify them according to content/scene.

In this post I describe a pythonic way to:
1) automatically download images based on input coordinate (lat, long)
2) extract a set features from each image
3) classify each image into groups
4) display the results as a dendrogram

This first example, 2) is achieved using a very simple means, namely the image histogram of image values. This doesn't take into account any texture similarities or connected components, etc. Nonetheless, it does a reasonably good job at classifying the images into a number of connected groups, as we shall see.

In subsequent posts I'm going to play with different ways to cluster images based on similarity, so watch this space!

User inputs: give a lat and long of a location, a path where you want the geo-tagged photos, and the number of images to download

	#I tried this out on Horseshoe Bend, a much-photographed site in Glen Canyon, AZ. The approximate coordinates are
	lt=36.879444
	ln= -111.513889
	#I wanted to limit the analysis to 37 images:
	numimages=37
	#Here I'm just setting the path to be the present working directory:
	path=os.getcwd()
	#Get the images from panoramio. I've hard-wired it to search for images within lt + lt+0.1 and ln + ln+0.01, but you can specify any search radius you like.
	url = "http://www.panoramio.com/map/get_panoramas.php?order=popularity&set=public&from=0&to=%s&minx=%s&miny=%s&maxx=%s&maxy=%s&size=medium" % (str(numimages),str(ln),str(lt),str(ln+.01),str(lt+.1))

view raw geotag.py hosted with ❤ by GitHub

First, import the libraries you'll need. Then interrogates the website and downloads the images.

	try: import simplejson as json
	except ImportError: import json
	import os
	import urllib, urlparse
	import glob
	import numpy as np
	import Image
	from itertools import combinations
	import simplejson as json
	import pylab as mpl
	c = urllib.urlopen(url)
	j =json.loads(c.read())
	imurls=[]
	for im in j['photos']:
	imurls.append(im['photo_file_url'])
	for url in imurls:
	image = urllib.URLopener()
	image.retrieve(url,os.path.basename(urlparse.urlparse(url).path))

view raw auto-class 1.py hosted with ❤ by GitHub

As you can see the resulting images are a mixed bag. There's images of the river bend, the road, the desert, Lake Powell and other random stuff. My task is to automatically classify these images so it's easier to pull out just the images of the river bend

The following compiles a list of these images. The last part is to sort which is not necessary but doing so converts the list into a numpy array which is.  Then clustering of the images is achieved using Jan Erik Solem's rather wonderful book 'Programming Computer Vision with Python'. The examples from which can be downloaded here. Download this one, then this bit of code does the clustering:

	imlist=[]
	for infile in glob.glob( os.path.join(path, '*.jpg') ):
	imlist.append(infile)
	imlist=np.sort(imlist)
	features = np.zeros([len(imlist),512])
	for i,f in enumerate(imlist):
	im = np.array(Image.open(f))
	h, edges = np.histogramdd(im.reshape(-1,3),8,normed=True, range=[(0,255),(0,255),(0,255)])
	features[i] = h.flatten()
	import hcluster
	tree = hcluster(features)

view raw auto-class 2.py hosted with ❤ by GitHub

The approach taken is to use hierarchical clustering using a simple euclidean distance function. This bit of code does the dendogram plot of the images:

draw_dendrogram(tree,imlist,filename='horseshoe.png')

view raw dendro.py hosted with ❤ by GitHub

which in my case looks like this (click on the image to see the full extent):

It's a little small, but you can just about see (if you scroll to the right) that it does a good job at clustering images of the river bend which look similar. The single-clustered, or really un-clustered, images to the left are those of the rim, road, side walls, etc which don't look anything like the river bend. 

Next, reorder the image list in terms of similarity distance (which increases in both the x and y directions of the dendrogram above)

	order = tree.get_cluster_elements()
	slist=imlist[order]
	#If i'm only interested in images of the river bend, you can see at a glance that I can remove the final cluster of 4 images, and then take the last 2 sub-clusters (comprising 10 images) so this gives me those images:
	a=4
	slist=slist[:len(order)-a]
	a=10
	slist=slist[-a:]
	#Finally, I create a plot of only those final 2 clusters
	mpl.figure()
	for i in range(10):
	im=Image.open(slist[i])
	mpl.subplot(2,5,i+1)
	mpl.imshow(np.array(im))
	mpl.axis('equal')
	mpl.axis('off')
	mpl.savefig('riverbend.png')
	mpl.close()

view raw geotag2.py hosted with ❤ by GitHub

Which gives me:

As you can see they are all images of the river bend, except the 2nd from the left on the top row, which is a picture of a shrub. Interestingly, the pattern of a circular shrub surrounded by a strip of sand is visually similar to the horse shoe bend!!

However, we don't want to include it with images of the river, which is why a more sophisticated method that image histogram is required to classify and group similar image ... the subject of a later post.

Alpha Shapes in Python

Alpha shapes include convex and concave hulls. Convex hull algorithms are ten a penny, so what we're really interested in here in the concave hull of an irregularly or otherwise non-convex shaped 2d point cloud, which by all accounts is more difficult.

The function is here:

	def get_alpha_shape(infile,radius):
	command = "%s -A -aa %s -r -m10000000 -oN -oFpoints < %s" % (hull_path, str(radius), infile)
	print >> sys.stderr, "Running command: %s" % command
	retcode = subprocess.call(command, shell=True)
	results_file = open("points-alf")
	results_file.next()
	results_indices = [[int(i) for i in line.rstrip().split()] for line in results_file]
	results_file.close()
	return results_indices

view raw alpha.py hosted with ❤ by GitHub

The above is essentially the same wrapper as posted here, except a different way of reading in the data, and providing the option to specify a probe radius. It uses Ken Clarkson's C-code, instructions on how to compile are here

An example implementation:

	t = linspace(0.6,5.7,500)
	C = 2*np.vstack([cos(t),sin(t)] )
	C=C+np.random.rand(2,500)

view raw alpha2.py hosted with ❤ by GitHub

The above data generation is translated from the example in a matlab alpha shape function . Then:

	hull_path = "./hull"
	infile='C.txt'
	with open(infile, 'w') as f:
	np.savetxt(f, C.T, delimiter=' ', fmt="%0.7f %0.7f\n")
	radius=1
	indices = get_alpha_shape(infile,radius)
	alpha=C.T[indices]
	plot(C[0],C[1],'.'), plot(alpha[1:].T[0],alpha[1:].T[1],'r')
	savefig("C_radius="+str(radius)+".png")
	close()

view raw using alpha.py hosted with ❤ by GitHub

Which produces (points are blue dots, alpha shape is red line):

and on a larger data set (with radius=1, this is essentially the convex hull):

and a true concave hull:

Patchwork quilt in Python

Here's a fun python function to make a plot of [x,y] coordinate data as a patchwork quilt with user-defined 'complexity'. It initially segments the data into 'numsegs' number of clusters by a k-means algorithm. It then takes each segment and creates 'numclass' sub-clusters based on the euclidean distance from the centroid of the cluster. Finally, it plots each subcluster a different colour and prints the result as a png file.

inputs:
'data' is a Nx2 numpy array of [x,y] points
'numsegs' and 'numclass' are integer scalars. The greater these numbers the greater the 'complexity' of the output and the longer the processing time

	import numpy as np
	from scipy.cluster.vq import kmeans,vq
	from pylab import *
	def patchwork(data,numsegs,numclass):
	# computing K-Means with K = numsegs
	centroids,_ = kmeans(data,numsegs)
	# assign each sample to a cluster
	idx,_ = vq(data,centroids)
	# loop through number of classes
	for k in range(numsegs):
	# get data in kth segment
	datXY=data[idx==k,:]
	# get distances of each point from centroid
	dat1=np.sqrt((datXY[:,0]-centroids[k,0])**2 + (datXY[:,1]-centroids[k,1])**2)
	# get numclass new centroids within kth segment
	cents,_ = kmeans(dat1,numclass)
	# assign each to subcluster
	i,_ = vq(dat1,cents)
	# loop through and plot each
	for p in range(numclass):
	plot(datXY[i==p,0],datXY[i==p,1],'.')
	savefig("outputs"+str(numsegs)+'_'+str(numclass)+'_res.png')
	close()

view raw patchwork.py hosted with ❤ by GitHub

Some example outputs in increasing complexity:

numsegs=10, numclass=5

numsegs=15, numclass=15

numsegs=20, numclass=50

Compiling Clarkson's Hull on Fedora

A note on compiling Ken Clarkson's hull program for efficiently computing convex and concave hulls (alpha shapes) of point clouds, on newer Fedora

First, follow the fix described here which fixes the pasting errors given if using gcc compiler.

Then, go into hullmain.c and replace the line which says:

	FILE *INFILE, *OUTFILE, *DFILE = stderr, *TFILE;
	//to read:
	FILE *INFILE, *OUTFILE, *TFILE, *DFILE;
	//and in main() add the following lines
	FILE *DFILE;
	DFILE = stderr;

view raw hullcompile.c hosted with ❤ by GitHub

at the end of the function declarations (before the first while loop)

Finally, rewrite the makefile so it reads as below It differs significantly from the one provided.

	CC = gcc
	AR = ar
	OBJS = hull.o ch.o io.o rand.o pointops.o fg.o
	HDRS = hull.h points.h pointsites.h stormacs.h
	SRC = hull.c ch.c io.c rand.c pointops.c fg.c
	PROG = hull
	BINDIR = /usr/bin
	LIBDIR = /usr/lib64 # or just /usr/lib on a 32 bit machine?
	LIB = $(LIBDIR)/lib$(PROG).a
	all : $(PROG) rsites
	cp $(PROG) $(BINDIR)/.
	cp rsites $(BINDIR)/.
	$(OBJS) : $(HDRS)
	hullmain.o : $(HDRS)
	$(PROG) : $(OBJS) hullmain.o
	$(CC) $(OBJS) hullmain.o -o $(PROG) -lm
	$(AR) rcv $(LIB) $(OBJS)
	rsites : rsites.c
	$(CC) -o rsites rsites.c -lm
	clean :
	-rm -f $(OBJS) hullmain.o core a.out $(PROG)

view raw makefile hosted with ❤ by GitHub

Notice how the CFLAGS have been removed. Compile as root (sudo make)

Create Custom Launchers in Fedora

I just finished installed SmartGit on my new OS, like I did the last one

SmartGit is a simple bash script that I had to run from the command line. This time I decided I wanted a program launcher from my applications menu and my 'favourites bar' (or whatever that thing on the left hand side of the desktop is called!) 

So first I moved the smartgit folder from ~/Downloads/ to /usr/local/ and made the .sh script executable

	sudo cp -R ~/Downloads/smartgit-3_0_11/ /usr/local/
	chmod +x /usr/local/smartgit-3_0_11/bin/smartgit.sh

view raw clf1.sh hosted with ❤ by GitHub

Then to create a launcher I made a new .desktop file:

sudo touch /usr/local/share/applications/smartgit.desktop

view raw clf2.sh hosted with ❤ by GitHub

and added the following:

	Encoding=UTF-8
	Version=1.0
	Type=Application
	Exec=/usr/local/smartgit-3_0_11/bin/smartgit.sh
	Name=SMARTGIT
	Comment=custom launcher for smartgit.sh
	NoDisplay=false
	Categories=Application;Programming;Utility;
	Icon=/usr/local/smartgit-3_0_11/bin/smartgit-128.png
	[Desktop Entry]

view raw clf3.sh hosted with ❤ by GitHub

Of course the 'Icon' can be anything you like - I chose the one which comes with the program out of sheer non-inventiveness! In Applications it will appear under All and Accessories. Right click and 'Add to Favourites' and it will appear on your side bar.

Create a diary for your daily workflow

You can modify your ~./bash_history file so it doens't store duplicates and stores more lines than the default (only 1000):

	export HISTCONTROL=erasedumps
	export HISTSIZE=100000
	shopt -s histappend

view raw diary1.sh hosted with ❤ by GitHub

However, sometimes you want to go back to your history file and see what commands you issued on a particular day .... one solution is to make a periodic backup with a timestamp using cron.

First, make a directory to put your files in, and make a script which is going to do the work:

	cd
	touch backup_history.sh
	mkdir history

view raw diary2.sh hosted with ❤ by GitHub

Open the shell script and copy the following inside

	nano backup_history.sh
	DATE=$(date +"%Y%m%d%H%M")
	cp ~/.bash_history ~/history/$DATE.log

view raw diary3.sh hosted with ❤ by GitHub

close, then in the terminal issue

	sudo nano /etc/crontab
	#and copy the following
	@hourly ~/./backup_history.sh

view raw cron.sh hosted with ❤ by GitHub

admittedly hourly is a little excessive - go here is you are unfamiliar with cron scheduling definitions and modify as you see fit

Installing Mb-System on Fedora 17

MB-System is the only fully open-source software I am aware of for swath sonar data including multibeam and sidescan backscatter. As such, I'm really hoping using it works out for me!

Compiling it on a fresh install of Fedora 17 wasn't that straight-forward (took me about 5 hours to figure all this out), so I thought I'd put together a little how-to for the benefit of myself, and hopefully you.

I've broken it down into a few stages so some sort of semblance of order can be imposed on this random collection of installation notes.

Step 1: download the software from ftp server here

In the terminal:

	cd ~/Downloads
	tar xvzf MB-System.tar.gz

view raw mb1.sh hosted with ❤ by GitHub

For me this created a folder called 'mbsystem-5.3.1982'

Step 2: install pre-requisites

A) Generic Mapping Tools (GMT). You need to go to the Downloads tab, then click on the INSTALL FORM link. This will take you to a page where you can download the installation script, and fill in an online form of your parameter settings which you can submit and save the resulting page as a textfile which becomes the input to the bash script. Sounds confusing, but it's not - the instructions on the page are adequate. I opted to allow GMT to install netCDF for me. Then in the terminal I did this:

	sudo mkdir /usr/local/GMT
	sudo sh install_gmt.sh GMT[45]param.txt
	export NETCDFHOME=/home/me/Downloads/netcdf-3.6.3
	export PATH=/usr/local/GMT/bin:$PATH

view raw mb2.sh hosted with ❤ by GitHub

B) Perl parallel forkmanager:

sudo yum install perl-Parallel-ForkManager.noarch

view raw mb3.sh hosted with ❤ by GitHub

C) X11

sudo yum groupinstall "X Software Development"

view raw mb4.sh hosted with ❤ by GitHub

D) openmotif

sudo yum install openmotif-clients.x86_64 openmotif.x86_64 openmotif-demos.x86_64

view raw mb5.sh hosted with ❤ by GitHub

E) fftw

sudo yum install fftw.x86_64 fftw-devel.x86_64

view raw mb6.sh hosted with ❤ by GitHub

F) ghostview - I had to install this indirectly using kde3:

sudo yum install kdebase3.x86_64

view raw mb6b.sh hosted with ❤ by GitHub

(Note to Fedora - why oh why oh why are B, C, and F above not installed by default!?)

G)  OTPSnc tidal prediction software: download from here

untar, and cd to the directory

first double check that ncdump and ncgen are installed (which ncdump ncgen)

then edit the makefile so it reads:

	FC = gfortran
	NCLIB = /usr/lib64
	NCINCLUDE = /usr/include
	NCLIBS= -lnetcdf -lnetcdff

view raw mb7.sh hosted with ❤ by GitHub

then in the terminal issue:

	make extract_HC
	make extract_predict_tide

view raw mb8.sh hosted with ❤ by GitHub

Hopefully this compiles without errors, then I moved them to a executable directory:

	sudo mv extract_HC /usr/bin/
	sudo mv predict_tide /usr/bin/

view raw mb9.sh hosted with ❤ by GitHub

H) openGL libraries (again I had to do this indirectly, this time through mesa):

sudo yum install mesa*

view raw mb10.sh hosted with ❤ by GitHub

Step 3: prepare mbsystem makefiles

cd mbsystem-5.3.1982/

You have to go in and point install_makefiles to where all your libraries are. This is time-consuming and involves a lot of ls, which, and whereis!

Here's a copy of the lines I edited in my install_makefiles parameters:

	$MBSYSTEM_HOME = "/home/me/Downloads/mbsystem-5.3.1982";
	$OS = "LINUX";
	$CFLAGS = "-Wall -g -I/usr/include/X11";
	$LFLAGS = "-Wall -lm -bind_at_load -Wl,-dylib_file,/usr/lib64/libGL.so";
	$NETCDFLIBDIR = "/usr/local/lib";
	$NETCDFINCDIR = "/usr/include";
	$GMTLIBDIR = "/usr/local/GMT/lib";
	$GMTINCDIR = "/usr/local/GMT/include";
	$LEVITUS = "$MBSYSTEM_HOME/share/annual";
	$PROJECTIONS = "$MBSYSTEM_HOME/share/Projections.dat";
	$FFTWLIBDIR = "/usr/lib64";
	$FFTWINCDIR = "/usr/include";
	$MOTIFINCDIR = "/usr/include/openmotif/";
	$MOTIFLIBS = "-L/usr/lib64/openmotif -L/usr/lib64/X11 -lXm -lXt -lX11";
	$LEVITUS = "$MBSYSTEM_HOME/share/annual";
	$OTPSDIR = "/usr/bin";
	$CC = "gcc";

view raw mb12.sh hosted with ❤ by GitHub

Then in the terminal

Step 3: install mbsystem

first I had to do this (but you may not need to)

	sudo cp /usr/lib64/libfftw3.so /usr/lib64/libfftw3.a
	#then
	make all
	#assuming all worked ... I copied the man and bin files into somewhere more logical/permanent:
	sudo mkdir /usr/local/mbsystem
	sudo mkdir /usr/local/mbsystem/bin
	sudo cp ~/Downloads/mbsystem-5.3.1982/bin/* /usr/local/mbsystem/bin
	sudo cp -R ~/Downloads/mbsystem-5.3.1982/man/ /usr/local/mbsystem/

view raw mb14.sh hosted with ❤ by GitHub

I then updated my ~/.bashrc so the computer can find all these lovely new files:

	MANPATH=$MANPATH:/usr/local/GMT/man:/usr/local/mbsystem/man
	export MANPATH
	PATH=$PATH:/usr/local/GMT/bin:/usr/local/mbsystem/bin
	export PATH

view raw mb15.sh hosted with ❤ by GitHub

Test the install by issuing

	man gmt
	man mbfilter

view raw mb16.sh hosted with ❤ by GitHub

you should see the manual pages come up. Right, now down to the business of multibeam analysis!

em1 to ethX on Fedora 17 (how to install matlab if you don't have eth)

I have just moved from Ubuntu 12.04 to Fedora 17. No reason other than necessity for my new job. I'm still happy with Ubuntu and even don't mind Unity!

First thing I wanted to do, of course, was install matlab. It turns out I had a bit of work to do to get it installed on my Dell Precision machine! So here's a post to help out anyone who has the same problem.

If you don't have some eth listed when you ifconfig (I had the NIC listed as em1), matlab will install but you can't activate it. The problem is described on their website  

Basically, Fedora udev names NICs (network cards) embedded on the motherboard as em1, em2 (etc) rather then eth0, eth1 (etc). More details here

The two solutions listed on the Mathworks for changing the names back to eth do not work in my Fedora install. The reason why the first one doesn't work is because there is no /etc/iftab file. the reason why the second one works is because udev rules are not listed in etc (they are in lib instead and have different names and structures)

Thankfully, Fedora has built in a mechanism to modify the names of NICs post install. You do it by editing the following file

/lib/udev/rules.d/71-biosdevname.rules

so the line which reads

GOTO="netdevicename_end"

commented by default, is uncommented

Then you need to edit
/etc/default/grub

by adding the following line to the bottom:

biosdevname=0

Then I backed up my ifcfg em1 network script by issuing

cp /etc/sysconfig/network-scripts/ifcfg-em1 /etc/sysconfig/network-scripts/ifcfg-em1.old

then I renamed it eth0 (the number can be anything, I believe, between 0 and 9, and matlab will figure it out when it activates):

mv /etc/sysconfig/network-scripts/ifcfg-em1 /etc/sysconfig/network-scripts/ifcfg-eth0

and edited ifcfg-eth0 so the line which read

DEVICE="em1"

to read:

DEVICE="eth0"

This is all the information the machine needs, when it boots up, to change the naming convention from em1 back over to eth0.

Upon reboot, ifconfig shows eth0 where em1 was before (dmesg | grep eth would also show the change), and matlab installs perfectly.

A satisfyingly geeky end to a Tuesday!

Downloading Photobucket albums with Matlab and wget

Here's a script to take the pain out of downloading an album of photos from photobucket.

It uses matlab to search a particular url for photo links (actually it looks for the thumbnail links and modifies the string to give you the full res photo link), then uses a system command call to wget to download the photos All in 8 lines of code! 

Ok to start you need the full url for your album. The '?start=all' bit at the end is important because this will show you all the thumbnails on a single page. Then simply use regexp to find the instances of thumbnail urls. Then loop through each link, strip the unimportant bits out, remove the 'th_' from the string which indicates the thumbnail version of the photo you want, and call wget to download the photo. Easy-peasy!

	s = urlread('http://s1270.photobucket.com/albums/SERVERNAME/YOURUSERNAME/YOURALBUMNAME/?start=all');
	indices = regexp(s,'pbthumburl=');
	for k=1:length(indices)
	tmp=strtrim(s(indices(k):indices(k)+100));
	tmp=tmp(regexp(tmp,'"','once')+1:end);
	tmp=tmp(1:regexp(tmp,'"','once')-1);
	system(['wget ',tmp(1:regexp(tmp,'th_')-1),tmp(regexp(tmp,'th_')+3:end)]);
	end

view raw photobucket.m hosted with ❤ by GitHub

Alternative to bwconvhull

Matlab's bwconvhull is a fantastic new addition to the image processing library. It is a function which, in its full capacity, will return either the convex hull of the entire suite of objects in a binary image, or alternatively the convex hull of each individual object. I use it, in the latter capacity, for filling multiple large holes (holes large enough that imfill has no effect) in binary images containing segmented objects 

For those who do not have a newer Matlab version with bwconvhull included, I have written this short function designed to be an alternative to the usage

 P= bwconvhull(BW,'objects');

and here it is:

	function P = bwconvhull_alt(BW)
	% usage: P is a binary image wherein the convex hull of objects are returned, BW is the input binary image
	% P= bwconvhull_alt(BW);
	warning off all
	s=regionprops(logical(BW),'ConvexImage','BoundingBox');
	P=zeros(size(BW));
	for no=1:length(s)
	P(s(no).BoundingBox(2):s(no).BoundingBox(2)+s(no).BoundingBox(4)-1,...
	s(no).BoundingBox:s(no).BoundingBox(1)+s(no).BoundingBox(3)-1)=s(no).ConvexImage;
	end
	warning on all

view raw bwconvhull_alt.m hosted with ❤ by GitHub

It uses the IP toolbox function regionprops, around for some time, to find the convex hull images of objects present, then inscribes them onto the output image using the bounding box coordinates. Simples!

Craigslist Part 2: using matlab to plot your search results in Google Maps

In 'Craigslist Part 1' I demonstrated a way to use matlab to automatically search craigslist for properties with certain attributes. In this post I show how to use matlab and python to create a kml file for plotting the results in Google Earth. Download the matlab googleearth toolbox from here Add google earth toolbox to path. Import search results data and get map links (location strings). Loop through each map location and string the address from the url, then geocode to obtain coordinates.

	addpath(genpath('~/googleearth'))
	dat=importdata('myFile.csv');
	maps=dat.textdata(:,2);
	prices=dat.data;
	for i=1:size(maps,1)
	address=char(maps(i,:));
	address=address(regexp(address,'loc')+3:end);
	% strip address of unnneccessary stuff
	tmp = strrep(address, '%', ' ');
	tmp = strrep(tmp, '+', ' ');
	tmp = strrep(tmp, '+', ' ');
	tmp(isstrprop(tmp, 'digit'))=[];
	space=isspace(tmp);
	letters=find(isletter(tmp)==1);
	letters = letters(floor(gradient(letters))==1);
	space(letters)=1;
	% trim single letters but preserve location of spaces to make address
	% readable
	clean_address=strtrim(tmp(space));
	% geocode address to lat,long coordinate. Geocode generates and runs a
	% python script to do the conversion
	[g1,g2] = geocode(clean_address);
	% if any empty or unconverted addresses, call them nans
	if isempty(g1)
	g1=NaN; g2=NaN;
	end
	lat(i)=g1; lon(i)=g2;
	end

view raw craig2.m hosted with ❤ by GitHub

The function 'geocode' writes and executes a python script to geocode the addresses (turn the address strings into longitude and latitude coordinates). The python module may be downloaded here

	function [lat,lon] = geocode(address)
	fid=fopen('test.py','wt');
	fprintf(fid,'%s\n','from pygeocoder import Geocoder');
	fprintf(fid,'%s\n',['results = Geocoder.geocode(''',deblank(address),''')']);
	fprintf(fid,'%s\n','print results[0].coordinates');
	fclose(fid);
	[stat1,stat2]=system('python test.py');
	stat2=deblank(stat2);
	stat2=stat2(2:end-1);
	lon=str2num(deblank(stat2(regexp(stat2,',')+1:end)));
	lat=str2num(deblank(stat2(1:regexp(stat2,',')-1)));

view raw geocode.m hosted with ❤ by GitHub

Once we have the coordinates, we then need to get rid of nans and outliers (badly converted coordinates due to unreadable address strings). Use the google earth toolbox to build the kml file. Finally, run google earth and open the kml file using a system command:

	f=isnan(lat);
	lat(f)=[];
	lon(f)=[];
	maps(f)=[];
	prices(f)=[];
	f=find(lat>(mean(lat)+2*std(lat)));
	lat(f)=[]; lon(f)=[]; maps(f)=[]; prices(f)=[];
	f=find(lat<(mean(lat)-2*std(lat)));
	lat(f)=[]; lon(f)=[]; maps(f)=[]; prices(f)=[];
	f=find(lon>(mean(lon)+2*std(lon)));
	lat(f)=[]; lon(f)=[]; maps(f)=[]; prices(f)=[];
	f=find(lon<(mean(lon)-2*std(lon)));
	lat(f)=[]; lon(f)=[]; maps(f)=[]; prices(f)=[];
	kmlStr = ge_scatter(lon(:),lat(:),...
	'marker','*',...
	'markerEdgeColor','FFFF00FF',...
	'markerFaceColor','80FF00FF',...
	'markerScale',1e-3);
	ge_output('scatter.kml',[kmlStr])
	system(['google-earth "',pwd,filesep,'scatter.kml" &'])

view raw craig2b.m hosted with ❤ by GitHub

The above is a simple scatter plot which only shows the location of the properties and not any information about them. Next shows a more complicated example where the points are plotted with labels (the asking price) and text details (the google map links) in pop-up boxes First each coordinate pair is packaged with the map and name tags. Concatenate the strings for each coordinate and make a kml file. Finally, run google earth and open the kml file using a system command:

	iconStr = 'http://maps.google.com/mapfiles/kml/pal2/icon10.png';
	kmlStr=cell(1,length(maps));
	for i=1:length(maps)
	kmlStr{i} = ge_point_new(lon(i),lat(i),0,...
	'iconURL',iconStr,...
	'iconColor','FF0080FF',...
	'description',maps{i},...
	'name',['$',num2str(prices(i))]);
	end
	tmp=[];
	for i=1:length(maps)
	tmp=[tmp,eval(['kmlStr{',num2str(i),'}'])];
	end
	ge_output('points.kml',...
	ge_folder('points',tmp));
	system(['google-earth "',pwd,filesep,'points.kml" &'])

view raw craig2c.m hosted with ❤ by GitHub

Broadcom wireless on Ubuntu 12.04

I just upgraded to Ubuntu 12.04 and my broadcom wireless card was not being picked up. It took me some time to figure out how to do this, and the forums give a lot of tips which didn't work for me, so I thought I'd give it a quick mention.

	#In the terminal type
	sudo rfkill unblock all
	#then if you type
	rfkill list all
	#It should show something like:
	1: brcmwl-0: Wireless LAN
	    Soft blocked: no
	    Hard blocked: no
	#then type
	echo "blacklist hp-wmi" | sudo tee -a /etc/modprobe.d/blacklist.conf

view raw broadcom.sh hosted with ❤ by GitHub

Then disconnect your wired connection and reboot. I found it is crucial that you disconnect any wired connections.

 Upon reboot, go to System Settings - Hardware - Additional Drivers

 It should then pick up the broadcom proprietary driver. Install it, then reboot, and you should be back up and running!

Craigslist Part 1: Using matlab to search according to your criteria and retrieve posting urls

I have a certain fondness for using matlab for obtaining data from websites. It simply involves understanding the way in which urls are constructed and a healthy does of regexp. In this example, I've written a matlab script to search craiglist for rental properties in Flagstaff with the following criteria: 1) 2 bedrooms 2) has pictures 3) within the price range 500 to 1250 dollars a month 4) allows pets 5) are not part of a 'community' housing development. This script will get their urls, prices, and google maps links, and write the results to a csv file. First, get the url and find indices of the hyperlinks only, then get only the indices of urls which are properties. Then get only the indices of urls which have a price in the advert description; find the prices; strip the advert urls from the string; sort the prices lowest to highest; and get rid of urls which mention a community housing development. Finally, get the links to google maps in the adverts if there are some, and write the filtered urls, maps, and prices to csv file.

	s = urlread('http://flagstaff.craigslist.org/search/apa?query=&srchType=T&minAsk=500&maxAsk=1250&bedrooms=2&addThree=wooof&hasPic=1');
	% find the hyperlinks
	ind=regexp(s,'href');
	useind=[];
	for k=1:length(ind)
	if ~isempty(regexp(s(ind(k):ind(k)+100),'apa', 'once'))
	useind=[useind;k];
	end
	end
	use=[];
	for k=1:length(useind)
	if ~isempty(regexp(s(ind(useind(k)):ind(useind(k))+200),'\$', 'once'))
	use=[use;ind(useind(k))];
	end
	end
	prices=cell(1,length(use));
	for k=1:length(use)
	prices{k}=str2double(s(use(k)+regexp(s(use(k):use(k)+200),'\$'):...
	use(k)+regexp(s(use(k):use(k)+200),'\$')+3));
	end
	prices=cell2mat(prices);
	urls=cell(1,length(use));
	for k=1:length(use)
	urls{k}=s(use(k)+regexp(s(use(k):use(k)+200),'http')-1:...
	use(k)+regexp(s(use(k):use(k)+200),'.html')+3);
	end
	urls=char(urls');
	[a,b]=sort(prices);
	urls_sortprice=urls(b,:);
	filter_out_community=[];
	for k=1:length(use)
	tmp=urlread(urls_sortprice(k,:));
	if ~isempty(regexpi(tmp,'community '))
	filter_out_community=[filter_out_community;k];
	end
	end
	urls_sortprice_no_community=urls_sortprice;
	urls_sortprice_no_community(filter_out_community,:)=[];
	prices_no_community=a;
	prices_no_community(filter_out_community)=[];
	maps=cell(1,size(urls_sortprice_no_community,1));
	for k=1:size(urls_sortprice_no_community,1)
	tmp=urlread(urls_sortprice_no_community(k,:));
	maps{k}=tmp(regexp(tmp,'http://maps.google'):regexp(tmp,'">google map')-1);
	end
	fid=fopen('myFile.csv','wt');
	for i=1:length(maps)
	fprintf(fid,'%s\n',[urls_sortprice_no_community(i,:),', ',maps{i},', ',...
	num2str(prices_no_community(i))]);
	end
	fclose(fid);

view raw craig1.m hosted with ❤ by GitHub

Save and execute in matlab, or perhaps add this line to your crontab for a daily update!

/usr/local/MATLAB/R2011b/bin/matlab -nodesktop -nosplash -r "craigslist_search;quit;"

Reduce file size/quality of pdf using ghostscript

gs -sDEVICE=pdfwrite -dCompatibilityLevel=1.4 -dPDFSETTINGS=/screen -dNOPAUSE -dQUIET -dBATCH -sOutputFile=output.pdf input.pdf

view raw reduce.sh hosted with ❤ by GitHub

Write and execute a Matlab script without opening Matlab

Here's an example of using cat to write your m-file and then execute from the command line. This particular example will just plot and print 10 random numbers.

	filename=my_script.m;
	cat > $filename << EOF
	figure
	plot(rand(1,10))
	print -dtiff test.tif
	close
	EOF
	chmod +x $filename
	/usr/local/MATLAB/R2011b/bin/matlab -nodesktop -nosplash -r "my_script;quit;"

view raw writem.sh hosted with ❤ by GitHub

SmartGit on Ubuntu 10.04

SmartGit is a graphical user interface for managing and creating git repositories for software version control

First, uninstall any git you already have (or thought you had)

 sudo apt-get remove git  

Then, get the latest git

 sudo aptitude build-dep git-core  
 wget http://git-core.googlecode.com/files/git-1.7.9.1.tar.gz  
 tar xvzf git-1.7.9.1.tar.gz  
 cd git-1.7.9.1  
 ./configure  
 make  
 sudo make install  
 cd ..  
 rm -r git-1.7.9.1 git-1.7.9.1.tar.gz  

Next, you need Java Runtime Environment 7 which you can download from here

and then follow the installation instructions here

Download SmartGit from here, then

 tar xvzf smartgit-generic-2_1_7.tar.gz  
 cd smartgit-2_1_7  
 cd bin  

Then edit the smartgit shell script

 gedit smartgit.sh  

updating SMARTGIT_JAVA_HOME to where JRE is installed (e.g. /usr/lib/jvm/jrel1.7.0_03). Now you're ready to run

 ./smartgit.sh  

As part of the installation you will be prompted for the location of git, which you can find out by issuing in another terminal:

 whereis git  

Ellipse Fitting to 2D points, part 3: Matlab

The last of 3 posts presenting algorithms for the ellipsoid method of Khachiyan for 2D point clouds. This version in Matlab will return: 1) the ellipse radius in x, 2) the ellipse radius in y, 3) the centre in x, 4) the centre in y, and 5) orientation of the ellipse, and coordinates of the ellipse

	% P is a matrix of [x,y] points
	P=[x,y];
	K = convhulln(P);
	K = unique(K(:));
	PK = P(K,:)';
	[d N] = size(PK);
	Q = zeros(d+1,N);
	Q(1:d,:) = PK(1:d,1:N);
	Q(d+1,:) = ones(1,N);
	% initializations
	count = 1;
	err = 1;
	u = (1/N) * ones(N,1); % 1st iteration
	tolerance=.01;
	while err > tolerance,
	X = Q * diag(u) * Q'; % X = \sum_i ( u_i * q_i * q_i') is a (d+1)x(d+1) matrix
	M = diag(Q' * inv(X) * Q); % M the diagonal vector of an NxN matrix
	[maximum j] = max(M);
	step_size = (maximum - d -1)/((d+1)*(maximum-1));
	new_u = (1 - step_size)*u ;
	new_u(j) = new_u(j) + step_size;
	count = count + 1;
	err = norm(new_u - u);
	u = new_u;
	end
	% (x-c)' * A * (x-c) = 1
	% It computes a dxd matrix 'A' and a d dimensional vector 'c' as the center
	% of the ellipse.
	U = diag(u);
	% the A matrix for the ellipse
	A = (1/d) * inv(PK * U * PK' - (PK * u)*(PK*u)' );
	% matrix contains all the information regarding the shape of the ellipsoid
	% center of the ellipse
	c = PK * u;
	[U Q V] = svd(A);
	r1 = 1/sqrt(Q(1,1));
	r2 = 1/sqrt(Q(2,2));
	v = [r1 r2 c(1) c(2) V(1,1)]';
	% get ellipse points
	N = 100; % number of points on ellipse
	dx = 2*pi/N; % spacing
	theta = v(5); %orinetation
	Rot = [ [ cos(theta) sin(theta)]', [-sin(theta) cos(theta)]']; % rotation matrix
	Xe=zeros(1,N); Ye=zeros(1,N); % pre-allocate
	for i = 1:N
	ang = i*dx;
	x = v(1)*cos(ang);
	y = v(2)*sin(ang);
	d1 = Rot*[x y]';
	Xe(i) = d1(1) + v(3);
	Ye(i) = d1(2) + v(4);
	end

view raw ellipsefit.m hosted with ❤ by GitHub

Ellipse Fitting to 2D points, part 2: Python

The second of 3 posts presenting algorithms for the ellipsoid method of Khachiyan for 2D point clouds. This version in python will return: 1) the ellipse radius in x, 2) the ellipse radius in y x is a 2xN array [x,y] points

	def fitellipse( x, **kwargs ):
	x = asarray( x )
	## Parse inputs
	## Default parameters
	kwargs.setdefault( 'constraint', 'bookstein' )
	kwargs.setdefault( 'maxits', 200 )
	kwargs.setdefault( 'tol', 1e-5 )
	if x.shape[1] == 2:
	x = x.T
	centroid = mean(x, 1)
	x = x - centroid.reshape((2,1))
	## Obtain a linear estimate
	if kwargs['constraint'] == 'bookstein':
	## Bookstein constraint : lambda_1^2 + lambda_2^2 = 1
	z, a, b, alpha = fitbookstein(x)
	## Add the centroid back on
	z = z + centroid
	return z
	def fitbookstein(x):
	'''
	function [z, a, b, alpha] = fitbookstein(x)
	FITBOOKSTEIN Linear ellipse fit using bookstein constraint
	lambda_1^2 + lambda_2^2 = 1, where lambda_i are the eigenvalues of A
	'''
	## Convenience variables
	m = x.shape[1]
	x1 = x[0, :].reshape((1,m)).T
	x2 = x[1, :].reshape((1,m)).T
	## Define the coefficient matrix B, such that we solve the system
	## B *[v; w] = 0, with the constraint norm(w) == 1
	B = hstack([ x1, x2, ones((m, 1)), power( x1, 2 ), multiply( sqrt(2) * x1, x2 ), power( x2, 2 ) ])
	## To enforce the constraint, we need to take the QR decomposition
	Q, R = linalg.qr(B)
	## Decompose R into blocks
	R11 = R[0:3, 0:3]
	R12 = R[0:3, 3:6]
	R22 = R[3:6, 3:6]
	## Solve R22 * w = 0 subject to norm(w) == 1
	U, S, V = linalg.svd(R22)
	V = V.T
	w = V[:, 2]
	## Solve for the remaining variables
	v = dot( linalg.solve( -R11, R12 ), w )
	## Fill in the quadratic form
	A = zeros((2,2))
	A.ravel()[0] = w.ravel()[0]
	A.ravel()[1:3] = 1 / sqrt(2) * w.ravel()[1]
	A.ravel()[3] = w.ravel()[2]
	bv = v[0:2]
	c = v[2]
	## Find the parameters
	z, a, b, alpha = conic2parametric(A, bv, c)
	return z, a, b, alpha
	def conic2parametric(A, bv, c):
	'''
	function [z, a, b, alpha] = conic2parametric(A, bv, c)
	'''
	## Diagonalise A - find Q, D such at A = Q' * D * Q
	D, Q = linalg.eig(A)
	Q = Q.T
	## If the determinant < 0, it's not an ellipse
	if prod(D) <= 0:
	raise RuntimeError, 'fitellipse:NotEllipse Linear fit did not produce an ellipse'
	## We have b_h' = 2 * t' * A + b'
	t = -0.5 * linalg.solve(A, bv)
	c_h = dot( dot( t.T, A ), t ) + dot( bv.T, t ) + c
	z = t
	a = sqrt(-c_h / D[0])
	b = sqrt(-c_h / D[1])
	alpha = atan2(Q[0,1], Q[0,0])
	return z, a, b, alpha

view raw ellipsefit.py hosted with ❤ by GitHub

Ellipse Fitting to 2D points, part 1: R

The first of 3 posts presenting algorithms for the ellipsoid method of Khachiyan for 2D point clouds. This version in R will return: 1) the ellipse radius in x, 2) the ellipse radius in y, 3) the centre in x, 4) the centre in y, and 5) orientation of the ellipse, and coordinates of the ellipse

	#P is an array of points [x,y]
	P=matrix(c(x,y),nc=2)
	K=convhulln(Re(P), options="QJ")
	K=sort(unique(c(K)))
	PK = t(P[K,])
	d= dim(PK)[1]
	N= dim(PK)[2]
	Q = matrix(0,d+1,N)#create a matrix, with dim, filled with zeros.
	Q[1:d,] = PK[1:d,1:N]
	Q[d+1,] = matrix(1,N)# a matrix filled with ones.
	# initializations
	count = 1
	err = 1
	u = (1/N)* matrix(1,N) # 1st iteration
	tolerance=.01
	while (err > tolerance){
	X=Q%*% diag(c(u),length(u))%*%t(Q);X
	M = diag(t(Q) %*% solve(X) %*% Q);M #; % M the diagonal vector of an NxN matrix
	maximum = max(M);maximum
	j=which.max(M);j
	step_size = (maximum - d -1)/((d+1)*(maximum-1));step_size
	new_u = (1 - step_size)* u ;new_u
	new_u[j] = new_u[j] + step_size;new_u
	count = count + 1
	err=sqrt(sum((new_u - u)^2));err
	u = new_u;u
	}
	# compute a dxd matrix 'A' and a d dimensional vector 'c' as the center of the ellipse.
	U=diag(u[1],length(u))
	diag(U)=u
	# the A matrix for the ellipse
	A = (1/d) * solve(PK %*% U %*% t(PK) - (PK %*% u) %*% t(PK %*% u))
	center = PK %*% u # % center of the ellipse
	Q = svd(A)$d
	Q=diag(Q,length(Q))
	U = svd(A)$u
	V = svd(A)$v
	r1 = 1/sqrt(Q[1,1]);
	r2 = 1/sqrt(Q[2,2]);
	# v contains the 1) radius in x, 2) radius in y, 3) centre x, 4) centre y, and 5) orientation of the ellipse
	v = matrix(c(r1, r2, center[1], center[2], V[1,1]));
	# get ellipse points
	N = 100;
	dx = 2*pi/N;
	theta = v[5];
	R = matrix(c(cos(theta),sin(theta), -sin(theta), cos(theta)),2)
	Xe=vector(l=N)
	Ye=vector(l=N)
	for(i in 1:N){
	ang = i*dx;ang
	x = v[1]*cos(ang);x
	y = v[2]*sin(ang);y
	d1 = R %*% matrix(c(x,y));
	Xe[i] = d1[1] + v[3];
	Ye[i] = d1[2] + v[4];
	}

view raw ellipsefit.r hosted with ❤ by GitHub

Non-linear least squares fit forced through a point

Sometimes you want a least-squares polynomial fit to data but to constrain it so it passes through a certain point. Well ... Inputs: x, y (data); x0 and y0 are the coordinates of the point you want it to pass through, and n is the degree of polynomial

	function [p,yhat] = lsq_fit_nonlin_force_thru_point(x,y,x0,y0,n);
	x = x(:); %reshape the data into a column vector
	y = y(:);
	% 'C' is the Vandermonde matrix for 'x'
	V(:,n+1) = ones(length(x),1,class(x));
	for j = n:-1:1
	V(:,j) = x.*V(:,j+1);
	end
	C = V;
	% 'd' is the vector of target values, 'y'.
	d = y;
	%%
	% There are no inequality constraints in this case, i.e.,
	A = [];
	b = [];
	%%
	% We use linear equality constraints to force the curve to hit the required point. In
	% this case, 'Aeq' is the Vandermoonde matrix for 'x0'
	Aeq = x0.^(n:-1:0);
	% and 'beq' is the value the curve should take at that point
	beq = y0;
	%%
	p = lsqlin( C, d, A, b, Aeq, beq );
	%%
	% We can then use POLYVAL to evaluate the fitted curve
	yhat = polyval( p, linspace(0,max(x),100) );

view raw nonlinlstsq.m hosted with ❤ by GitHub

Matlab function to generate a homogeneous 3-D Poisson spatial distribution

This is based on Computational Statistics Toolbox by W. L. and A. R. Martinez, which has a similar algorithm in 2D. It uses John d'Errico's inhull function Inputs: N = number of data points which you want to be uniformly distributed over XP, YP, ZP which represent the vertices of the region in 3 dimensions Outputs: x,y,z are the coordinates of generated points

	function [x,y,z] = csbinproc3d(xp, yp, zp, n)
	x = zeros(n,1);
	y = zeros(n,1);
	z = zeros(n,1);
	i = 1;
	minx = min(xp);
	maxx = max(xp);
	miny = min(yp);
	maxy = max(yp);
	minz = min(zp);
	maxz = max(zp);
	cx = maxx-minx;
	cy = maxy - miny;
	cz = maxz - minz;
	while i <= n
	xt = rand(1)*cx + minx;
	yt = rand(1)*cy + miny;
	zt = rand(1)*cz + minz;
	k = inhull([xt,yt,zt],[xp,yp,zp]);
	if k == 1
	% it is in the region
	x(i) = xt;
	y(i) = yt;
	z(i) = zt;
	i = i+1;
	end
	end

view raw csbinproc3d.m hosted with ❤ by GitHub

3D spectral noise in Matlab with specified autocorrelation

Here is a function to generate a 3D matrix of 1/frequency noise Inputs are 1)rows, 2) cols (self explanatory I hope), 3) dep (the 3rd dimension - number of images in the stack) and 4) factor, which is the degree of autocorrelation in the result. The larger this value, the more autocorrelation an the bigger the 'blobs' in 3D Enjoy!

	imfft=fftshift(fftn(randn(cols,rows,dep)));
	mag=abs(imfft);
	phase=imfft./mag;
	[xi,yi,zi]=meshgrid(1:rows,1:cols,1:dep);
	radius=sqrt(xi.^2+yi.^2+zi.^2);
	radius(cols/2+1,rows/2+1,dep/2+1)=1;
	filter=1./(radius.^factor);
	im=real(ifftn(fftshift(filter.*phase)));
	%im=ifftn(fftshift(filter.*phase),'symmetric');
	im=rescale(im,0,1);
	im=im./sum(im(:));
	im=imresize(im(1:2:cols,1:2:rows,1:dep),2);

view raw specnoise.m hosted with ❤ by GitHub

Computational geometry with Matlab and MPT

For computational geometry problems there is a rather splendid tool called the Multi Parametric Toolbox for Matlab Here are a few examples of using this toolbox for a variety of purposes using 3D convex polyhedra represented as polytope arrays. In the following examples, P is always a domain bounded 3D polytope array

	V=volume(P);
	% return a vector of the number of sides in each polytope
	[H,K]=double(P);
	[numsides, Cncols] = cellfun(@size, H);
	% return a vector of polytope radii
	[xyz,Rad] = chebyball(P);
	% find the extrema of each polytope in P
	for ip=1:length(P)
	E(ip)=extreme(P(ip));
	end
	%Slice through a 3-dimensional volume. Here, I'm slicing through 3D polytope array P to return the number of polytopes in each slice (N) and the mean diameter of polytopes in each slice (M). Slices are made at arbitrary values in the space [0,1]
	M=[]; N=[];
	for ii=[0:.01:.2,.25:.05:.8]
	[Pcut, sliced] = slice(P,2,ii);
	[xc,rc] = chebyball(Pcut); M=[M;mean(rc.*2)]; N=[N;length(rc)];
	end %ii
	% here I find the indices where of polytopes whose centres are in the space bounded by x=min and y=max
	index=[];
	for ip=1:length(P)
	E=extreme(P(ip));
	if(any((r(1,ip)+E(:,1))<min) any((r(2,ip)+e(:,2))="" ||=""><min) any((r(3,ip)+e(:,3))="" ||=""><min)) continue="" elseif(any((r(1,ip)+e(:,1))="">max) || any((r(2,ip)+E(:,2))&gt;max) || any((r(3,ip)+E(:,3))&gt;max))
	continue
	else
	index=[index;ip];
	end
	% Here I reduce each 3D polytope in array P by factor fct to give new polytope array Q
	Q = [];
	for ip=1:length(P)
	E=extreme(P(ip));
	ey=E(:,2); ex=E(:,1); ez=E(:,3);
	ex = fct*ex+(1-fct)*mean(ex); % expand x by factor fct
	ey = fct*ey+(1-fct)*mean(ey); % expand y
	ez = fct*ez+(1-fct)*mean(ez); % expand y
	Er=[ex,ey,ez];
	Q = [Q,hull(Er)];
	end
	% Finally, here's a piece of code which will a) convert a 3D polytope array into a 3D matrix, then b) display slices of that volume (left side, back, and the centre line in both the x and y orientation) % first, slice P at small increments and build up matrix im
	im=zeros(501,601,length([0:.01:1]));
	counter=1;
	for i=[0:.01:1]
	Pcut=slice(Q,2,i);
	plot(Pcut)
	axis([0 1 0 1])
	axis off
	print -dtiffnocompression tmp.tiff
	I=double(rgb2gray(imread('tmp.tiff')));
	Ib=(rescale(I,0,1)~=1);
	Ib=Ib(200:700,400:1000);
	im(:,:,counter)=Ib;
	counter=counter+1;
	end
	close all
	delete tmp.tiff
	% resize by 50% (helps speed things up)
	im_50percent=imageresize(im,.5,.5);
	im=im_50percent;
	% create a grid
	[x,y,z]=meshgrid([1:size(im,2)],[1:size(im,1)],[1:size(im,3)]);
	xmin = min(x(:));
	ymin = min(y(:));
	zmin = min(z(:));
	xmax = max(x(:));
	ymax = max(y(:));
	zmax = max(z(:));
	% get the first slice
	hslice = surf(linspace(xmin,xmax,100),...
	linspace(ymin,ymax,100),...
	zeros(100));
	rotate(hslice,[-1,0,0],-90)
	xd = get(hslice,'XData');
	yd = get(hslice,'YData');
	zd = rescale(get(hslice,'ZData'),zmin,zmax);
	delete(hslice)
	% then build more slices on top
	h = slice(x,y,z,im,xd,yd,zd);
	set(h,'FaceColor','flat',...
	'EdgeColor','none',...
	'DiffuseStrength',5)
	colormap gray
	hold on
	hslice = surf(linspace(xmin,xmax,100),...
	linspace(ymin,ymax,100),...
	zeros(100));
	rotate(hslice,[0,-1,0],-90)
	xd = get(hslice,'XData');
	yd = get(hslice,'YData');
	zd = rescale(get(hslice,'ZData'),zmin,zmax);
	delete(hslice)
	h = slice(x,y,z,im,xd,yd,zd);
	set(h,'FaceColor','flat',...
	'EdgeColor','none',...
	'DiffuseStrength',5)
	hx = slice(x,y,z,im,xmax,[],[]);
	set(hx,'FaceColor','flat','EdgeColor','none')
	hy = slice(x,y,z,im,[],ymax,[]);
	set(hy,'FaceColor','flat','EdgeColor','none')
	% hz = slice(x,y,z,im,round(xmax/2),[],[]);
	% set(hz,'FaceColor','flat','EdgeColor','none')
	axis tight
	box on
	camproj perspective
	lightangle(-45,45)
	set(gcf,'Renderer','zbuffer')
	warning on all

view raw mpt_examples.m hosted with ❤ by GitHub

Restoring Nvidia graphics on Ubuntu

If you're an idiot like me and do this:

sudo apt-get --purge remove nvidia-*

You'll soon realise your display looks a little crappy, because you've turned off your graphics card drivers. Not to worry! You restore like this:

sudp apt-get install nvidia-*

then go into /etc/modprobe.d/blacklist.conf and remove any 'blacklists' to nvidia drivers, then restart the x-server using

sudo /etc/init.d/gdm restart

Normal service resumed, and you are free to make a new balls-up!

Wrapping and interfacing compiled code into more flexible and user-friendly programs using Zenity and Bash

Subtitles include:

1. How to reconstruct the three-dimensional shape of an object in a scene in a very inexpensive way using shadow-casting
2. How to write and execute a matlab script from a bash script

I've written a bash script - shadowscan.sh - which is a wrapper for a set of programs for the reconstruction 3D surface of an object in a series of photographs using the ingenious shadow-casting algorithm of Jean-Yves Bouguet and Pietro Perona

The script relies heavily on the GUI-style interfacing tools of Zenity and acts as an interface and wrapper for the algorithm as implemented by these guys. It will also write and execute a matlab script for visualising the outputs.

This script should live in a directory with the following compiled programs, all available here:
1. ccalib (performs the camera calibration to retrieve the intrinsic parameters of camera and desk)
2. find_light (performs the light calibration to find light-source coordinates)
3. desk_scan (does the 3d reconstruction)
4. merge (merges 2 3d reconstructions, e.g. from different angles)

and the following matlab script:
1. selectdata.m (available here)

When the makefile is run, it creates the compiled programs in the ./bin directory. It's important to consult the README in the main folder as well as within the 'desk_scan' folder to have a clear idea of what the program does and what the input parameters mean.

My program will check if the camera calibration and light calibration files exist in the folder and if not will carry them out accordingly then several reconstructions are carried out, and the outputs merged. Finally matlab is called from the command line to run 'p_lightscan_data.m' to visualize the results.

when the script finishes there will be the following output files:
1. xdim.txt - max dimension to scan to
2. params - camera calibration parameters
3. data - camera calibration data
4. light.txt - light calibration parameters
5. sample.pts - [x,y,z] of object, in coordinates relative to light source
6. sample_lightscan1.tif - filled colour contour plot of reconstructed object
7. sample_lightscan2.tif - mesh plot of reconstructed object

I've tried to keep it flexible. It will run only the elements it needs to. For example, if the camera or lighting has been carried out before for a different object it will pick up on that. Also, the program can be invoked by passing it two arguments (path to folder where sample images are, and what sequence of contrasts to use for the reconstruction), which avoids having to input with the graphical menus. Example usage:
1. bash shadowscan.sh
2. bash shadowscan.sh /home/daniel/Desktop/shadow_scanning/sample/bottle 30.35.40.45

	defined()
	{
	[ ${!1-X} == ${!1-Y} ]
	}
	args=("$@") # get arguments passed to the program
	# if no arguments, IMAGES and CNTRST will be empty otherwise will be filled with the arguments
	defined args && IMAGES=${args[0]} # first = sample image location e.g. /home/daniel/Desktop/shadow_scanning/sample/
	defined args && CNTRST=${args[1]} # second = contrst values - no more than four! e.g. 30.35.40.45
	# if IMAGES empty, display flash info, and prompt user to direct program to folder
	if [ -z $IMAGES ] ; then
	zenity --info --title="SHADOWSCAN" --text "This program will create a 3D reconstruction of an object using images of the object scanned by a shadow. It is a bash/matlab wrapper program for code created by http://www.cs.columbia.edu/~pblaer/projects/photo/ which implements the shadow-casting algorithm of Jean-Yves Bouguet and Pietro Perona (http://www.vision.caltech.edu/bouguetj/ICCV98/)"
	# folder where sample images are
	IMAGES=$(zenity --file-selection --directory --title="Select directory where sample images are")
	fi
	if [ ! -f ./param ] # if param does not exist, do cam calib
	then
	# folder where checkerboard images are
	CAL_CAM=$(zenity --file-selection --title="Select checkerboard image" --text "")
	# get x-dimension size (number of pixels), subtract 20, and write to file
	convert $CAL_CAM -format "%[fx: min(w,h)-20]" info: > xdim.txt
	# input size of square in calibration photo # 23
	SQ_HGHT=$(zenity --title="Camera Calibration" --text "input square size" --entry)
	./ccalib $CAL_CAM $SQ_HGHT
	fi
	if [ ! -f ./light.txt ] # if light.txt does not exist, do light calib
	then
	# folder where light calibration images are
	CAL_LIGHT=$(zenity --file-selection --directory --title="Select directory where light calibration images are")
	# input height of the pencil in the light calibration photos # 76
	OBJ_HGHT=$(zenity --title="Light Calibration" --text "input object height" --entry)
	zenity --info --title="Instructions" --text "For each image, click first on the shadow tip, then on the base"
	./find_light param $OBJ_HGHT $CAL_LIGHT/*.pgm > light.txt
	fi
	# this is the reference at the bottom of the image. This
	# should be a row where at no time does the shadow pass
	# over an object, it always should project onto the desk plane.
	if [ ! -f ./xdim.txt ] ; then # if xdim.txt does not exist
	BTM=$(zenity --title="reference value" --text "at the bottom of the image" --entry)
	else
	BTM=$(head xdim.txt)
	fi
	# this is the reference at the top of the image. This
	# should be a row where at no time does the shadow pass
	# over an object, it always should project onto the desk plane.
	TOP=20
	# if CNTRST empty, prompt user to select which CNTRST to use
	if [ -z $CNTRST ] ; then
	# user selects contrast levels (1 to 4) to be processed from a limited list
	CNTRST=$(zenity --list --text "Select one or more contrast values to process" --title="" --checklist --column "Pick" --column "options" TRUE "30" TRUE "35" FALSE "40" FALSE "45" --separator=".")
	fi
	# $Var are contrast values, 0 is desk plane
	if [ ${#CNTRST} -eq 2 ]; then # 1 contrast
	Var1=`echo $CNTRST | awk -F\. '{print $1}'`
	./desk_scan light.txt param $Var1 0 $TOP $BTM sample.pts $IMAGES/*.pgm
	else
	if [ ${#CNTRST} -eq 5 ] ; then # 2 contrasts
	Var1=`echo $CNTRST | awk -F\. '{print $1}'`
	Var2=`echo $CNTRST | awk -F\. '{print $2}'`
	./desk_scan light.txt param $Var1 0 $TOP $BTM sample1.pts $IMAGES/*.pgm
	./desk_scan light.txt param $Var2 0 $TOP $BTM sample2.pts $IMAGES/*.pgm
	./merge sample1.pts sample2.pts sample.pts
	rm sample1.pts sample2.pts
	else
	if [ ${#CNTRST} -eq 8 ] ; then # 3 contrasts
	Var1=`echo $CNTRST | awk -F\. '{print $1}'`
	Var2=`echo $CNTRST | awk -F\. '{print $2}'`
	Var3=`echo $CNTRST | awk -F\. '{print $3}'`
	./desk_scan light.txt param $Var1 0 $TOP $BTM sample1.pts $IMAGES/*.pgm
	./desk_scan light.txt param $Var2 0 $TOP $BTM sample2.pts $IMAGES/*.pgm
	./desk_scan light.txt param $Var3 0 $TOP $BTM sample3.pts $IMAGES/*.pgm
	# merge the 3 data sets
	./merge sample1.pts sample2.pts sample12.pts
	./merge sample1.pts sample3.pts sample13.pts
	./merge sample2.pts sample3.pts sample23.pts
	./merge sample12.pts sample13.pts sample12_13.pts
	./merge sample12_13.pts sample23.pts sample.pts
	rm sample1.pts sample2.pts sample3.pts sample12.pts sample13.pts sample23.pts sample12_13.pts
	else
	if [ ${#CNTRST} -eq 11 ] ; then # 4 contrasts
	Var1=`echo $CNTRST | awk -F\. '{print $1}'`
	Var2=`echo $CNTRST | awk -F\. '{print $2}'`
	Var3=`echo $CNTRST | awk -F\. '{print $3}'`
	Var4=`echo $CNTRST | awk -F\. '{print $4}'`
	./desk_scan light.txt param $Var1 0 $TOP $BTM sample1.pts $IMAGES/*.pgm
	./desk_scan light.txt param $Var2 0 $TOP $BTM sample2.pts $IMAGES/*.pgm
	./desk_scan light.txt param $Var3 0 $TOP $BTM sample3.pts $IMAGES/*.pgm
	./desk_scan light.txt param $Var4 0 $TOP $BTM sample4.pts $IMAGES/*.pgm
	# merge the 4 data sets
	./merge sample1.pts sample2.pts sample12.pts
	./merge sample3.pts sample4.pts sample34.pts
	./merge sample12.pts sample34.pts sample.pts
	rm sample1.pts sample2.pts sample3.pts sample4.pts sample12.pts sample34.pts
	fi
	fi
	fi
	fi
	# once completed, in matlab run 'p_lightscan_data.m' to visualize the results
	# http://ubuntuforums.org/showthread.php?t=825042
	filename=p_lightscan_data.m;
	cat > $filename << EOF
	% p_lightscan_data
	% example matlab code to clean up and plot data generated from shadowscan.sh
	% Daniel Buscombe March 2011
	close all;clear all;clc
	pts=importdata('sample.pts');
	plot(pts.data(:,1),pts.data(:,2),'k.')
	view(2)
	[pl,xs,ys] = selectdata('sel','rect');
	close
	x=pts.data(pl,1); y=pts.data(pl,2); z=pts.data(pl,3);
	f=find(z==0); x(f)=[]; y(f)=[]; z(f)=[];
	[xi,yi]=meshgrid([min(x):.1:max(x)],[min(y):.1:max(y)]);
	zi=griddata(x,y,z,xi,yi);
	figure
	contourf(xi,yi,-zi-min(-zi(:))), colorbar
	print -dtiff -r300 sample_lightscan1.tif
	close
	figure
	mesh(xi,yi,-zi-min(-zi(:)))
	print -dtiff -r300 sample_lightscan2.tif
	close
	EOF
	chmod +x $filename
	/opt/matlab/bin/matlab -nodesktop -nosplash -r "p_lightscan_data;quit;"

view raw shadowscan.sh hosted with ❤ by GitHub

Program for spatial correlogram of an image, in C

This program will read an image, convert it to greyscale, crop it to a square, and calculate it's spatial (1D) autocorrelation function (correlogram), printing it to screen.

It uses OpenCV image processing libraries - see here and here

Compilation instructions (in Ubuntu):

 sudo apt-get install libcv-dev libcv1 libcvaux-dev libcvaux1 libhighgui-dev libhighgui1  
 pkg-config --libs opencv; pkg-config --cflags opencv  
 PKG_CONFIG_PATH=/usr/lib/pkgconfig/:${PKG_CONFIG_PATH}  
 export PKG_CONFIG_PATH  
 gcc `pkg-config --cflags --libs opencv` -o Rspat Rspat.c  

Usage in the terminal:

 ./Rspat IMG_0208.JPG   

The program:

	#include <stdio.h>
	#include "cv.h"
	#include "highgui.h"
	int main( int argc, char** argv )
	{
	time_t t1,t2; // initiate timer
	(void) time(&amp;t1);
	// define some variables
	IplImage *im = 0;
	/* check for supplied argument */
	if ( argc <> \n" );
	return 1;
	}
	// load image
	im = cvLoadImage( argv[1], CV_LOAD_IMAGE_GRAYSCALE );
	CvSize sz = cvGetSize(im); // get image dimensions
	int st=(sz.width/2)-(sz.height/2);
	// sets the Region of Interest
	cvSetImageROI(im, cvRect(st, 0, sz.height, sz.height));
	// create destination image
	IplImage *img1 = cvCreateImage(cvGetSize(im),
	im->depth,
	im->nChannels);
	// copy subimage > crop
	cvCopy(im, img1, NULL);
	// reset ROI
	cvResetImageROI(im);
	int M = img1->width; //integer image dimensions
	int N = img1->height;
	int m; // initiate loop counters
	int n;
	CvScalar s1; // initiate temporary storage of scalars
	CvScalar s2;
	CvScalar img1_avg = cvAvg( img1, NULL ); // find average of img1
	int l;
	// initiate main loop
	for( l=1; l<100 img2 =" cvCreateImage(cvGetSize(img1),">depth,
	img1->nChannels);
	// copy subimage > create shifted
	cvCopy(img1, img2, NULL);
	// reset ROI
	cvResetImageROI(img1);
	CvScalar img2_avg = cvAvg( img2, NULL ); // find average of img2
	// reset counters
	double sum_img1_img2 = 0;
	double sum_img1_2 = 0;
	double sum_img2_2 = 0;
	for( m=0; m<M-l ; ++m ){
	for( n=0; n<N-l; ++n ){
	s1=cvGet2D(img1,m,n); // get the (m,n) pixel value of img1
	s2=cvGet2D(img2,m,n); // get the (m,n) pixel value of img2
	sum_img1_img2 = sum_img1_img2 + (s1.val[0]-img1_avg.val[0])*(s2.val[0]-img2_avg.val[0]);
	sum_img1_2 = sum_img1_2 + (s1.val[0]-img1_avg.val[0])*(s1.val[0]-img1_avg.val[0]);
	sum_img2_2 = sum_img2_2 + (s2.val[0]-img2_avg.val[0])*(s2.val[0]-img2_avg.val[0]);
	}
	}
	double corr = sum_img1_img2/sqrt(sum_img1_2*sum_img2_2); // autocorrelation coefficient, lag l
	printf("lag%i,corr=%f\n",l,corr);
	(void) time(&t2);
	return 0;
	}

view raw Rspat.c hosted with ❤ by GitHub

Controlling a serial device with Python

This program to sends a command to serial device, receives data from the device, converts it to numeric format, and graphs the data in real time. It is written in Python using wxPython for the display.

The device is connected to /dev/ttyS0 with a baud rate of 19200, 8-bit data, no parity, and 1 stopbit. You need to edit the line 'serial.Serial' with the relevant parameters to use this code for your device. Editing 'DATA_LENGTH' will change the rate at which the graph scrolls forward. Also, you'd need to change the 'command' that the device understands as 'send me data' (here it is the s character - line 'input'). Knowledge of the format of the returned data is also essential (here I turn ascii into 2 integers) line 'struct.unpack' - see here for help). Otherwise this should work with any device. Tested on Ubuntu 9.10 and 10.4 using serial port and USB-serial converters.

	import os
	import pprint
	import sys
	import wx
	# use wx with mpl is with the WXAgg backend
	import matplotlib
	matplotlib.use('WXAgg')
	from matplotlib.figure import Figure
	from matplotlib.backends.backend_wxagg import FigureCanvasWxAgg as FigCanvas, NavigationToolbar2WxAgg as NavigationToolbar
	import numpy as npy
	import pylab
	import serial
	import time
	import struct
	# set connection:
	ser = serial.Serial('/dev/ttyUSB0', 19200, bytesize=serial.EIGHTBITS, parity=serial.PARITY_NONE, stopbits=1, timeout=0.1, xonxoff=1, rtscts=0)
	ser.open() # open connection
	ser.isOpen() # should be true
	DATA_LENGTH = 120 #seconds
	REDRAW_TIMER_MS = 60 # milliseconds
	def getData():
	input = 's'
	# send 's' the character to the device
	ser.write(input)
	out = ''
	# give device time to answer
	time.sleep(0.02)
	while ser.inWaiting() > 0:
	out += ser.read(1) #request data (serial.read() only returns 8 bit ASCII strings)
	#out += ser.read(ser.inWaiting())
	if out != '':
	# read in a length-2 string of 8-bit characters and then use the struct module to convert the strings into integers
	response=struct.unpack(">bb", out)
	# use only first integer in tuple (2nd is empty)
	print response[0]
	return response[0]
	class GraphFrame(wx.Frame):
	# the main frame of the application
	def __init__(slf):
	wx.Frame.__init__(slf, None, -1, "Sensor response", size=(800,600))
	slf.Centre()
	slf.data = []
	slf.create_main_panel()
	slf.redraw_timer = wx.Timer(slf)
	slf.Bind(wx.EVT_TIMER, slf.on_redraw_timer, slf.redraw_timer)
	slf.redraw_timer.Start(REDRAW_TIMER_MS)
	def create_main_panel(slf):
	slf.panel = wx.Panel(slf)
	slf.init_plot()
	slf.canvas = FigCanvas(slf.panel, -1, slf.fig)
	slf.hbox1 = wx.BoxSizer(wx.HORIZONTAL)
	slf.vbox = wx.BoxSizer(wx.VERTICAL)
	slf.vbox.Add(slf.canvas, 1, flag=wx.LEFT | wx.TOP | wx.GROW)
	slf.panel.SetSizer(slf.vbox)
	def init_plot(slf):
	slf.dpi = 100
	slf.fig = Figure((3.0, 3.0), dpi=slf.dpi)
	slf.axes = slf.fig.add_subplot(111)
	slf.axes.set_axis_bgcolor('black')
	slf.axes.set_title('My data', size=12)
	slf.ymin = -50000 # initalise with large values
	slf.ymax = 50000
	pylab.setp(slf.axes.get_xticklabels(), fontsize=8)
	pylab.setp(slf.axes.get_yticklabels(), fontsize=8)
	# plot the data as a line series, and save the reference
	# to the plotted line series
	slf.plot_data = slf.axes.plot(
	slf.data,
	linewidth=2,
	color="white",
	)[0]
	#count=count+1
	#if count>120:
	# slf.Destroy()
	def draw_plot(slf):
	# redraws the plot
	xmax = len(slf.data) if len(slf.data) > DATA_LENGTH else DATA_LENGTH
	xmin = xmax - DATA_LENGTH
	ymin = slf.ymin # ymin and ymax set based on previous entries
	ymax = slf.ymax #50000
	slf.axes.set_xbound(lower=xmin, upper=xmax)
	slf.axes.set_ybound(lower=ymin, upper=ymax)
	pylab.setp(slf.axes.get_xticklabels(), visible=True)
	slf.plot_data.set_xdata(npy.arange(len(slf.data)))
	slf.plot_data.set_ydata(npy.array(slf.data))
	slf.canvas.draw()
	def on_redraw_timer(slf, event):
	newData = getData()
	slf.data.append(newData)
	slf.draw_plot()
	data=slf.data
	if len(data)>=6:
	# get last 5 values
	datanew=[data[len(slf.data)-5],data[len(slf.data)-4],data[len(slf.data)-3],data[len(slf.data)-2],data[len(slf.data)-1]]
	slf.ymin = (min(datanew))-(min(datanew)*10)/100 # min of last 5, minus 10%
	slf.ymax = (max(datanew))+(max(datanew)*10)/100 # max of last 5, plus 10%
	else:
	slf.ymin = (min(data))-(min(data)*10)/100 # min minus 10%
	slf.ymax = (max(data))+(max(data)*10)/100 # max plus 10%
	def on_exit(slf, event):
	slf.Destroy()
	def flash_status_message(slf, msg, flash_len_ms=1500):
	slf.statusbar.SetStatusText(msg)
	slf.timeroff = wx.Timer(slf)
	slf.Bind(
	wx.EVT_TIMER,
	slf.on_flash_status_off,
	slf.timeroff)
	slf.timeroff.Start(flash_len_ms, oneShot=True)
	def on_flash_status_off(slf, event):
	slf.statusbar.SetStatusText('')
	if __name__ == '__main__':
	app = wx.PySimpleApp()
	app.frame = GraphFrame()
	app.frame.Show()
	app.MainLoop()

view raw serialcontrol.py hosted with ❤ by GitHub

Using cron to backup Thunderbird emails

My bash script (which I've saved as backup_emails.sh in my home folder and made executable) looks a little like this:

 set -x   
  x=`date +%y%m%d.%H%M`   
  tar zcf thunderb-mail-${x}.tgz ~/.mozilla-thunderbird   
  mv thunderb-mail-${x}.tgz /media/MY_HARDDRIVE/email_backups/   

Then invoke in cron using:

 crontab -e   

add something like:

 30 4 * * * ~/./backup_emails.sh  

which does the backup at 04.30 every day

Make Google Earth Faster on Ubuntu

In terminal type:

 export LIBGL_ALWAYS_INDIRECT=true  

Batch GIMP script for auto-sharpen, white-balance and colour enhance

Like you, my point-and-shoot camera in auto-ISO mode often, if not always, gives me pictures which need a little adjustment to realise their full potential.

In GIMP, I usually use the 'auto white balance' option which usually gives great results. Often, I also use the 'auto color enhance' option and, occasionally, I use a simple filter to sharpen the features in the image.

Here's a GIMP script which will automate the process for a folder of pictures

	(define (batch-auto-fix pattern
	radius
	amount
	threshold)
	(let* ((filelist (cadr (file-glob pattern 1))))
	(while (not (null? filelist))
	(let* ((filename (car filelist))
	(image (car (gimp-file-load RUN-NONINTERACTIVE
	filename filename)))
	(drawable (car (gimp-image-get-active-layer image))))
	(plug-in-unsharp-mask RUN-NONINTERACTIVE
	image drawable radius amount threshold)
	(gimp-levels-stretch drawable)
	(plug-in-color-enhance RUN-NONINTERACTIVE
	image drawable)
	(gimp-file-save RUN-NONINTERACTIVE
	image drawable filename filename)
	(gimp-image-delete image))
	(set! filelist (cdr filelist)))))

view raw fix hosted with ❤ by GitHub

There are four input parameters. The first is the 'pattern' to find the files, e.g. "*.tiff" gives you all the tiff files in the directory. The next three inputs refer only to the sharpening filter. I use the (confusingly named) 'unsharp' filter which works by comparing using the difference of the image and a gaussian-blurred version of the image.

'radius' (pixels) of the blur (>1, suggested 5)
'amount' refers to the strength of the effect (suggested 0.5)
'threshold' refers to the pixel value beyond which the effects are applied (0 - 255, suggested 0)

These require a little bit of trial and error, depending on the type of picture. In the command line, the script is called like this:

 gimp -i -b '(batch-auto-fix "*.JPG" 5.0 0.5 0)' -b '(gimp-quit 0)'  

An obligatory example. Before:




and after:






Enjoy!

Convert movies to be compatible with ZEN X-Fi media player

Creative - another company which doesn't provide support for non-Windows users. Get with the times, losers!

The ZEN player needs a specific format of movie (bit-rate and aspect ratio). Here's how to do it with mencoder (available for free on all platforms)

 mencoder 'MyMovie.avi' -oac mp3lame -lameopts cbr:mode=2:br=96 -af resample=44100 -srate 44100 -ofps 20 -ovc lavc -lavcopts vcodec=mpeg4:mbd=2:cbp:trell:vbitrate=300 -vf scale=320:240 -ffourcc XVID -o MyMovie4zen.avi  

Granular Therapy

Here's a relaxing movie I just made of spinning modelled sand grains:



I'm not going to divulge the details of the model of how to create the sand grains (a current research project of mine). I can tell you I used these rather excellent functions for MATLAB:

Volume interpolation Crust algorithm, described more here

A toolbox for Multi-Parametric Optimisation

I used 'recordmydesktop' to capture the grains in my MATLAB figure window. It created an OGV format video (which I renamed to 'beautiful_grains.ogv'). I then used the following code to decode the video into jpegs

 mkdir out_pics  
 mplayer -ao null -ss 00:00:03 -endpos 50 beautiful_grains.ogv -vo jpeg:outdir=out_pics  

I needed to crop the pictures to just the grains, then turn the cropped pictures back into a video. I used this:

 cd out_pics  
 for image in *.jpg  
 do  
 # crop the image  
 convert $image -crop 600x600+350+135 grains%05d.tiff  
 done  
 # make movie from the cropped tiffs  
 ffmpeg -r 25 -i grains%05d.tiff -y -an spinning_grains_small.avi  

et voila! Now sit back with a cuppa and relax to those spinning grains.

Updated Power Spectrum from an Image

A while back I posted on using ImageMagick to do a full frequency-decomposition on an image, and create a power spectrum. It had several inaccuracies. Here's an update with how to do it better

First, you need to install a version of IM which is HDRI enabled (see here)

# test for valid version of Image Magick

 im_version=`convert -list configure | \  
  sed '/^LIB_VERSION_NUMBER /!d; s//,/; s/,/,0/g; s/,0*\([0-9][0-9]\)/\1/g'`  
 [ "$im_version" -lt "06050407" ] && errMsg "--- REQUIRES IM VERSION 6.5.4-7 OR HIGHER ---"  

# test for hdri enabled

 hdri_on=`convert -list configure | grep "enable-hdri"`  
 [ "$hdri_on" = "" ] && errMsg "--- REQUIRES HDRI ENABLED IN IM COMPILE ---"  

#read $file, load and converting to greyscale

 convert $file -colorspace gray in_grey.tiff  

# crop image to smallest dimension

 min=`convert in_grey.tiff -format "%[fx: min(w,h)]" info:`  

# crop image, and convert to png (why?)

 convert in_grey.tiff -background white -gravity center -extent ${min}x${min} in_grey_crop.png  

# get data, find mean as proportion, convert mean to percent, and de-mean image

 data=`convert in_grey_crop.png -verbose info:`  
 mean=`echo "$data" | sed -n '/^.*[Mm]ean:.*[(]\([0-9.]*\).*$/{ s//\1/; p; q; }'`  
 m=`echo "scale=1; $mean *100 / 1" | bc`  
 convert in_grey_crop.png -evaluate subtract ${m}% in_grey_crop_dm.png  

# pad to power of 2

 p2=`convert in_grey_crop_dm.png -format "%[fx:2^(ceil(log(max(w,h))/log(2)))]" info:`  
 convert in_grey_crop_dm.png -background white -gravity center -extent ${p2}x${p2} in_grey_crop_dm_pad.png  

# fft

 convert in_grey_crop_dm_pad.png -fft in_ft.png  

# get spectrum

 convert in_ft-0.png -contrast-stretch 0 -evaluate log 10000 in_spec.png  

... but it's still not quite right. Any tips?

Turn your OGV movie section into a cropped animated GIF

An example from a recent project

 # make directory for output jpg  
 mkdir out_pics   
 # section to the ogv video file into jpegs  
 mplayer -ao null -ss 00:00:03 -endpos 50 beautiful_movie.ogv -vo jpeg:outdir=out_pics  
 cd out_pics  
 # loop over each image and crop into desired rectangle  
 for image in *.jpg  
 do  
 convert $image -crop 600x600+350+135 "$image crop.tiff"  
 done  
 # turn into an animated gif  
 convert -page 600x600+0+0 -coalesce -loop 0 *.tiff beautiful_animation.gif  

Miscellaneous useful bash

For my benefit more than anything - things I use over and over again, but do not have the mental capacity to remember

# number of files in current directory

 ls -1R | grep .*.<file extension> | wc -l  

# check if you have a program installed

 aptitude show <name of program> | grep State  

# only list the (sub-)directories in your current directory

 ls -l | grep ^d  

# remove backup files

 find ./ -name '*~' | xargs rm  

# use zenity get the directory where the file you want is

 zenity --info --title="My Program" --text="Select  
 which directory contains the file(s) you want"  
 export target_dir=`zenity --file-selection --directory`  

# list all txt files in this directory

 x=$( ls $target_dir | grep .*.txt |sort|zenity --list --title="Select" --column="coordinate file" )  

# get this directory

 export raw_dir=`zenity --file-selection --directory`  

# removes .txt (for example, for later use as start of filename)

 str="$x"  
 Var1=`echo $str | awk -F\. '{print $1}'`  

# gets the rest of the file name

 str="$Var1"  
 Var2=`echo $str | awk -F\- '{print $2}'`   

# turn string into (indexable) array

 e=$(echo $mystring | sed 's/\(.\)/\1 /g')  
 b=($e)  

# get max and min dimensions of an image

 min=`convert image1.tiff -format "%[fx: min(w,h)]" info:`  
 max=`convert image1.tiff -format "%[fx: max(w,h)]" info:`  

# take away 1 from both

 min=$((min-1))   
 max=$((max-1))