Numerics

Numerics provides approximate answers to real world problems.

Toolpath Design and Milling Simulation
Tracking on Homogeneous Manifolds
Sphere-fit to a set of points
Triangulations of elemental shapes in Matlab
Solving Soduko in Matlab
There never will be any good, general methods for
solving systems of more than one nonlinear equation.
Press, Teukolsky, Vetterling, Flannery

Tracking on Homogeneous Manifolds

We present a computation to yield the transformation that matches best two given sets of landmarks on a linearized homogeneous manifold. The method allows to restrict the set of feasible transformations in a way that is most relevant in practical applications. If the homogeneous manifold is a vector space, the transformation is optimal in the limit.

Tracking of a differential-drive robot

Input: transitions of correlated landmarks perceived by a differential-drive robot after a short lapse of time

Output: linear transformation that is consistent with the transitions and the mobility of the robot: driving forward and turning, but no slipping to the side

Time is nature's way of keeping
everything from happening at once.
Allen Stewart Königsberg
Tracking on Homogeneous Manifolds tracking_on_hom... 180 kB
Examples and illustrations (Mathematica 5.2) * tracking.nb 35 kB
* includes an implementation of the Campbell-Baker-Hausdorff formula

Our computation is an alternative to the conventional algorithms in tracking (that might utilize quaternions, or Euler-angles). Our method minimizes the sum of squared errors between the pairs of landmarks in the match. The number of variables is precisely the dimensions of the mobility of the robot. To solve the non-linear equations, Newtons iteration works well.

Tracking of a robot in 3d

Input: correlated landmarks in 3d

Output: in this example, the robot is free to rotate, but can translate only along the first axis

In the paper, Meshless Deformations Based on Shape Matching by M. Mueller, B. Heidelberger, M. Teschner, M. Gross an optimal rotation is found via the polar decomposition of a matrix. This motivated me to make a tracking algorithm using the Campbell-Baker-Hausdorff series in the context of homogeneous manifolds.

If everything seems under control,
you're just not going fast enough.
Mario Andretti

Sphere-fit to a set of points


Example for n=2 and m=5

Given a set of points, we find the center and radius of a sphere that fits the points well. The method is extracted from code of Izhak Bucher.

In dimensional Euclidean space, points on a sphere with center and radius are characterized by

Equivalently,

where the summations are over the coordinates . Now, we are given a set of points

To obtain the center of the fitting sphere , we solve the following system of linear equations in the least square sense

Then, the radius of the sphere is

If the matrix above does not have maximum rank, we define .


Example for n=2 and m=3

As desired, the radius of the sphere is invariant under affine transformation of the points, while the center of the sphere undergoes the same affine transformation. The method is implemented by the following function in Matlab, download. 

function [c,r]=spherefit(X,draw)                % code by jph
%syntax:   [C,R]=SPHEREFIT(X,DRAW)
% X is of size [m n] and contains m points of dimension n.
% C is the center, and R is the radius of the fitting sphere.
% If DRAW is specified and n==2, the result is plotted.
%example:  [c,r]=spherefit(randn(7,2),1)

n=size(X,2);
A=[2*X ones(size(X,1),1)];

if rank(A)<n+1
  c=[];
  r=Inf;
  return;
end

a=A\sum(X.*X,2);
c=a(1:end-1);
r=sqrt(a(end)+c'*c);

if nargin>1&&n==2
  plot(X(:,1),X(:,2),'r*'); 
  hold on; axis equal;
  plot(c(1),c(2),'bx');
  t=linspace(0,2*pi);
  plot(r*cos(t)+c(1),r*sin(t)+c(2),'b');  
  hold off;
end

The following example corresponds to the two dimensional scenario.

To infinity and beyond.
Buzz Lightyear

Triangulations of elemental shapes in Matlab

We release Matlab code to create triangulations of four basic shapes in three dimensions: cube, cylinder, sphere, torus. The functions are 

[T,X,Y,Z]=tricube(n);
[T,X,Y,Z]=tricylinder(n);
[T,X,Y,Z]=trisphere(n);
[T,X,Y,Z]=tritorus(R/r,n);

The parameter n controls the granularity. The generation of the torus requires an additional parameter: the ratio of the big and small radius R/r. To visualize the outputs simply use the Matlab command 

trimesh(T,X,Y,Z);
Triangulations of elemental shapes in Matlab * trishapes.zip 50 kB
* for a demo run tridemo.m in Matlab. The triangulation of a sphere only works on Windows.
Nullius in verba

Solving Soduko in Matlab

Sudoku is a solitaire game involving combinatorics. For the rules visit sudoku.com. According to some mathematician and his computer enumeration there are 9!×722×27×27,704,267,971 different sudoku puzzles - the last factor being prime.

We provide a simple Matlab function sudoku.m to find all solutions to a given sudoku puzzle. Furthermore, the code decides if the puzzle is fiendish, i.e. if branching is necessary.

The code can be naively modified to generate your own puzzles: From an initially full array, successivly delete an entry while the puzzle implies yet uniquely the initially array. Good puzzles are those with few entries given. The current record is 17, i.e. there is a puzzle with 17 initial entries that yet has a unique solution.

Here is the listing of sudoku.m: 

function SOLUTIONS=sudoku(PUZZLE)
%Retrieves all solutions to the sudoku given by PUZZLE.
% PUZZLE    is a 9x9 matrix with 0's as the missing entries
% SOLUTIONS is a 1xN cell where SOLUTIONS{i} is the i-th solution.
%
%Example: P=[0 4 0 0 0 5 1 0 0
%            5 6 0 0 0 3 0 0 2
%            0 0 2 6 0 8 0 0 4
%            7 0 5 0 0 0 0 3 0
%            0 0 0 0 0 0 0 0 0
%            0 1 0 0 0 0 5 0 7
%            2 0 0 9 0 6 4 0 0
%            3 0 0 8 0 0 0 2 6
%            0 0 6 4 0 0 0 8 0];
%         S=sudoku(P);
%         then S is a 1x1 cell and
%         S{1}=[8 4 7 2 9 5 1 6 3
%               5 6 1 7 4 3 8 9 2
%               9 3 2 6 1 8 7 5 4
%               7 8 5 1 6 4 2 3 9
%               4 2 3 5 7 9 6 1 8
%               6 1 9 3 8 2 5 4 7
%               2 5 8 9 3 6 4 7 1
%               3 7 4 8 5 1 9 2 6
%               1 9 6 4 2 7 3 8 5];

global SOLUTIONS
SOLUTIONS=cell(0);

jisudoku(PUZZLE)
function jisudoku(PUZZLE)
global SOLUTIONS

A=cell(9); num=zeros(9); %---initialize---
[X,Y]=meshgrid([0 0 0 3 3 3 6 6 6],[0 0 0 1 1 1 2 2 2]);
T=X+Y+1;
I=repmat(find(T==1),[1 9])+repmat([0 3 6 27 30 33 54 57 60],[9 1]);

%---list available elements for free entries,
%   i.e. unique in row/column/field
for c1=1:numel(PUZZLE)
  x=ceil(c1/9);
  y=mod(c1-1,9)+1;
  if ~PUZZLE(y,x)
    A{c1}=1:9;
    A{c1}=setdiff(A{c1},PUZZLE(:,x));
    A{c1}=setdiff(A{c1},PUZZLE(y,:));
    A{c1}=setdiff(A{c1},PUZZLE(I(:,T(y,x))));
  end
  num(c1)=numel(A{c1});
end
  
if any(any(~num&~PUZZLE))%---free entry but no available => quit!  
  return
end
  
[y,x]=find(num==1);      %---exactly one available => put!
if length(x)
  for c1=1:length(x)
    PUZZLE(y(c1),x(c1))=A{y(c1),x(c1)};
  end
  jisudoku(PUZZLE);
  return
end

cnt=0;
for c1=1:9               %---forall 3x3 fields
  u=[];
  for c2=I(:,c1)'        %----collect availables 
    u=[u A{c2}];
  end
  [u,U]=jiunique(u);  
  if length(intersect(1:9,[u PUZZLE(I(:,c1))']))<9
    %---union in field of available & occupied != 1:9 => quit!
    return
  end
  for c3=find(U==1)      %---entry with unique available => put!
    for c2=I(:,c1)'
      if ismember(u(c3),A{c2})
        PUZZLE(c2)=u(c3);
        cnt=cnt+1;        
      end
    end
  end
end

if cnt
  jisudoku(PUZZLE);
  return
end

if all(all(PUZZLE))      %---found solution
  disp('solution')
  SOLUTIONS{numel(SOLUTIONS)+1}=PUZZLE;
else  
  [mv,mp]=min(reshape(num+10*~~PUZZLE,[1 81]));
  disp('fiendish')
  for c1=A{mp}
    PUZZLE(mp)=c1;
    jisudoku(PUZZLE);
    PUZZLE(mp)=0;
  end
end

function [u,U,I,J]=jiunique(x)
  [u,I,J]=unique(x);
  o=ones(size(x));
  U=full(sparse(o,J,o));
Ich habe so lange ein Motivationsproblem,
bis ich ein Zeitproblem habe.