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