Android Question area of polygon

Discussion in 'Android Questions' started by RB Smissaert, Jun 6, 2019.

  1. derez

    derez Expert Licensed User

    Thanks for testing, emexes ! The java source for the Surfacearea2 method is following, if you find anything wrong just point it out:
    Code:
    /**
          * The methods calculates the area on the Earth
    's spherical enclosed inside a polygon.
          * n is the polygon's size.
          * GindxLat and GindxLong are the indices of the parameters Lat and Long in the array.
          * The returned area 
    is in square Km.
          */
         
    public double SphericalArea2(double [][] Poligon, int n,
                 int GindxLat, int GindxLong)
            {
                double meanlat = 
    0;
                double [][] points =  new double[n][
    2] ; // new array because the content must change
                
    for (int i = 0;i<n;i++)
                {
                    points[i][
    0] = Poligon[i][GindxLat]* deg2rad ; // allocation by the index
                    points[i][
    1] = Poligon[i][GindxLong]* deg2rad ;
                    meanlat += points[i][
    0];
                
    }
                meanlat = meanlat/n ;
                return Spherical(points , n, 0, 1, meanlat);
            }
         
         private double Spherical(double [][] points , int n, int UindxX, int UindxY, double meanlat)
           {
                 double sum1,sum2,Rs,phis  ;
             
                double es = 0.0067394967423 ;  // eccentricity
                double twoPi = 2 * Math.PI ;
                phis = Math.atan(Math.tan(meanlat)/(1+es));
                Rs = 6378.137 / Math.pow((1+es*Math.sin(phis)*Math.sin(phis)),0.5);  // km
                // find angles for point 0
                double f1 = Bearing(points[0][0],points[0][1],points[1][0],points[1][1]);
                double f2 = Bearing(points[0][0],points[0][1],points[n-1][0],points[n-1][1]);
                sum1 = (f1 + twoPi - f2)%twoPi;
                sum2 = (f2 + twoPi - f1)%twoPi;
                // find angles for point n-1
                f1 = Bearing(points[n-1][0],points[n-1][1],points[0][0],points[0][1]);
                f2 = Bearing(points[n-1][0],points[n-1][1],points[n-2][0],points[n-2][1]) ;
                sum1 += (f1 + twoPi - f2)%twoPi;
                sum2 += (f2 + twoPi - f1)%twoPi;
                // find angles for the intermediate points
                for (int i = 1;i<n-1;i++)
                {
                f1 = Bearing(points[i][0],points[i][1],points[i+1][0],points[i+1][1]);
                f2 = Bearing(points[i][0],points[i][1],points[i-1][0],points[i-1][1]);
                sum1 += (f1 + twoPi - f2)%twoPi;
                sum2 += (f2 + twoPi - f1)%twoPi;
                }
                sum1 =  Math.min(sum1,sum2) ; // the minimun is the internal angles sum
                return Math.abs((Math.abs(sum1)-(n-2)*Math.PI)) * Rs*Rs ;// in sqr km
           }
           
     // This method uses arguments in Radians and returns Radians !
            private double Bearing(double LatR1,double LongR1,
                                   double LatR2,double LongR2)
            {
                double  range,bearing;
                double L = LongR1 - LongR2 ;
                double D = (Math.sin(LatR1) * Math.sin(LatR2)) + 
                           (Math.cos(LatR1) * Math.cos(LatR2) *Math.cos(L ));
                range = Math.acos(D);
                double C = (Math.sin(LatR2) - (Math.sin(LatR1) * D)) / 
                           (Math.cos(LatR1) * Math.sin(range));
                if (Math.sin(L) < 0) 
                    {bearing = Math.acos(C) ;}
                else 
                    {bearing = (2*Math.PI - Math.acos(C))  ;}
                return  bearing ;
            }
    The code can be easily translated to enable b4a debug.
     
    Last edited: Jun 9, 2019
  2. emexes

    emexes Well-Known Member Licensed User

    That's what I thought, too... I've translated it, but it's coming back with massively huge areas which I think has something to do with this bit:
    Code:
    1            sum1 += (f1 + twoPi - f2)%twoPi;
    2            sum2 += (f2 + twoPi - f1)%twoPi;
    3            }
    4            sum1 =  Math.min(sum1,sum2) ; // the minimun is the internal angles sum
    5            return Math.abs((Math.abs(sum1)-(n-2)*Math.PI)) * Rs*Rs ;// in sqr km
    where, as I understand it, on lines 1 and 2: sum1 and sum2 are constrained to 0..twoPi by the % mod operation,
    and thus on line 4, sum1 must thus also similarly be within 0..twoPi,
    but then on line 5 when we subtract (n-2)*Pi, I end up with a large number is then multiplied by the Rs squared

    But I've been using nav.GeoDistanceBearing to calculate the bearing; now that you've posted the Bearing function, I'll use that instead... like perhaps I got the units or zero base angle wrong.

    The only thing I've worked out thus far is that it doesn't handle repeated points, which is fair enough because it would be pretty difficult to calculate a bearing from one position to the same position.
     
  3. emexes

    emexes Well-Known Member Licensed User

    which is probably this line:
    Code:
    double C = (Math.sin(LatR2) - (Math.sin(LatR1) * D)) / (Math.cos(LatR1) * Math.sin(range));
    where if either (or both!) of cos(LatR1) or sin(range) evaluate to zero, then it's game over.
     
  4. derez

    derez Expert Licensed User

    I'm still debugging the translation... of course the NaN comes from division by range = 0
     
  5. RB Smissaert

    RB Smissaert Well-Known Member Licensed User

    I wouldn't post a complex polygon as I know it can't handle that.
    I think if the code in the library falls over repeated points then it should just check for these points and skip them.
    I added some code to my app to clear repeated points.

    Calculating radius from latitude is quite simple and fast, so added that. I take the latitude from the calculated centroid of the area.
    With that the code to calculate the area, posted by Klaus seems fine and has the benefit that there is no library, so you can see what is happening.
    Sofar all working fine.

    Code:
    Sub GetPolygonAreaSimple(lstPoints As ListAs Double
     
        
    Dim i As Int
        
    Dim Lat1, Lat2, Lng1, Lng2 As Double
        
    Dim dArea, LatScale, LngScale As Double
        
    Dim dRadius As Double
        
    Dim gpCentroid As OSMDroid_GeoPoint

        
    'Get the areas.
        Dim oGeoPoint1 As OSMDroid_GeoPoint = lstPoints.Get(0)
     
    'dNavA = 6378137
     If bCalculateRadius Then 'this looks the better way
      gpCentroid = GetAreaCentroid(lstPoints)
      dRadius = GetEarthRadius(gpCentroid.Latitude)
      LatScale =     
    cPI * dRadius / 180        ' converts latitude degrees to m
     Else
      LatScale =     
    cPI * dNavA / 180        ' converts latitude degrees to m
     End If

        LngScale = LatScale * 
    CosD(oGeoPoint1.Latitude) ' converts longitude degrees to m
        Lat1 = oGeoPoint1.Latitude * LatScale
        Lng1 = oGeoPoint1.Longitude * LngScale
     
        
    For i = 1 To lstPoints.Size - 1
            
    Dim oGeoPoint2 As OSMDroid_GeoPoint = lstPoints.Get(i)
            Lat2 = oGeoPoint2.Latitude * LatScale
            Lng2 = oGeoPoint2.Longitude * LngScale
            dArea = dArea + (Lng2 - Lng1) * (Lat2 + Lat1) / 
    2
            Lat1 = Lat2
            Lng1 = Lng2
      
    'General.RunLog("GetPolygonAreaSimple, " & i & ": " & oGeoPoint2.Latitude & "," & oGeoPoint2.Longitude)
        Next
        
     
    Return Round2(Abs(dArea / 1000000), 4)
     
    End Sub

    Sub GetEarthRadius(dLatitude As Double) As Double
     
     
    Dim dRadius As Double
     
    Dim dRadiusP As Double = 6356752 'at poles
     Dim dRadiusE As Double = 6378137 'at equator
     
     
     dRadius = 
    Sqrt((Power(dRadiusE * dRadiusE * Cos(dLatitude), 2) + Power(dRadiusP * dRadiusP * Sin(dLatitude), 2)) / _
        (
    Power(dRadiusE * Cos(dLatitude), 2) + Power(dRadiusP * Sin(dLatitude), 2)))
        
     
    Return dRadius
     
    End Sub

    Sub GetAreaCentroid(lstPoints As ListAs OSMDroid_GeoPoint

        
    Dim i As Long
        
    Dim second_factor As Double
        
    Dim polygon_area As Double
        
    Dim X As Double
        
    Dim Y As Double
     
    Dim oGeoPointCenter As OSMDroid_GeoPoint

        
    For i = 0 To lstPoints.Size - 2
      
    Dim oGeoPoint1 As OSMDroid_GeoPoint = lstPoints.Get(i)
      
    Dim oGeoPoint2 As OSMDroid_GeoPoint = lstPoints.Get(i + 1)
            second_factor = _
            oGeoPoint1.Longitude * oGeoPoint2.Latitude - _
                            oGeoPoint2.Longitude * oGeoPoint1.Latitude
            X = X + (oGeoPoint1.Longitude + oGeoPoint2.Longitude) * _
                second_factor
            Y = Y + (oGeoPoint1.Latitude + oGeoPoint2.Latitude) * _
                second_factor
        
    Next

        
    'Divide by 6 times the polygon's area.
        polygon_area = GetPolygonArea(lstPoints)
     General.RunLog(
    "GetAreaCentroid, polygon_area: " & polygon_area)
     
        X = X / 
    6 / polygon_area
        Y = Y / 
    6 / polygon_area
     
     General.RunLog(
    "GetAreaCentroid, Y: " & Y)
     General.RunLog(
    "GetAreaCentroid, X: " & X)

        
    'If the values are negative, the polygon is
        'oriented counterclockwise. Reverse the signs.
     '(this works for me, but maybe wrong for other locations)
     If Y < 0 Then
      X = -X
      Y = -Y
     
    End If
     
     oGeoPointCenter.Initialize(Y, X)

     
    Return oGeoPointCenter

    End Sub

    'this is not the actual area in km2
    Sub GetPolygonArea(lstPoints As ListAs Double

        
    Dim i As Int
        
    Dim dArea As Double

        
    'Get the areas.
        For i = 0 To lstPoints.Size - 2
      
    Dim oGeoPoint1 As OSMDroid_GeoPoint = lstPoints.Get(i)
      
    Dim oGeoPoint2 As OSMDroid_GeoPoint = lstPoints.Get(i + 1)
            dArea = dArea + _
                    (oGeoPoint2.Longitude - oGeoPoint1.Longitude) * _
                    (oGeoPoint2.Latitude + oGeoPoint1.Latitude) / 
    2
     
    Next

     
    Return Abs(dArea)

    End Sub

    RBS
     
  6. emexes

    emexes Well-Known Member Licensed User

    My translation from Java to B4A thus far is below. I did trip over another little anomaly in the Bearing function, which could be of my own doing by caching intermediate Sin/Cos values into temporary variables, probably loses a few low-end bits of precision in the process.

    But my main problem is still that I am still getting large numbers out of Abs(Abs(MinSum) - (N - 2) * cPI) on the Spherical() return line.

    My other looming problem is that it is nighttime here and I am starting to go a bit fuzzy, so I'm letting go until tomorrow.

    :)

    Code:
    Sub Bearing(LatR1 As Double, LongR1 As Double, LatR2 As Double, LongR2 As Double) As Double
     
        
    Dim SinLatR1 As Double = Sin(LatR1)
        
    Dim SinLatR2 As Double = Sin(LatR2)
        
    Dim CosLatR1 As Double = Cos(LatR1)
        
    Dim CosLatR2 As Double = Cos(LatR2)
     
        
    Dim L As Double = LongR1 - LongR2
        
    Dim D As Double = (SinLatR1 * SinLatR2) + (CosLatR1 * CosLatR2 * Cos(L))
        
    Dim Range As Double = ACos(D)
        
    Dim C As Double = (SinLatR2 - (SinLatR1 * D)) / (CosLatR1 * Sin(Range))
     
        
    'sometimes C falls very slightly outside of the range -1..1
        'eg C = -1.0000000000000997 and thus ACos(C) = NAN
     
        
    If C < -1 Then
            C = -
    1
        
    else if C > 1 Then
            C = 
    1
        
    End If
                   
        
    Dim ReturnValue As Double
     
        
    If Sin(L) < 0 Then
            ReturnValue = 
    ACos(C)
        
    Else
            ReturnValue = 
    2 * cPI - ACos(C)
        
    End If
     
        
    Return ReturnValue

    End Sub
     
    Sub Spherical(Points(,) As Double , N As Int, MeanLat As Double) As Double

        
    Dim es As Double = 0.0067394967423    'eccentricity
        Dim twoPi As Double = 2 * cPI
        
    Dim phis As Double = ATan(Tan(MeanLat) / (1 + es))
        
    Dim Sinphis As Double = Sin(phis)
        
    Dim Rs As Double = 6378.137 / Sqrt(1 + es * Sinphis * Sinphis)
        Rs = 
    6378.137    'close enough for time being
     
        
    Dim Sum1 As Double = 0
        
    Dim Sum2 As Double = 0
     
        
    For I = 0 To N - 1
            
    Dim J As Int = (I + 1Mod N
            
    Dim K As Int = (I + 2Mod N
       
            
    Dim F1 As Double = Bearing(Points(j, 0), Points(j, 1), Points(k, 0), Points(k, 1))
            
    Dim F2 As Double = Bearing(Points(j, 0), Points(j, 1), Points(i, 0), Points(i, 1))
       
            Sum1 = (Sum1 + F1 + twoPi - F2) 
    Mod twoPi
            Sum2 = (Sum2 + F2 + twoPi - F1) 
    Mod twoPi
        
    Next
     
        
    Dim MinSum As Double = Min(Sum1, Sum2)    'the minimum is the internal angles sum
     
        
    Return Abs(Abs(MinSum) - (N - 2) * cPI) * Rs * Rs    'in square km

    End Sub

    Sub SphericalArea2(Poligon(,) As Double, N As Int) As Double

        
    Dim MeanLat As Double = 0

        
    Dim Points(N, 2As Double    'new Array because the content must change
     
        
    For I = 0 To N - 1
            Points(I, 
    0) = Poligon(I, 0) * cPI / 180
            Points(I, 
    1) = Poligon(I, 1) * cPI / 180
            MeanLat = MeanLat + Points(I, 
    0)
        
    Next
     
        MeanLat = MeanLat / N
     
        
    Return Spherical(Points, N, MeanLat)
     
    End Sub
     
    Last edited: Jun 9, 2019
  7. emexes

    emexes Well-Known Member Licensed User

    <thumbs up> x 5

    Yea, stick with the flat area calculations. It'll be nice to get the spherical calculations working, but it's more a labour of love than a degree of precision actually required in the real world. Even NASA suggested that straight lines are often close enough, and that a line had to be over a hundred kilometres long before it was out by a pixel on their display compared to an absolutely correct line.

    If it's good enough for NASA, then it's good enough for me, and presumably good enough for your front garden landscape planning or whatever your mission is.

    :)

    Except for O-rings, let's not follow NASA's lead on those...

    :-/
     
  8. derez

    derez Expert Licensed User

    Checking the rectangle case in SphericalArea I found that a slight difference is required between two points Lat or Long:
    Code:
    Sub test2
        
        
    Dim p(4,2As Double

        p(
    0,0) = 50
        p(
    0,1) = 30.000001
        p(
    1,0) = 51
        p(
    1,1) = 30
        p(
    2,0) = 51
        p(
    2,1) = 31
        p(
    3,0) = 50.00001
        p(
    3,1) = 31
        
    Dim lst As List
        lst.Initialize
        
    For i = 0 To 3
            lst.Add(p(i,
    0) & "," & p(i,1))
            
    Log(lst.Get(i))
        
    Next

        
    Log(nav.SphericalArea2(p,4,0,1))
        
    Log(nav.SphericalArea(lst))
    End Sub
    I'll try a simplified bearing calculation and see if it makes a change.
     
  9. derez

    derez Expert Licensed User

    emexes, your note in the code drew my attention and I found that it causes all the problem. I added two lines to eliminate these cases and everything works fine.
    Here is the complete translated code, as an app in b4j:
    Code:
    Sub Process_Globals
    '    Private fx As JFX
    '    Private MainForm As Form
        Private deg2rad As Double = cPI/180

    End Sub

    Sub AppStart (Form1 As Form, Args() As String)
    '    MainForm = Form1
        'MainForm.RootPane.LoadLayout("Layout1") 'Load the layout file.
    '    MainForm.Show
        test2
    End Sub

    Sub test2

        
    Dim points(4,2As Double

        points(
    0,0) = 52.597316
        points(
    0,1) = -2.150112
        points(
    1,0) = 52.597316
        points(
    1,1) = -2.147188
        points(
    2,0) = 52.594396
        points(
    2,1) = -2.147188
        points(
    3,0) = 52.594396
        points(
    3,1) = -2.150112
        
    Log("Area = " & SphericalArea2(points,4,0,1))
      
    End Sub
      

    'Return true to allow the default exceptions handler to handle the uncaught exception.
    Sub Application_Error (Error As Exception, StackTrace As StringAs Boolean
        
    Return True
    End Sub


    '  The methods calculates the area on the Earth's spherical enclosed inside a polygon.
    '  n Is the polygon's size.
    '  GindxLat And GindxLong are the indices of the parameters Lat And Long in the Array.
    '  The returned area Is in square Km.
    Public Sub  SphericalArea2(Poligon(,) As Double,  n As Int, _
        GindxLat 
    As Int ,  GindxLong As Int) As Double

        
    Dim  meanlat As Double = 0
        
    Dim  points(n,2As Double ' new Array because the content must change
        For i = 0 To n-1

            points(i,
    0) = Poligon(i,GindxLat)* deg2rad ' allocation by the index
            points(i,1) = Poligon(i,GindxLong)* deg2rad
            meanlat = meanlat + points(i,
    0)
        
    Next
        meanlat = meanlat/n
        
    Return Spherical(points , n, meanlat)
    End Sub


    Sub  Spherical( points(,) As Double ,n As Int, meanlat As Double) As Double
      
        
    Dim sum1,sum2,phis As Double
     
        
    Dim es = 0.0067394967423  As Double  ' eccentricity
        Dim twoPi = 2 * cPI As Double
        phis = 
    ATan(Tan(meanlat)/(1+es))
        
    Dim Rs As Double = 6378.137 /Power((1 + es * Sin(phis) * Sin(phis)), 0.5' km

        
    ' find angles For point 0
        Dim f1 As Double = Bearing(points(0,0), points(0,1), points(1,0), points(1,1))
        
    Dim f2 As Double = Bearing(points(0,0),points(0,1),points(n-1,0),points(n-1,1))
        sum1 = (f1 + twoPi - f2) 
    Mod twoPi
        sum2 = (f2 + twoPi - f1) 
    Mod twoPi
        
    ' find angles For point n-1
        f1 = Bearing(points(n-1,0),points(n-1,1),points(0,0),points(0,1))
        f2 = Bearing(points(n-
    1,0),points(n-1,1),points(n-2,0),points(n-2,1))
        sum1 = sum1 + (f1 + twoPi - f2) 
    Mod twoPi
        sum2 = sum2 + (f2 + twoPi - f1) 
    Mod twoPi
        
    ' find angles For the intermediate points
        For i = 1 To n-2
            f1 = Bearing(points(i,
    0),points(i,1),points(i+1,0),points(i+1,1))
            f2 = Bearing(points(i,
    0),points(i,1),points(i-1,0),points(i-1,1))
            sum1 = sum1 + (f1 + twoPi - f2) 
    Mod twoPi
            sum2 = sum2 + (f2 + twoPi - f1) 
    Mod twoPi
        
    Next

        sum1 =  
    Min(sum1,sum2) ' the minimun Is the internal angles sum
        Return Abs((Abs(sum1)-(n-2)*cPI)) * Rs*Rs ' in sqr km
    End Sub
     
        
            
    ' This method uses arguments in Radians And returns Radians !
     Sub Bearing( LatR1 As Double, LongR1 As Double, LatR2 As Double, LongR2 As Double) As Double
          
        
    Dim  range, Bearing1 As Double
        
    Dim L As Double = LongR1 - LongR2
        
    Dim D As Double= (Sin(LatR1) * Sin(LatR2)) + (Cos(LatR1) * Cos(LatR2) *Cos(L ))

        D = 
    Max(Min(D,1),-1)   ' Addition
        range = ACos(D)
        
    Dim C As Double = (Sin(LatR2) - (Sin(LatR1) * D)) /(Cos(LatR1) * Sin(range))
        C = 
    Max(Min(C,1),-1)   ' Addition
        If (Sin(L) < 0Then
            Bearing1 = 
    ACos(C)
        
    Else
            Bearing1 = (
    2*cPI - ACos(C))
        
    End If

        
    Return  Bearing1
    End Sub
    When I'll remember how to update a library in Eclipse I'll update the library...
    The area of the small rectangle is much smaller than 1 kmr :
    Area = 0.06396440387541677
     
  10. derez

    derez Expert Licensed User

    Library updated to 1.2 !
     
    emexes likes this.
  11. RB Smissaert

    RB Smissaert Well-Known Member Licensed User

    Thanks, will try it out.
    Won't have to deal with duplicate points anymore.

    RBS
     
  12. RB Smissaert

    RB Smissaert Well-Known Member Licensed User

    Just tested and seems to work well.
    SphericalArea2 fine with small areas, eg 0.015 km2.
    Did you look at the accuracy of the 3 mentioned algo's?

    RBS
     
  13. derez

    derez Expert Licensed User

    Please reload the library, I made some more corrections of the same problem in other methods.
     
  14. RB Smissaert

    RB Smissaert Well-Known Member Licensed User

    Will do.

    RBS
     
  15. emexes

    emexes Well-Known Member Licensed User

    I had a similar discovery myself (a good sleep works wonders), in my translation of your code to:
    Code:
    Sum1 = (Sum1 + F1 + twoPi - F2) Mod twoPi
    Sum2 = (Sum2 + F2 + twoPi - F1) 
    Mod twoPi
    the opening bracket needs to move along one step, ie:
    Code:
    Sum1 = Sum1 + (F1 + twoPi - F2) Mod twoPi
    Sum2 = Sum2 + (F2 + twoPi - F1) 
    Mod twoPi
    and now I am getting results that very nearly match:

    12380.233223394947 from Navigation library
    12380.233251685318 from B4A translation

    which I'm almost happy with. Will try running same code with B4J (ie x86 FPU, JVM) instead of on Android device (ie ARM FPU, Android JVM) - I doubt that will fix it, but stranger things have happened. Also, I changed the order in which the polygon angles are summed, so I'll try reverting to your way (first angle, last angle, then the middle angles).

    What could possibly go wrong?

    :)
     
  16. emexes

    emexes Well-Known Member Licensed User

    12380.233251685318 from B4J translation, so that's kinda reassuring.
     
  17. emexes

    emexes Well-Known Member Licensed User

    Changng order of polygon angles summing doesn't change anything, which is very reassuring ;-)

    12380.233251685318 from B4A translation using original summing order
     
  18. emexes

    emexes Well-Known Member Licensed User

    Righto, final result:

    12380.233223394947 from Navigation library 1.1
    12380.233251685318 from my translation to B4X (sums polygon angles in the order that they are presented by the array)
    12380.233251685318 from your translation to B4X (sums polygon angles in first, last, remainder order)

    I'm prepared to say that the library result difference sure looks like a single-vs-double floating point issue. If pressed, my backup excuse is that the output can only be as good as the input, so we would expect accuracy in the order of the 7 digits of the supplied earth radius precision.
     
    Last edited: Jun 10, 2019
  19. emexes

    emexes Well-Known Member Licensed User

    Postscript: yeah, I should just let this go, instead of chasing down a less-than-a-millionth-of-a-percent discrepancy, but... :-?

    compiled and ran original Java code, still got very slightly different answer to that from Navigation library:

    12380.233223394947 from Navigation library 1.1
    12380.233251685318 from my translation to B4X (sums polygon angles in the order that they are presented by the array)
    12380.233251685318 from your translation to B4X (sums polygon angles in first, last, remainder order)
    12380.233251685318 from compiling and running original Java code on a Windows computer

    so... I don't know. The usual suspects would be:

    1/ a floating point literal in the library code is being default cast to single precision rather than double precision
    2/ the Math.pow(x. 0.5) comprises two precision-losing operations eg exp(log(x) * 0.5) rather than one square-root operation
    3/ deg2rad is single, not double (although I doubt that, given that Math.PI is a double, and the compiler promptly points out the loss of precision)
     
  20. emexes

    emexes Well-Known Member Licensed User

    BTW I tried this, and my problem test cases were magically fixed too.

    :)
     
    derez likes this.
Loading...
  1. This site uses cookies to help personalise content, tailor your experience and to keep you logged in if you register.
    By continuing to use this site, you are consenting to our use of cookies.
    Dismiss Notice