See also Intersection of two line segments.
Code: Select all
; PB 5.42, cross-platform
EnableExplicit
Structure PointD
x.d
y.d
EndStructure
Structure CircleD
x.d ; x-coordinate of center
y.d ; y-coordinate of center
r.d ; radius
EndStructure
Procedure.i Intersect_Circles (*c1.CircleD, *c2.CircleD, *ps1.PointD, *ps2.PointD)
; in : c1, c2: two circles
; out: ps1, ps2 : intersection point(s), if existing
; return value: number of intersection points
; (-1 if the circles are coincident)
;
; [after Paul Bourke:
; <http://stackoverflow.com/questions/3349125/circle-circle-intersection-points>, 2016-04-17]
Protected.d dx, dy, d, a, h
Protected pp.PointD
dx = *c2\x - *c1\x
dy = *c2\y - *c1\y
d = Sqr(dx*dx + dy*dy) ; distance between the centers of the circles
If d > *c1\r + *c2\r ; the circles are separate
ProcedureReturn 0
EndIf
If d < Abs(*c1\r - *c2\r) ; one circle is contained within the other
ProcedureReturn 0
EndIf
If d = 0 And *c1\r = *c2\r
ProcedureReturn -1 ; the circles are coincident
EndIf
a = (*c1\r * *c1\r - *c2\r * *c2\r + d*d) / (2*d)
h = Sqr(*c1\r * *c1\r - a * a)
pp\x = *c1\x + a * (*c2\x - *c1\x) / d
pp\y = *c1\y + a * (*c2\y - *c1\y) / d
*ps1\x = pp\x + h * (*c2\y - *c1\y) / d
*ps1\y = pp\y - h * (*c2\x - *c1\x) / d
*ps2\x = pp\x - h * (*c2\y - *c1\y) / d
*ps2\y = pp\y + h * (*c2\x - *c1\x) / d
; In theory the following line should be "If h = 0",
; but we have to take floating point inaccuracy into account here.
If h < 0.000001
ProcedureReturn 1
Else
ProcedureReturn 2
EndIf
EndProcedure
CompilerIf #PB_Compiler_IsMainFile
; -- Demo
Define.CircleD c1, c2
Define.PointD ps1, ps2
Define.d a11, a12, a21, a22
Define.i numIs
Define.i offX = 250, offY = 150
; The two circles
c1\x = 80 : c1\y = 100 : c1\r = 60
c2\x = 133 : c2\y = 85 : c2\r = 40
If OpenWindow(0, 0, 0, 500, 500, "Circle-circle intersection") = 0
MessageRequester("Fatal error", "Program terminated.")
End
EndIf
CanvasGadget(0, 0, 0, 500, 500)
If StartVectorDrawing(CanvasVectorOutput(0))
; a1) Show original circles
VectorSourceColor(RGBA(0, 0, 255, 255))
AddPathCircle(c1\x, c1\y, c1\r)
AddPathCircle(c2\x, c2\y, c2\r)
StrokePath(2)
; a2) Calculate and show intersection points
numIs = Intersect_Circles(@c1, @c2, @ps1, @ps2)
Debug Str(numIs) + " intersection point(s)"
If numIs > 0
VectorSourceColor(RGBA(0, 0, 0, 255))
AddPathCircle(ps1\x, ps1\y, 5)
If numIs = 2
AddPathCircle(ps2\x, ps2\y, 5)
EndIf
FillPath()
EndIf
VectorSourceColor(RGBA(255, 0, 0, 255))
; b) Apply FillPath() to both complete circles
AddPathCircle(c1\x+offX, c1\y, c1\r)
AddPathCircle(c2\x+offX, c2\y, c2\r)
FillPath()
; Central angle of circle #1
a11 = Degree(ATan2(ps1\x-c1\x, ps1\y-c1\y))
a12 = Degree(ATan2(ps2\x-c1\x, ps2\y-c1\y))
; Central angle of circle #2
a21 = Degree(ATan2(ps1\x-c2\x, ps1\y-c2\y))
a22 = Degree(ATan2(ps2\x-c2\x, ps2\y-c2\y))
; c) Show intersection
c1\y + offY
c2\y + offY
AddPathCircle(c1\x, c1\y, c1\r, a11, a12) ; overlapping arc of circle #1
AddPathCircle(c2\x, c2\y, c2\r, a22, a21) ; overlapping arc of circle #2
FillPath()
; d) Show union
AddPathCircle(c1\x+offX, c1\y, c1\r, a12, a11) ; separate arc of circle #1
AddPathCircle(c2\x+offX, c2\y, c2\r, a21, a22) ; separate arc of circle #2
FillPath()
; e) Show c1 - c2
c1\y + offY
c2\y + offY
AddPathCircle(c1\x, c1\y, c1\r, a12, a11) ; separate arc of circle #1
AddPathCircle(c2\x, c2\y, c2\r, a22, a21) ; overlapping arc of circle #2
FillPath()
; f) Show c2 - c1
AddPathCircle(c1\x+offX, c1\y, c1\r, a11, a12) ; overlapping arc of circle #1
AddPathCircle(c2\x+offX, c2\y, c2\r, a21, a22) ; separate arc of circle #2
FillPath()
StopVectorDrawing()
EndIf
Repeat
Until WaitWindowEvent() = #PB_Event_CloseWindow
CompilerEndIf