find position using three beacons

General FreeBASIC programming questions.
paul doe
Moderator
Posts: 1832
Joined: Jul 25, 2017 17:22
Location: Argentina
Contact:

Re: find position using three beacons

Post by paul doe »

BasicCoder2 wrote: Nov 26, 2024 19:58
paul doe wrote: Nov 25, 2024 0:59
BasicCoder2 wrote: Nov 24, 2024 23:34 It only uses two beacons which is sufficient unless both are at 180 degrees apart. In practice the viewer can use any two beacons that have a difference of less than 180 degrees.
Indeed, two beacons suffices in this case. Two lines will eventually intersect in a 2D plane, unless they're parallel or coincident. You might want to keep the third beacon, and add code to select the two most appropriate for the calculation...
After spending perhaps too much time looking for a trig solution I have to bow down to your math savvy and agree. The intersect function is an easy simple and elegant solution.
Perhaps to minimize the error (albeit in a crude way), you can compute the three intersection points and then average their positions...

Here's an interesting paper on the subject:
https://www.researchgate.net/publicatio ... tersection
dodicat
Posts: 8171
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: find position using three beacons

Post by dodicat »

Hi basiccoder2.
This is the same as plane or coastal sailing on a ship.
The azimuths (bearings) of two known landmarks are taken via the compass and azimuth mirror.
On the chart, the bearings are taken from the compass rose on the chart and moved to the landmarks via parallel ruler, two lines are drawn and the intersection of these lines marks the position of the ship.
Plane sailing means you regard your immediate area as a plane and not a surface of a sphere.
For trans ocean sailing, then you have to consider you are actually traversing a sphere, so the navigation is more complex.
Anyway, using fb, and a cartesian type screen, and realizing that not all random bearings with random landmarks are viable.
press space to refresh and esc to end.

Code: Select all


Type Point
    As double x,y
End Type

Type Line
    As Point s,f
End Type

Function isleft(L As Line,p As Point) As Long
    Return  -Sgn((L.s.x-L.f.x)*(p.y-L.f.y)-(p.x-L.f.x)*(L.s.y-L.f.y))
End Function

Function intersects(L1 As Line,L2 As Line) As Long
    If isleft(L2,L1.s) = isleft(L2,L1.f) Then Return 0 
    If isleft(L1,L2.s) = isleft(L1,L2.f) Then Return 0
    Return -1
End Function


Function intersections(l1 As Line,l2 As Line,_out As Point) As Long
    Dim As Point p1=l1.s
    Dim As Point p2=l1.f
    Dim As Point p3=l2.s
    Dim As Point p4=l2.f
    Var x12=p1.x-p2.x
    Var x34=p3.x-p4.x
    Var y12=p1.y-p2.y
    Var y34=p3.y-p4.y
    Var c=x12 * y34 - y12 * x34
    If (Abs(c) < 0.01) Then 
        Return 0
    Else
        Var a = p1.x * p2.y - p1.y * p2.x
        Var b = p3.x * p4.y - p3.y * p4.x
        Var x = (a * x34 - b * x12) / c
        Var y = (a * y34 - b * y12) / c
        _out=Type(x,y)
    End If
    Return 1
End Function

Function GetLine(x As Long,y As Long,angle As Double,length As Double)As Line
    Dim As Double x2,y2
    angle=angle*.0174532925199433  
    x2=x+length*Cos(angle)
    y2=y+length*Sin(angle)
    Dim As Line _out=Type<Line>((x,y),(x2,y2))
    Return _out
End Function

Function getpos(p1 As Point,b1 As Double,p2 As Point,b2 As Double) As Point
    Var z1=GetLine(p1.x,p1.y,360-(b1-90)+180,900)
    Var z2=GetLine(p2.x,p2.y,360-(b2-90)+180,900) 
    Dim As Point pt
    Var p=intersections(z1,z2,pt)
    If p And intersects(z1,z2)Then Return Type(pt.x,pt.y) Else Return Type(0,0)
End Function

Screen 19
Window (0,0)-(799,599)
randomize
#define range(f,l) Rnd*((l)-(f))+(f)
#define onscreen(p) p.x>0 and p.x<800 and p.y>0 and p.y<600

Dim As Point mypos,p1,p2
Dim As Double bearingp1,bearingp2
Do
    
    Do
        p1=Type(range(200,600),range(200,400))
        p2=Type(range(200,600),range(200,400)) 
        bearingp1=Rnd*360
        bearingp2=Rnd*360
        mypos=getpos(p1,bearingp1,p2,bearingp2)
    Loop Until mypos.x<>0 And mypos.y<>0 And onscreen(mypos)
    Circle (p1.x,p1.y),5,2,,,,f
    Circle (p2.x,p2.y),5,6,,,,f
    Circle(mypos.x,mypos.y),5,7,,,,f
    Line(mypos.x,mypos.y)-(p1.x,p1.y)
    Line(mypos.x,mypos.y)-(p2.x,p2.y)
    Color 2
    Print "bearing of green ";bearingp1
    Color 6
    Print "bearing of red   ";bearingp2
    Color 15
    Draw String(mypos.x+2,mypos.y),Str(Int(mypos.x))+","+Str(Int(mypos.y))
    
    Sleep
    Cls
Loop Until Inkey=Chr(27)
'do


 
dafhi
Posts: 1714
Joined: Jun 04, 2005 9:51

Re: find position using three beacons

Post by dafhi »

fantastic work, dodicat. will use for future projects

i made this for reference

Code: Select all

' 2024 Dec 2 update: show projected intersection to reflect thread topic

Type Point
    As double x,y
End Type

Type Line
    As Point s,f
End Type

Function isleft(L As Line,p As Point) As Long
    Return  -Sgn((L.s.x-L.f.x)*(p.y-L.f.y)-(p.x-L.f.x)*(L.s.y-L.f.y))
End Function

Function intersects(L1 As Line,L2 As Line) As Long
    If isleft(L2,L1.s) = isleft(L2,L1.f) Then Return 0 
    If isleft(L1,L2.s) = isleft(L1,L2.f) Then Return 0
    Return -1
End Function


Function intersections(l1 As Line,l2 As Line,_out As Point) As Long
    Dim As Point p1=l1.s
    Dim As Point p2=l1.f
    Dim As Point p3=l2.s
    Dim As Point p4=l2.f
    Var x12=p1.x-p2.x
    Var x34=p3.x-p4.x
    Var y12=p1.y-p2.y
    Var y34=p3.y-p4.y
    Var c=x12 * y34 - y12 * x34
    If (Abs(c) < 0.01) Then 
        Return 0
    Else
        Var a = p1.x * p2.y - p1.y * p2.x
        Var b = p3.x * p4.y - p3.y * p4.x
        Var x = (a * x34 - b * x12) / c
        Var y = (a * y34 - b * y12) / c
        _out=Type(x,y)
    End If
    Return 1
End Function

sub DrawLine( L as Line, col as ulong )
    line (L.s.x,L.s.y)-(L.f.x,L.f.y), col
end sub

#define range(f,l) Rnd*((l)-(f))+(f)


#define min( a, b)    iif( (a)<(b), (a), (b) )
#define max( a, b)    iif( (a)>(b), (a), (b) )

function clamp( in as double, hi as double = 1, lo as double = 0) as double
    return min( max(in, lo), hi ) '' 2023 June 12
End Function

    function triwave( i as single ) as single
      return abs( i - int(i) - .5 ) - .25 ' by Stonemonkey
    end function

    function _cchsv(h as single, s as single, v as single) as ubyte ' 2024 July 24
      var wave_hgt = s * v
      return 255.499 * (wave_hgt * (clamp(triwave(h)*6+.5)-1) + v)
    end function

function hsv( h as single=0, s as single=1, v as single=1, a as ubyte = 255 ) as ulong ' 2024 May 21
      return rgba( _
    _cchsv( h + 0/3, s,v ), _
    _cchsv( h + 2/3, s,v ), _
    _cchsv( h + 1/3, s,v ), a )
end function


var w = 800
var h = 600

screenres w,h, 32

Dim As Point  mypos
dim as line   L1,L2

    dim as double t_end = timer + 25

randomize

Do
    for i as long = 1 to 2
        dim as long bor = 120
        L1.s=Type( range(bor,w - bor), range(bor, h - bor) )
        L1.f=Type( range(bor,w - bor), range(bor, h - bor) )
        L2.s=Type( range(bor,w - bor), range(bor, h - bor) )
        L2.f=Type( range(bor,w - bor), range(bor, h - bor) )
        dim as ulong col = hsv(rnd, .67, .5 * (1+rnd))
        Drawline L1, col
        Drawline L2, col
        intersections L1,L2, mypos
        circle (mypos.x, mypos.y), 3.5, iif( intersects( L1,L2 ), rgb(255,255,255), col ),,,,f
    next
    Sleep 1200
    Cls
Loop Until Inkey=Chr(27) or timer > t_end

locate 2,2
? "fin!"

sleep 1500
dodicat
Posts: 8171
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: find position using three beacons

Post by dodicat »

Thanks Dafhi.
I can only imagine Aussie, we've just had a big storm here, the sky is permanently laden grey, the temperature barely above zero, daylight ends about four and starts about 8 the next morning.
I watch Aussie gold hunters and Aussie opal diggers on the telly, only to get a glimpse of the glimmering bush heat and greenery.
So, among triangulating points for a robot, and messing about in his shed, basiccoder2 has the Great Australian bush at his doorstep.
If triangulation (3 points) are used, then a slight code adjustment for basiccoder2's robot.
Bearing in mind that the triangulation uses the average of three positions (as paul doe talked about earlier)
press a key to refresh, escape to end.

Code: Select all


Type Point
    As Double x,y
End Type

Type Line
    As Point s,f
End Type

Function isleft(L As Line,p As Point) As Long
    Return  -Sgn((L.s.x-L.f.x)*(p.y-L.f.y)-(p.x-L.f.x)*(L.s.y-L.f.y))
End Function

Function intersects(L1 As Line,L2 As Line) As Long
    If isleft(L2,L1.s) = isleft(L2,L1.f) Then Return 0 
    If isleft(L1,L2.s) = isleft(L1,L2.f) Then Return 0
    Return -1
End Function


Function intersections(l1 As Line,l2 As Line,_out As Point) As Long
    _out=Type(0,0)
    Dim As Point p1=l1.s
    Dim As Point p2=l1.f
    Dim As Point p3=l2.s
    Dim As Point p4=l2.f
    Var x12=p1.x-p2.x
    Var x34=p3.x-p4.x
    Var y12=p1.y-p2.y
    Var y34=p3.y-p4.y
    Var c=x12 * y34 - y12 * x34
    If (Abs(c) < 0.01) Then 
        Return 0
    Else
        Var a = p1.x * p2.y - p1.y * p2.x
        Var b = p3.x * p4.y - p3.y * p4.x
        Var x = (a * x34 - b * x12) / c
        Var y = (a * y34 - b * y12) / c
        _out=Type(x,y)
    End If
    Return 1
End Function

Function GetLine(x As Long,y As Long,angle As Double,length As Double)As Line
    Dim As Double x2,y2
    angle=angle*.0174532925199433  
    x2=x+length*Cos(angle)
    y2=y+length*Sin(angle)
    Dim As Line _out=Type<Line>((x,y),(x2,y2))
    Return _out
End Function

Function getbearing(p As Point,x As Point) As Double
    Const pi=4*Atn(1)
    Dim As Line z1=(p,x)
    Var a2=(Atan2((z1.s.y-z1.f.y),(z1.s.x-z1.f.x))*180/pi)-90
    Return Iif(-a2<0,a2+90,-a2)
End Function

Function getpos(p1 As Point,b1 As Double,p2 As Point,b2 As Double,p3 As Point,b3 As Double) As Point
    Var z1=GetLine(p1.x,p1.y,360-(b1-90)+180,200)
    Var z2=GetLine(p2.x,p2.y,360-(b2-90)+180,200)
    Var z3=GetLine(p3.x,p3.y,360-(b3-90)+180,200)
    If intersects(z1,z2) And intersects(z2,z3) And intersects(z3,z1) Then
        
        Dim As Long pp,flag
        Dim As Point pt1,pt2,pt3
        
        pp=intersections(z1,z2,pt1)
        If pp then flag+=1
        
        pp=intersections(z2,z3,pt2)
        If pp Then flag+=1
        
        pp=intersections(z3,z1,pt3)
        If pp Then flag+=1
        
        If flag=3 Then
            Pset(pt1.x,pt1.y),4'optional
            Line -(pt2.x,pt2.y),4'optional
            Line -(pt3.x,pt3.y),4'optional
            Return Type((pt1.x+pt2.x+pt3.x)/3,(pt1.y+pt2.y+pt3.y)/3)
        Else 
            Return Type(0,0)
        End If
        
    End If
End Function

Screen 19
Window (0,0)-(799,599)
Randomize
#define range(f,l) Rnd*((l)-(f))+(f)
#define onscreen(p) p.x>0 and p.x<800 and p.y>0 and p.y<600

Dim As Point mypos,p1,p2,p3,x
Dim As Double bearingp1,bearingp2,bearingp3

Do
    
    Do
        Cls
        
        p1=Type(range(100,700),range(100,500))'three fixed beacons
        p2=Type(range(100,700),range(100,500)) 
        p3=Type(range(100,700),range(100,500))
        
        x=Type(range(100,700),range(100,500))'a position
        bearingp1=getbearing(p1,x)
        bearingp2=getbearing(p2,x)
        bearingp3=getbearing(p3,x)
        mypos=getpos(p1,bearingp1,p2,bearingp2,p3,bearingp3)
    Loop Until mypos.x<>0 And mypos.y<>0 And onscreen(mypos)
    
    Circle (p1.x,p1.y),5,2,,,,f
    Circle (p2.x,p2.y),5,3,,,,f
    Circle (p3.x,p3.y),5,4,,,,f
    Circle(mypos.x,mypos.y),1,15,,,,f
    Line(mypos.x,mypos.y)-(p1.x,p1.y)
    Line(mypos.x,mypos.y)-(p2.x,p2.y)
    Line(mypos.x,mypos.y)-(p3.x,p3.y)
    Color 2
    Print "bearing of green  ";bearingp1
    Color 3
    Print "bearing of blue   ";bearingp2
    Color 4
    Print "bearing of red    ";bearingp3
    Color 15
    Print "true random position "
    Print x.x,x.y
    Draw String(mypos.x-80,mypos.y),"calculated "+Str((mypos.x))+" , "+Str((mypos.y))
    
    Sleep
    
Loop Until Inkey=Chr(27)
'do


 
dodicat
Posts: 8171
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: find position using three beacons

Post by dodicat »

Thanks for that impressive picture from your childhood basiccoder2.
You shouldn't have taken your post down, people should be reminded that we are destroying our planet.
dafhi
Posts: 1714
Joined: Jun 04, 2005 9:51

Re: find position using three beacons

Post by dafhi »

people know. they just need solutions.

in real life i never ask for money for what i do. it's my solution and it comes naturally
paul doe
Moderator
Posts: 1832
Joined: Jul 25, 2017 17:22
Location: Argentina
Contact:

Re: find position using three beacons

Post by paul doe »

dodicat wrote: Dec 12, 2024 10:31 You shouldn't have taken your post down, people should be reminded that we are destroying our planet.
Indeed, what happened there, BasicCoder2? Unfortunately we don't make backups so often, so those posts are lost forever, and the following conversation is left 'hanging'. Any particular reason for the deletion?
BasicCoder2
Posts: 3942
Joined: Jan 01, 2009 7:03
Location: Australia

Re: find position using three beacons

Post by BasicCoder2 »

dodicat
paul doe

I judged it not relevant to a FreeBASIC coders forum. Perhaps it would be relevant to a nature lovers forum. So to sum it up, yes, I do have what is left of the Great Australian bush at my doorstep. There are still some amazing natural wonders still remaining all over the world.

Dodicat, I seem to remember that you are living in Scotland? I think once there was a geographical report on where all the FreeBASIC users lived?
paul doe
Moderator
Posts: 1832
Joined: Jul 25, 2017 17:22
Location: Argentina
Contact:

Re: find position using three beacons

Post by paul doe »

BasicCoder2 wrote: Dec 13, 2024 3:27 I judged it not relevant to a FreeBASIC coders forum.
Yeah but that iguana was cute nonetheless...
Post Reply