POV-Ray : Newsgroups : povray.general : Intersection of three spheres? : Re: Intersection of three spheres? Server Time
8 Aug 2024 20:23:53 EDT (-0400)
  Re: Intersection of three spheres?  
From: John VanSickle
Date: 20 Oct 2000 23:16:54
Message: <39F10BE1.FEE72508@erols.com>
Greg M. Johnson wrote:
> 
> If I have three spheres that intersect, how can I find the two points
> at which they intersect?

This macro will return one of the two intersection points you want, or
bomb out with an error if the spheres do not intersect in the way you
want; if you don't want it to bomb, then you can change it to do
something else.

I tested it and it seems to work.  There may be a simpler algorithm.

// START OF MACRO CODE

#macro IntersectThreeSpheres(pA,rA,pB,rB,pC,rC,fW)
// pX = center of sphere X
// rX = radius of sphere X
// fW= flag controlling which of two possible values to use
//     >0 if you want the upper one, <0 if the lower

#local vB=pB-pA;
#local vC=pC-pA;
#local vF=pC-pB;

#local vN=vcross(vC,vB);

#if(vlength(vN)=0) #error "Centers of circles lie on one line.\n" #end

#local vN=vnormalize(vN);

#local dAB=vlength(vB);
#local dAC=vlength(vC);
#local dBC=vlength(vF);

#if(dAB>rA+rB | rA>dAB+rB | rB>rA+dAB)
  #error "First and second spheres do not intersect.\n"
#end

#if(dAC>rA+rC | rA>dAC+rC | rC>rA+dAC)
  #error "First and third spheres do not intersect.\n"
#end

#if(dBC>rB+rC | rB>dBC+rC | rC>rB+dBC)
  #error "Second and third spheres do not intersect.\n"
#end

#local dKAB=(rA*rA-rB*rB+dAB*dAB)/2/dAB;
#local dKAC=(rA*rA-rC*rC+dAC*dAC)/2/dAC;
#local vB=vnormalize(vB);
#local vC=vnormalize(vC);

#local sD=vdot(vN,vcross(vB,vC));

#local sN=vdot(vN,pA);
#local sB=vdot(vB,pA)+dKAB;
#local sC=vdot(vC,pA)+dKAC;

#local pM=<vdot(<sN,vN.y,vN.z>,vcross(<sB,vB.y,vB.z>,<sC,vC.y,vC.z>))/sD,
           vdot(<vN.x,sN,vN.z>,vcross(<vB.x,sB,vB.z>,<vC.x,sC,vC.z>))/sD,
           vdot(<vN.x,vN.y,sN>,vcross(<vB.x,vB.y,sB>,<vC.x,vC.y,sC>))/sD>;

#debug concat("<",str(pM.x,0,3),",",str(pM.y,0,3),",",str(pM.z,0,3),">\n")

#local sD=vlength(pA-pM);
#local sD=rA*rA-sD*sD;
#if(sD<0) #error "Not the kind of intersection you're after.\n" #end
#local sD=sqrt(sD);

( (fW>0)?(pM+vN*sD):(pM-vN*sD) )

#end

//END OF MACRO CODE

Hope this helps,
John
(posted and e-mailed)
-- 
ICQ: 46085459


Post a reply to this message

Copyright 2003-2023 Persistence of Vision Raytracer Pty. Ltd.