Distance between a point and the intersection of two spheres

Hello,

Im trying to filter a sampled 3D human body signal. In sumary thanks to a motion capture system I get the coordinates of the shoulder (A) elbow (B) and wrist (C). the problem is that the distances mesured between AB (Humerus) and BC (Forearm) are not constant. To solve this I'm using the A and C samples as true and relocate the elbow. My idea is to see the points on a sphere of radius R1 (known humerus length) on the shoulder and R2 (known forearm length) on the wrist intersect. This intersection of the two spheres gives a circunference of possible possitions of the elbow. The point of this circunference which is closer to the sample B would be the true position of the elbow.

I've gotten to see the intersection petween both spheres but i couldn't get the points that create it. if anyone could help me it would be amazing.

Thanks in advance.

H=1.80;
humerus=0.179*H;
forearm=0.151*H;

theta=0:pi/80:2*pi;
shi=0:pi/80:pi;
r1=humerus;
c1=[0.2530    2.0740    0.2370];
[THETA, SHI]=meshgrid(theta,shi);
x1=c1(1)+r1*sin(SHI).*cos(THETA);
y1=c1(2)+r1*sin(SHI).*sin(THETA);
z1=c1(3)+r1*cos(SHI);
surf1=surf(x1,y1,z1);
r2=forearm;
c2=[0.2820    2.0920   -0.2870];
x2=c2(1)+r2*sin(SHI).*cos(THETA);
y2=c2(2)+r2*sin(SHI).*sin(THETA);
z2=c2(3)+r2*cos(SHI);
surf2=surf(x2,y2,z2);

Answer 1

[John D'Errico]

This is confusing, possibly because of how you are trying to solve it. I think you are working too hard.

The intersection between two spheres is a circle. Unless you are just trying to plot the spheres, there is no reason to generate them completely. Just find the equation of the circle. MAKE THINGS SIMPLE. Don't overly complicate them.

H=1.80;
humerus=0.179*H;
forearm=0.151*H;
r1=humerus;
c1=[0.2530    2.0740    0.2370];
r2=forearm;
c2=[0.2820    2.0920   -0.2870];

I presume that c1 is the center of sphere 1, and r1 the corresponding radius. Likewise, c2 and r2 are the center and radius of the second sphere.

First, do the spheres intersect? That is easy. I'll call D the distance between centers of the spheres, What mattes is the sum of the sphere radii.

D = norm(c1 - c2)
D =
      0.52511
      
r1 + r2
ans =
        0.594
        
(r1 + r2) > D
ans =
  logical
   1

If that sum exceeds the distance between spheres, then there is an intersection. If the sum of radii is EXACTLY equal to the distance between spheres, then the two spheres touch at a single point.

We also need to consider the case where one sphere lies entirely inside the other? How can that happen? For that to happen, the radius of the larger sphere must be at least as large as the sum of the distance between spheres, and the SMALLER radius.

max(r1,r2) >= (D + min(r1+r2))
ans =
  logical
   0

Had that been true, then one sphere would live entirely insode the other.

Given ll that, what is the equation of the circle of intersection? Also easy enough.

The center of that circle must lie on the line segment drawn between the spheres. How far? Assume the distance of the center of the circle along that line segment is d1 from the center of sphere 1. Then a little algebra tells me that

d1 = 1/2*(D + (r1^2 - r2^2)/D)
d1 =
      0.29106

So, we now that the circle has center at a distance of 0.2916 units along the line segment, starting at the point c1.

circlecenter = c1 + d1*(c2 - c1)/norm(c2 - c1)
circlecenter =
      0.26907        2.084    -0.053446

The radius of the circle of intersection easy, it just comes from the Pythagorean theorem.

circleradius = sqrt(r1^2 - d1^2)
circleradius =
      0.13819

What else do you need to compute? The circle lives in the plane perpendicular to that line connecting the two sphere centers.

It is all just basic geometry, with an application or two of the Pythagorean theorem.

[JAVIER GARCÍA DÍAZ]

Hello!!

First of all thank you for the answer, You were right I didn't need to generate both spheres.

The next step would be to check the distance between the measured elbow and the circunference where it is suposed to be. The closest point of the circunference to the measured one will be where the elbow is.

The location of the measured point is [0.2990 2.1120 -0.0310].

Again, thank you very much. It was of great help.

[John D'Errico]

The distance between the measure point and the circumferance? Still confusingly stated, manking me guess.

I derived what is essentially a circle that lives around some center point, in a plane

Given the point [0.2990 2.1120 -0.0310], now are you asking to find the closest point on that circle?

A simplistic solution would just have you generate say 1000 (or even 10000) points on the circumference of the circle, then pick the point that was closest to your given point. That would surely be adequate for your purposes. It would not be that difficult to solve the problem analytically either. I suppose I would just find the closest point on the plane of the circle. Then I would just shift the point so it lies at the known radius of the circle as a distance from the center. (In my case, I could also use some code I wrote and hanve posted on the file exchange to find the closest point on a general curve. That is not necessary here though.)

Anyway, you have a point, I'll call it P, and a description of the circle in terms of a center, a radius, and the plane it lives in. I would just do this:

H=1.80;
humerus=0.179*H;
forearm=0.151*H;
r1=humerus;
c1=[0.2530    2.0740    0.2370];
r2=forearm;
c2=[0.2820    2.0920   -0.2870];

% as above
D = norm(c1 - c2);
d1 = 1/2*(D + (r1^2 - r2^2)/D);
circlecenter = c1 + d1*(c2 - c1)/norm(c2 - c1);
circleradius = sqrt(r1^2 - d1^2);

% The plane of the circle is spanned by the unit basis vectors:
V = null(c1 - c2)
V =
    -0.034279      0.99789
      0.99889     0.032416
     0.032416      0.05634
     
% points on the circle
N = 1000;
t = linspace(0,2*pi,N)';

circlepoints = circlecenter + circleradius*[cos(t), sin(t)]*V';

% elbonian location. (Must an elbow always live in elbonia? only for Dilbert.)
P = [0.2990   2.1120   -0.0310];

% the closest point on the circle. Again, this is only approximate,
% since I chose only a finite number of points.
[mindist,minind] = min(sqrt(sum((circlepoints - P).^2,2)))
mindist =
     0.097864
minind =
   137
   
circlepoints(minind,:)
ans =
      0.37005       2.1779    -0.044631
      

Thus the distance to the perimeter of the circle is .098, and the closest point on the circle is found as circlepoints(minind,:).

I doubt I really need to show how to find that more accurately for your application.