# PureBasic Forum

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

 All times are UTC + 1 hour

 Page 1 of 1 [ 5 posts ]
 Print view Previous topic | Next topic
Author Message
 Post subject: Intersection of two circlesPosted: Fri May 13, 2016 8:20 pm

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.

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
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

If StartVectorDrawing(CanvasVectorOutput(0))
; a1) Show original circles
VectorSourceColor(RGBA(0, 0, 255, 255))
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))
If numIs = 2
EndIf
FillPath()
EndIf

VectorSourceColor(RGBA(255, 0, 0, 255))

; b) Apply FillPath() to both complete circles
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

_________________
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

 Post subject: Re: Intersection of two circlesPosted: Fri May 13, 2016 9:16 pm
 PureBasic Expert

Joined: Sun Apr 12, 2009 6:27 am
Posts: 3599
I like this kind of staff
Thanks LJ

flags for your OpenWindow next time

_________________
Egypt my love

Top

 Post subject: Re: Intersection of two circlesPosted: Sat May 14, 2016 12:00 am
 Enthusiast

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

 Post subject: Re: Intersection of two circlesPosted: Sat May 14, 2016 6:50 am

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

 Post subject: Re: Intersection of two circlesPosted: Sat May 14, 2016 10:01 am

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.

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

 Display posts from previous: All posts1 day7 days2 weeks1 month3 months6 months1 year Sort by AuthorPost timeSubject AscendingDescending
 Page 1 of 1 [ 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 forumYou cannot reply to topics in this forumYou cannot edit your posts in this forumYou cannot delete your posts in this forum

Search for:
 Jump to:  Select a forum ------------------ PureBasic    Coding Questions    Game Programming    3D Programming    Assembly Programming    The PureBasic Editor    The PureBasic Form Designer    General Discussion    Feature Requests and Wishlists    Tricks 'n' Tips Bug Reports    Bugs - Windows    Bugs - Linux    Bugs - Mac OSX    Bugs - IDE    Bugs - Documentation OS Specific    AmigaOS    Linux    Windows    Mac OSX Miscellaneous    Announcement    Off Topic Showcase    Applications - Feedback and Discussion    PureFORM & JaPBe    TailBite