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