Intersection of two circles
Posted: Fri May 13, 2016 8:20 pm
The procedure Intersect_Circles() calculates the intersection point(s) of two circles (if existing) by using pure mathematics. The demo code uses PB's new VectorDrawing library.
See also Intersection of two line segments.
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
Thank you for pointing to this issue!