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)