Trilateration issue

[MartinManganiello]

Hi, I am facing to a trilateration exercise, i'm doing this:

from pygeodesy.vector3d import trilaterate3d2, Vector3d
p1 = Vector3d(2, 2, 0)
p2 = Vector3d(3, 3, 0)
p3 = Vector3d(1, 4, 0)
r1 = 1
r2 = 1
r3 = 1.4142
trilaterate3d2(p1, r1, p2, r2, p3, r3)

When i call to trilaterate3d2 i get: pygeodesy.errors.IntersectionError: center1 (Vector3d(2, 2, 0)), center2 (Vector3d(3, 3, 0)), center3 (Vector3d(1, 4, 0)), radius1 (1), radius2 (1) or radius3 (1.4142): no intersection
But it should return P=(2,3,0)
What i'm doing wrong?

[mrJean1]

Which version of PyGeodesy are you using? In any case, the error means that the 3 circles do not intersect. It is numerically close, because for r3 = 1.41422, they do.

The distance between P1 and P3 is sqrt(5) or 2.236, those 2 circles do intersect. Same for P2 and P3. But there is no single intersection or overlap of all three circles.

The distance from (2, 3, 0) to P1 and P2 is 1. But to P3 it is sqrt(2) or 1.414213562373..., i.e. slightly larger than r3.

[MartinManganiello]

Thanks very much @mrJean1 ! One more question, Vector3D is an n-vector representing a normal to point on earth's surface. If my exercise is located in space, should I use the galactic coordinate system or like something? It is simply out of curiosity, your answer helped me a lot.

[mrJean1]

No, Vector3d is just a basic vector of x, y and z. It is used in several places in PyGeodesy mostly to simplify vector operations like addition, cross product, length, etc.

Take a look at the Nvector class in the sphericalNvector and ellipsoidalNvector modules. See also the Ned class in ellipsoidalNvector.

Also, local X, Y, Z vectors like ENU are handled by the EcefCartesian class in the ecef module. There are quite a few documentation mistakes in that class, sorry. An improved version is forthcoming in the next release.

Not sure which class is best for a galactic environment. Maybe sphericalNvector will suffice since at galactic distances the ellipsoidal/flattening accuracy offered by ellipsoidalNvector may not matter.

[MartinManganiello]

Hi @mrJean1, I want to show you this last example:

from pygeodesy.vector3d import trilaterate3d2, Vector3d
p1 = Vector3d(-500, -200, 0)  
p2 = Vector3d(100, -100, 0)  
p3 = Vector3d(500, 100, 0)  
r1 = 450 
r2 = 694.6221994724903  
r3 = 1011.1874208078342  
trilaterate3d2(p1, r1, p2, r2, p3, r3)

This gives me a no intersection error.
I drew the circles with matplotlib:

All the circles intersect the black point.
Also use this code to get the location of the point by trilateration:

def trilaterate(dataset: list) -> list:
    # dataset to variables
    x1, y1, r1 = dataset[0]
    x2, y2, r2 = dataset[1]
    x3, y3, r3 = dataset[2]

    # calculate system coefficients
    A = 2*x2 - 2*x1
    B = 2*y2 - 2*y1
    C = r1**2 - r2**2 - x1**2 + x2**2 - y1**2 + y2**2

    D = 2*x3 - 2*x2
    E = 2*y3 - 2*y2
    F = r2**2 - r3**2 - x2**2 + x3**2 - y2**2 + y3**2

    # calculate coordinates for system solution
    x = (C*E - F*B) / (E*A - B*D)
    y = (C*D - A*F) / (B*D - A*E)

    return x, y


SAMPLE_DATASET = [
    [
        -500, -200, 450
    ],
    [
        100, -100, 694.6221994724903
    ],
    [
        500, 100, 1011.1874208078342
    ]
]

coordinates = trilaterate(SAMPLE_DATASET)
print("Location:", coordinates)

Then i get the correct location (-500.0, 250.0).
I wish to know your opinion about this.

[mrJean1]

Function trilaterate3d2 uses a novel method to find the circle intersections (implemented in PyGeodesy with a Moore-Penrose pseudo-inverse, a null-space, numpy.linalg.svd and numpy.polynomial.polynomial.polyroots).

That does not seem to work well in this particular case. If r3 is rounded up to 1011.18742081 (8 decimals from 1011.1874208078342), then 2** intersections are found (Vector3d(-500.0, 250.0, -0.00382), Vector3d(-500.0, 250.0, 0.00382)). Other than that, I have no further explanation at this moment, but will look into this further. Perhaps adding some epsilon would help.

There are other several other trilaterate functions and LatLon.methods in PyGeodesy using different approaches, but those accept only lat- and longitudes.
__
**) Another issue is that trilaterate3d2 should not return 2 intersections with the same x and y values and opposite z-s if all given z-s are zero or the same.

[MartinManganiello]

Perfect, if I give a value of 1 to z then it returns the correct coordinates. Thank you very much

[mrJean1]

Good to know, but that shouldn't be.

In any case, the examples you provided will become new test cases once there is some resolution. Thank you for reporting the problems and preparing the examples.

[mrJean1]

Just FYI. Perturbing the radii passed to function trilaterate3d2 with up to 5 different epsilon values does produce the expected intersections. Those epsilon values** are 0, EPS, -EPS, EPSx4 and -EPSx4. An upcoming PyGeodesy release will include that.
__
**) The perturbation will only be applied if trilaterate3d2 argument eps is non-zero, which is the default.

[mrJean1]

PyGeodesy 21.04.04 is on GitHub and PyPi with the updated trilaterate3d2 function.

However, there are 2 minor pypy test failures for trilaterate3d2 to be addressed in the next version.

Please review and close this issue.

[mrJean1]

Just FYI, PyGeodesy 21.4.12 includes -among several other things- a new function called trilaterate2d2 trilaterating 3 circles in 2-D using the system equations you showed earlier plus optionally validating the trilaterated point.