How to find point closest to multiple circles?

Hi!

I want to solve the following problem using Matlab:

In the xy-plane, I have multiple circles, defined by their midpoint and their radius. Some of the circles overlap, others dont.

Now, I want to find a point A in the xy-plane that is closest to all circles.

So in other words: I want to find the point A for which the sum of distances to a point (it does not matter which) on each of the circles is minimal.

Does anyone know how to do this?

Edit:

To clarify my question: the circles can also be completely inside another circle. I added two images of two example situations. In both situations, I want to find point A.

Thanks!

Bas

[Jan]

What exactly is the distance to a set of circles? For one circle it is the set of points of the circle? Or does the distance to the center matter? For 2 circles it matters if they touch each other, overlap or have no contact. Then you get 1, 2 or 1 point, correct?

What is the solution, if all circles have the same center? Then the solution is the circle with the mean radius? This means, that you cannot be sure to get a unique point as solution.

If you have 3 circles, you get 3 distances for a certian point. What is the "closest" point then? The minimal geometric or arithmetic sum of distances? This sounds like the typical problem of a multi-objective minimization.

[Bas Wagemaker]

Hi!

Indeed, those are some interesting questions.

Today I spent a lot of time in investigating the intersection points of the circles. However, since some circles can indeed be inside another circle, a algorithm based on the intersection points did not seem like a very robust solution.

In my application, the circles will never be concentric or identical, which, in this case, seems a good thing. :)

[Image Analyst]

I have the same question as Jan. And your latest comment did not answer his or my questions.

What does distance of your point to a circle mean?

1. Distance to the center?
2. Or the perimeter points?

If it's the center, maybe just averaging the centers will be where you want it. If it's to the perimeter points, then be aware that the radii will affect the answer since there will be more perimeter points on a larger circle, so their distances will give a greater sum than for a small circle.

[Bas Wagemaker]

Hi!

I defined the distance that we are talking about as the distance between the point and one of the perimeter points of each circle. The distance to the center would be easier, since this is one single point per circle, instead of a set.

Since all of you have helped me a lot, I think you deserve some background information. The question is related to my Mechanical Engineering Masters thesis. In fact, the centers of the circles represent points at which one end of a mechanical spring is attached, for different positions of a lever arm. The length of the spring is determined by the required energy in the spring per position. The radius of each circle is equal to the required spring length in each position.

Ideally, all circles intersect in one point, which would be the optimal attachment point of the other end spring (this is why I defined the distance to be the distance between the point and the perimeter). However, in most cases, this ideal point does not exist, and I was therefore looking for the point that is closest to the non-existing intersection point.

Thanks again for your help!

Answer 1

[Matt J]

Let C be a 1x2xN matrix of the circle centers and R an 1x1xN vector of their radii. It is a low-dimensional, non-smooth minimization problem so, I would use fminsearch:

rownorm=@(z) sqrt(  sum(z.^2,2) );
fun=@(A) sum(  rownorm( (A-C) -  (A-C).*R./rownorm(A-C)) , 3 ); %edited
A = fminsearch(fun,A0(:).');

Luckily, it is only a 2D problem, so you can come up with a decent initial guess A0 using a brute force grid search, e.g.,

[X,Y]=ndgrid(linspace(-1,1,100)); %adjust grid limits depending on C and R
Agrid=[X(:),Y(:)];
[~,gridmin]=min( fun(Agrid) );
A0=Agrid(gridmin,:);

[Bas Wagemaker]

Thanks for your reply!

I tried your code (although I do not really understand it). The result:

The red dot represents point A but is not at the location I expect it to be (I expected it to be around (-0.07, 0.14).

I did a initial guess based on a physics formula, which should be (and is) quite close to where I expect point A to be.

Do you have any idea what causes this unexpected behaviour? Thanks!

This is my code:

clc
clear all
close all

%% Variables
mA=5;
rA=0.2;
L=0.5;
kA=1340;
L0spring=0.1;
g=9.81;
Phi=(0:pi/5:pi);

%% Coordinates
xpivot=0; % x coordinate of pivot
ypivot=0; % y coordinate of pivot 
xPivot=ones(size(Phi))*xpivot; % x coordinate of pivot (vector  for plotting)
yPivot=ones(size(Phi))*ypivot; % y coordinate of pivot (vector  for plotting)
xR=rA*sin(Phi); % x coordinate of point R (for all Phi's)
yR=rA*cos(Phi); % y coordinate of point R (for all Phi's)
xMass=L*sin(Phi); % x coordinate of mass (for all Phi's)
yMass=L*cos(Phi); % y coordinate of mass (for all Phi's)

%% Energy
% The energy balance is equal to: Umass+Uspring=CB
% The minimum value of CB is determined by the maximum value of Umass,
% because Uspring=CB-Umass and Uspring>0, thus: CB-Umass>0, thus: CB>Umass
UMass=mA*g*yMass; % Potential energy of mass (for all Phi's)
CBmin=max(UMass); % Minimum energy balance constant

% CB=(CBmin:CBmin/10:CBmin*20);
CB=CBmin;
for j=1:1:length(CB)
    cB=CB(j);
    % For every balancer position, a circle with midpoint equal to point R and
    % radius equal to the spring length is calculated
    % For each balancer position, the x and y coordinates of the points that
    % form the circle are calculated. In the columns of xCircle and yCircle, 
    % the x or y coordinates of the points that form thecircle 
    % corresponding to one balance position are stored.
    [xCircle,yCircle,LSpring]=circleCalculator(UMass,xR,yR,cB,kA,L0spring);
end

C(1,1,:)=xR;
C(1,2,:)=yR;
R(1,1,:)=LSpring;

xA=0;
yA=(mA*g*L)/kA; 
A0=[xA yA]; 

rownorm=@(z) sqrt(  sum(z.^2,2) );
fun=@(A) sum(  rownorm( A -  (A-C).*R./rownorm(A-C)) , 3 );
A = fminsearch(fun,A0(:).');

[Jan]

By the way, see the modern function vecnorm.

You assume, that the "distance to multiple circles" is defined as sum of the distances. Is this the only possible interpretation?

[Matt J]

Your code is incomplete. You have not provided circleCalculator, so we cannot run it. Alternatively, you could just attach a .mat file containing C and R.

[Bas Wagemaker]

Im sorry! Here is the code of circleCalculator:

function [xCircle, yCircle, LSpring] = circleCalculator(UMass,xR,yR,cB,kA,L0spring)
%UNTITLED2 Summary of this function goes here
%   Detailed explanation goes here

USpring=cB-UMass; % Potential energy in spring

%% Spring length
% Since USpring=1/2*k*dL^2, in which dL is the spring elongation, we can
% rewrite this to: dL=sqrt(2*USpring/k)
dLspring=sqrt(2*USpring/kA); % Elongation of spring for every position of the balancer
LSpring=L0spring+dLspring; % Spring length = (null-length)+elongation

%% Circles
% For every balancer position, a circle with midpoint equal to point R and
% radius equal to the spring length is calculated
alpha = (0:pi/50:2*pi)'; % angle to calculate full circle, 0 to 2pi
% For each balancer position, the x and y coordinates of the points that
% form the circle are calculated. In each column, the x or y coordinates of
% the points that form thecircle corresponding to one balance position 
% are stored
xCircle=LSpring.*cos(alpha)+xR; % x coordinates of the points that form the circles
yCircle=LSpring.*sin(alpha)+yR; % y coordinates of the poitns that form the circles

end

[Bas Wagemaker]

Here are the .mat files, thanks a lot!

[Matt J]

The red dot has smaller total distance than the point you say you expect.

K>> fun([-0.07, 0.14])  

ans =

    1.2066

K>> fun([-0.2862    0.1079])  %red dot

ans =

    0.7644

So, either your expectations are wrong, or fun does not represent your true distance metric.

[Bas Wagemaker]

Thanks!

Since I do not really understand your code (yet), I do not know whether fun represents the true distance metric.. But, visually, the red dot does not seem to be the closest point to all circles. Because, at the point where I expected the point A to be, all circles almost intersect. And at an intersection, the distance between two circles is zero, right?

[Matt J]

OK.I had a small bug. With this change to fun, it works now.

fun=@(A) sum(  rownorm( (A-C) -  (A-C).*R./rownorm(A-C)) , 3 );

[Bas Wagemaker]

Thanks, it works like a charm!

Could you please explain what the fun function does?

[Matt J]

You're welcome. Please Accept-click the answer, though, to certify that it addresses the problem.

fun is the thing you are trying to minimize. It computes the sum of the distances from A to your circles.

Answer 2

[Jim Riggs]

Personally, I would use a numerical searcing algorithm.

Overlay a grid on the problem space, then compute the distance function (sum of distance to all circles) from each grid point.

Select the minimum value, then refine the grid and search from the minimum point.

Continue until the desired numerical precision is schieved.

Distance function is the sum of distance from Point P, to each circle. If P is outside the circle, then the distance to that circle is the distance from P to the circle center, minus the circle radius, but if P is inside the circle (Distance from P to the center is less than the radius) then the distance to the circle is the radius minus the distance from P to the center.

[Jan]

The grid search will fail for a set of concentric circles, because you will get a circle as solution, not a single point.

[Bas Wagemaker]

Hi!

Thanks for your reply.

Indeed, while I was writing this question, I thought about the same solution. However, imagine that I define each circle by 100 points. And I have a grid of 10x10 points. Supposing I have six circles and I want to evaluate every point in my grid, I have to calculate the distance between each of the grid points and each point on the six circles. This results in 100*6*100=60.000 calculations, without any refinement steps. I do not know how fast Matlab can do this..

Do you know any efficient way to do this?

Thanks!

Answer 3

[Jan]

A simple grid search with loops - just for demonstration:

% Define circles:
nCircle = 10;
Center  = rand(nCircle, 2);
Radius  = rand(nCircle, 1) / 5;

% To show that result is not unique:
% nCircle = 3;
% Center  = repmat(0.5, nCircle, 2);
% Radius  = [0.5, 1, 1.5];

% Get distances on a 50x50 grid:
nGrid = 50;
gridv = linspace(0, 1, nGrid);
Dist = zeros(nGrid, nGrid);
axes('NextPlot', 'add');
view(3)
for ix = 1:nGrid
    x = gridv(ix);
    for iy = 1:nGrid
        y = gridv(iy);
        D = 0;
        for iC = 1:nCircle
            D = D + abs(sqrt((x - C(iC, 1))^2 + (y - C(iC, 2))^2) - R(iC));
        end
        Dist(ix, iy) = D;
    end
end

% Draw sum of distances:
surf(gridv, gridv, Dist.' - min(Dist(:)), 'FaceAlpha', 0.5);

% Draw the circles
v = linspace(0, 2*pi, 50);
for iC = 1:nCircle
    line(C(iC, 1) + sin(v) * R(iC), C(iC, 2) + cos(v) * R(iC));
end

Grafics for 3 concentric circles:

[Bas Wagemaker]

Thanks!

I like the visual representation of the distance! I will certainly look into your example code to implement it in mine.

Have a nice evening!