Intersection of two circles

Share your advanced PureBasic knowledge/code with the community.
Little John
Addict
Addict
Posts: 4519
Joined: Thu Jun 07, 2007 3:25 pm
Location: Berlin, Germany

Intersection of two circles

Post by Little John »

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.

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
Image
Last edited by Little John on Sat May 14, 2016 9:42 am, edited 1 time in total.
RASHAD
PureBasic Expert
PureBasic Expert
Posts: 4636
Joined: Sun Apr 12, 2009 6:27 am

Re: Intersection of two circles

Post by RASHAD »

I like this kind of staff :)
Thanks LJ

Well you please add (#PB_Window_SystemMenu| #PB_Window_ScreenCentered)
flags for your OpenWindow next time :P
Egypt my love
User avatar
BasicallyPure
Enthusiast
Enthusiast
Posts: 536
Joined: Thu Mar 24, 2011 12:40 am
Location: Iowa, USA

Re: Intersection of two circles

Post by BasicallyPure »

Hi,
Nice work.
Will it detect if the circles only touch at one point?
BasicallyPure
Until you know everything you know nothing, all you have is what you believe.
davido
Addict
Addict
Posts: 1890
Joined: Fri Nov 09, 2012 11:04 pm
Location: Uttoxeter, UK

Re: Intersection of two circles

Post by davido »

Very interesting. I'll have to study this carefully.
Works fine on my MacBook.

Thank you.
DE AA EB
Little John
Addict
Addict
Posts: 4519
Joined: Thu Jun 07, 2007 3:25 pm
Location: Berlin, Germany

Re: Intersection of two circles

Post by Little John »

Hi, thanks for all your kind words.
BasicallyPure wrote:Will it detect if the circles only touch at one point?
In theory, it did do so, but unfortunately I hadn't tested this situation. :-( Thank you for pointing to this issue!

I added a Debug line that tells the number of intersection points.

When I defined for instance the circles like this

Code: Select all

c1\x =  80 : c1\y = 100 : c1\r = 60
c2\x = 178.86859966642593615 : c2\y =  85 : c2\r = 40
I got 2 intersection points with the previous code.

Using instead

Code: Select all

c2\x = 178.86859966642593616
resulted in 0 intersection points.

When I used a value between the two values above, e.g.

Code: Select all

c2\x = 178.868599666425936165

then the compiler complained about "too many digits".
So it was not possible to get 1 intersection/touch point in this example with the previous version of the code.

Obviously, the problem is due to the well-known floating point inaccuracy.

In the procedure, I have now replaced

Code: Select all

If h = 0
with

Code: Select all

If h < 0.000001
which works fine in my tests.
Many thanks again, BP!
Post Reply