ClosestPointDistanceToClosestPointOnSegment fails for some pairs on globe... see video and code.. any ideas ?

See picture in second comment for problem.

See here for example pairs fail, near end of video:

https://www.youtube.com/live/-EsixpO9mBk?feature=share

// engineer who provided it to you for assistance. // Define some earth constants const MAJA = (6378137.0); // semi major axis of ref ellipsoid const FLAT = (1.0/298.2572235630); // flattening coef of ref ellipsoid. // These are derived values from the above WGS-84 values const ESQR = (FLAT (2.0-FLAT)); // 1st eccentricity squared const OMES = (1.0 - ESQR); // 1 minus eccentricity squared const EFOR = (ESQR ESQR); // Sqr of the 1st eccentricity squared const ASQR = (MAJA * MAJA); // semi major axis squared := nlmaja**

procedure ConvertECEFToLTP ( const ParaEcef : TEcef; var ParaLtp : TLtp ); var a0,a1,a2,a3,a4,b0,b1,b2,b3 : double; b,c0,opqk,qk,qkc,qks,f,fprm : double; k : integer; ytemp, xtemp : double; begin // Convert from ParaEcef (XYZ) to LTP (ParaLtp.Lat/ParaLtp.Lon/ParaLtp.Alt)

//-------------------------------------------
// b := (0*0 + 1*1) / semi_major_axis_squared
// c := (z*z) / semi_major_axis_squared
//-------------------------------------------
b :=
    (
        (ParaEcef.x * ParaEcef.x) +
        (ParaEcef.y * ParaEcef.y)
    ) / ASQR;
c0 :=
    (ParaEcef.z * ParaEcef.z) / ASQR;

//-------------------------------------------
// a0 := c * one_minus_eccentricity_sqr
// a1 := 2 * a0
// a2 := a0 + b - first_eccentricity_to_fourth
// a3 := -2.0 * first_eccentricity_to_fourth
// a4 := - first_eccentricity_to_fourth
//-------------------------------------------
a0 := OMES * c0;
a1 := 2.0 * a0;
a2 := a0 + b - EFOR;
a3 := -2.0 * EFOR;
a4 := - EFOR;

//-------------------------------------------
// b0 := 4 * a0, b1 := 3 * a1
// b2 := 2 * a2, b3 := a3
//-------------------------------------------
b0 := 4.0 * a0;
b1 := 3.0 * a1;
b2 := 2.0 * a2;
b3 := a3;

//-------------------------------------------
// Compute First Eccentricity Squared
//-------------------------------------------
qk := ESQR;

// for (k := 1; k <:= 3; k++) for k := 1 to 3 do begin qks := qk qk; qkc := qk qks; f := (a0 qks qks) + (a1 qkc) + (a2 qks) + (a3 qk) + a4; fprm := (b0 qkc) + (b1 qks) + (b2 qk) + b3; qk := qk - (f / fprm); end;

//-------------------------------------------
// Compute latitude, longitude, altitude
//-------------------------------------------
opqk := 1.0 + qk;
if
(
    (ParaEcef.x = 0.0) and
    (ParaEcef.y = 0.0)
) then
// on the earth's axis
begin
    // We are sitting EXACTLY on the earth's axis
    // Probably at the center or on or near one of the poles
    ParaLtp.Lon := 0.0; // as good as any other value
    if (ParaEcef.z >= 0.0) then
    begin
        ParaLtp.Lat := PI / 2; // ParaLtp.Alt above north pole
    end else
    begin
        ParaLtp.Lat := -PI / 2; // ParaLtp.Alt above south pole
    end;
end else
begin
    ytemp := opqk * ParaEcef.z;
    xtemp := sqrt
    (
        ( ParaEcef.x * ParaEcef.x ) +
        ( ParaEcef.y * ParaEcef.y )
    );
    ParaLtp.Lat := arctan2( ytemp, xtemp );
    ParaLtp.Lon := arctan2( ParaEcef.y, ParaEcef.x );
end;
ParaLtp.Alt :=
    (
        1.0 - OMES / ESQR * qk
    ) * MAJA *
    sqrt
    (
        (
            b / (opqk * opqk)
        ) + c0
    );

end; // ConvertECEFToLTP()

// written by Skybuck Flying ! ;) :) =D function ClosestPointDistanceToClosestPointOnSegment ( ParaPointX, ParaPointY, ParaPointZ : double; ParaStartX, ParaStartY, ParaStartZ : double; ParaEndX, ParaEndY, ParaEndZ : double; ParaRadius : double {= 6371e3} // (Mean) radius of earth (defaults to radius in metres). ) : double; var vPoint, vSegmentPoint1, vSegmentPoint2, vClosestPointOnSegment : TLatLonNvectorSpherical; begin vSegmentPoint1.FromXYZ2( ParaStartX, ParaStartY, ParaStartZ, ParaRadius ); vSegmentPoint2.FromXYZ2( ParaEndX, ParaEndY, ParaEndZ, ParaRadius ); vPoint.FromXYZ2( ParaPointX, ParaPointY, ParaPointZ, ParaRadius );

vClosestPointOnSegment := nearestPointOnSegment
(
    vPoint,
    vSegmentPoint1,
    vSegmentPoint2,
    ParaRadius
);

result := distanceTo
(
    vPoint, vClosestPointOnSegment, ParaRadius
);

end;

function ClosestPointDistanceToClosestPointOnSegmentV2 ( ParaPointX, ParaPointY, ParaPointZ : double; ParaStartX, ParaStartY, ParaStartZ : double; ParaEndX, ParaEndY, ParaEndZ : double; ParaRadius : double {= 6371e3} // (Mean) radius of earth (defaults to radius in metres). ) : double; var vPoint, vSegmentPoint1, vSegmentPoint2, vClosestPointOnSegment : TLatLonNvectorSpherical; vEcef : Tecef; vLtp : Tltp; begin // convert start x,y,z to lat lng vEcef.X := ParaStartX; vEcef.Y := ParaStartY; vEcef.Z := ParaStartZ;

ConvertECEFToLTP
(
    vEcef,
    vLtp
);

vSegmentPoint1.mLat := vLtp.Lat;
vSegmentPoint1.mLon := vLtp.Lon;

// convert end x,y,z to lat lang
vEcef.X := ParaEndX;
vEcef.Y := ParaEndY;
vEcef.Z := ParaEndZ;

ConvertECEFToLTP
(
    vEcef,
    vLtp
);

vSegmentPoint2.mLat := vLtp.Lat;
vSegmentPoint2.mLon := vLtp.Lon;

// ParaRadius ??

// convert point x,y,z to lat lng
vEcef.X := ParaPointX;
vEcef.Y := ParaPointY;
vEcef.Z := ParaPointZ;

ConvertECEFToLTP
(
    vEcef,
    vLtp
);

vPoint.mLat := vLtp.Lat;
vPoint.mLon := vLtp.Lon;

vClosestPointOnSegment := nearestPointOnSegment
(
    vPoint,
    vSegmentPoint1,
    vSegmentPoint2,
    ParaRadius
);

result := distanceTo
(
    vPoint, vClosestPointOnSegment, ParaRadius
);

end;

Any idea what could be wrong ? If you need more code let me know...

chrisveness / geodesy

ClosestPointDistanceToClosestPointOnSegment fails for some pairs on globe... see video and code.. any ideas ? #111