==============================================================================
function [d c r] = t71(S)
% S is a Nx2 or Nx3 matrix representing N points
N = size(S,1);
% find the diameter of S - the greatest distance between 2 points
d=0;
for i=1:N
for j=(i+1):N
newd = norm( S(i,:) - S(j,:) );
if newd > d
d = newd;
end
end
end
% find the center of S - the average of the points
c = zeros(1, size(S,2));
for i=1:N
c = c + S(i,:);
end
c = c/N;
% find the radius of S - the maximal distance from c to one of the points
r=0;
for i=1:N
newr = norm( S(i,:) - c );
if newr > r
r = newr;
end
end
==============================================================================
Sample runs:
============================================================================== >> A=rand(100,2); >> [d c r]=t71(a) d = 1.2627 c = 0.5290 0.5409 r = 0.6346 >> A=rand(100,3); >> [d c r]=t71(A) d = 1.4676 c = 0.5053 0.4343 0.5017 r = 0.7574 ==============================================================================Things to note here: in the unit square(cube) the diameter must be less than sqrt(2)=1.4142 (sqrt(3)=1.7321). The diameter is always less than or equal to twice the radius.
==============================================================================
function [ c r ] = t72(S)
% find the minimum area circle that bounds the points in the matrix S
% S is a N x 2 or N x 3 matrix
% c will be the center, r will be the radius
% the function t72b(x,S) computes the minimum radius of a circle,
% center x, that bounds the points in S. we want to minimize
% this, as a first guess use the "average" of the points
N=size(S,1);
c=0;
for i=1:N
c = c + S(i,:) ;
end
c=c/N;
[c r] = fminsearch( @(x) t72b(x,S) , c );
%let's make a plot of the center and the points in the 2d case
if size(S,2) == 2
figure(1)
hold on
plot( S(:,1), S(:,2), 'b.') % the points of S in blue
plot( c(1) , c(2) , 'r.') % the center in red
t=0:0.01:6.3;
x=c(1)+r*cos(t);
y=c(2)+r*sin(t);
plot(x,y,'g') % the circle in green
end
==============================================================================
function r = t72b(x,S)
N=size(S,1);
r=0;
for i=1:N
newr=norm( S(i,:) - x );
if newr > r
r = newr;
end
end
==============================================================================
Sample runs:
============================================================================= >> A=rand(100,3); >> [c r]=t72(A) c = 0.5167 0.4595 0.4710 r = 0.7275 >> A=rand(100,2); >> [c r]=t72(A) c = 0.4849 0.4942 r = 0.6668 ==============================================================================In the 2d case it also gives a plot:
==============================================================================
function [ min max av ] = t73(A,d,S)
% given a plane A(1) x + A(2) y + A(3) z = d
% and a Nx3 matrix S representing N points
% find the smallest, biggest and average distance of the points from the plane
N = size(S,1);
min = Inf;
max = 0;
av = 0;
q = sqrt(dot(A,A));
for i=1:N
newd = abs( dot(A,S(i,:)) - d )/q;
if newd < min
min = newd;
end
if newd > max
max = newd;
end
av = av + newd;
end
av = av/N;
==============================================================================
Sample runs:
============================================================================= >> S=rand(50,3); % 50 random points in the unit cube >> A=[1,0,0]; >> d=1; % the plane x=1 >> [min max av] = t73(A,d,S) min = 0.0436 max = 0.9440 av = 0.4706 % all sensible results >> A=[0,1,0]; >> d=3; % the plane y=3 >> [min max av] = t73(A,d,S) min = 2.0240 max = 2.9794 av = 2.5429 % all sensible results >> A=[1,1,1]; >> d=1; >> [min max av] = t73(A,d,S) % the plane x+y+z=1 min = 0.0287 max = 0.9218 av = 0.3573 ==============================================================================
Back to course homepage
Back to my teaching homepage
Back to my homepage