It is currently Mon May 25, 2020 8:47 am

All times are UTC + 1 hour




Post new topic Reply to topic  [ 5 posts ] 
Author Message
 Post subject: Intersection of two circles
PostPosted: Fri May 13, 2016 8:20 pm 
Offline
Addict
Addict

Joined: Thu Jun 07, 2007 3:25 pm
Posts: 3845
Location: Berlin, Germany
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:
; 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

_________________
Please excuse my flawed English. My native language is PureBasic.
Search
RSBasic's backups


Last edited by Little John on Sat May 14, 2016 9:42 am, edited 1 time in total.

Top
 Profile  
Reply with quote  
 Post subject: Re: Intersection of two circles
PostPosted: Fri May 13, 2016 9:16 pm 
Offline
PureBasic Expert
PureBasic Expert

Joined: Sun Apr 12, 2009 6:27 am
Posts: 3599
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


Top
 Profile  
Reply with quote  
 Post subject: Re: Intersection of two circles
PostPosted: Sat May 14, 2016 12:00 am 
Offline
Enthusiast
Enthusiast
User avatar

Joined: Thu Mar 24, 2011 12:40 am
Posts: 536
Location: Iowa, USA
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.


Top
 Profile  
Reply with quote  
 Post subject: Re: Intersection of two circles
PostPosted: Sat May 14, 2016 6:50 am 
Offline
Addict
Addict

Joined: Fri Nov 09, 2012 11:04 pm
Posts: 1758
Location: Uttoxeter, UK
Very interesting. I'll have to study this carefully.
Works fine on my MacBook.

Thank you.

_________________
DE AA EB


Top
 Profile  
Reply with quote  
 Post subject: Re: Intersection of two circles
PostPosted: Sat May 14, 2016 10:01 am 
Offline
Addict
Addict

Joined: Thu Jun 07, 2007 3:25 pm
Posts: 3845
Location: Berlin, Germany
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:
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:
c2\x = 178.86859966642593616

resulted in 0 intersection points.

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

Code:
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:
If h = 0

with

Code:
If h < 0.000001

which works fine in my tests.
Many thanks again, BP!

_________________
Please excuse my flawed English. My native language is PureBasic.
Search
RSBasic's backups


Top
 Profile  
Reply with quote  
Display posts from previous:  Sort by  
Post new topic Reply to topic  [ 5 posts ] 

All times are UTC + 1 hour


Who is online

Users browsing this forum: ar-s, FlatEarth and 16 guests


You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum

Search for:
Jump to:  

 


Powered by phpBB © 2008 phpBB Group
subSilver+ theme by Canver Software, sponsor Sanal Modifiye