How to interpolate ECEF coordinates on an WGS84 ellipsoid

Is there a direct method (not involving converting the coordinates to lat/lon) to interpolate between 2 ECEF coordinates (xyz) in order for the interpolated point to be located on the WGS84 ellispoid. The original 2 points are computed from geodetic coordinates.

Interpolating on a sphere seem obvious but I can't seem to derive a solution for the ellipsoid.

Thank you in advance.

Let assume you got 2 points p0(x,y,z) and p1(x,y,z) and want to interpolate some p(t) where t=<0.0,1.0> between the two.

you can:

1. rescale your ellipsoid to sphere

   simply like this:

   const double mz=6378137.00000/6356752.31414; // [m] equatoreal/polar radius of Earth
   p0.z*=mz;
   p1.z*=mz;
   

   now you got Cartesian coordinates refering to spherical Earth model.

2. interpolate

   simple linear interpolation would do

   p(t) = p0+(p1-p0)*t
   

   but of coarse you also need to normalize to earth curvature so:

   r0 = |p0|
   r1 = |p1|
   
   p(t) = p0+(p1-p0)*t
   r(t) = r0+(r1-r0)*t
   
   p(t)*=r/|p(t)|
   

   where |p0| means length of vector p0.

3. rescale back to ellipsoid

   by dividing with the same value

   p(t).z/=mz
   

This is simple and cheap but the interpolated path will not have linear time scale.

Here C++ example:

void XYZ_interpolate(double *pt,double *p0,double *p1,double t)
    {
    const double  mz=6378137.00000/6356752.31414;
    const double _mz=6356752.31414/6378137.00000;
    double p[3],r,r0,r1;
    // compute spherical radiuses of input points
    r0=sqrt((p0[0]*p0[0])+(p0[1]*p0[1])+(p0[2]*p0[2]*mz*mz));
    r1=sqrt((p1[0]*p1[0])+(p1[1]*p1[1])+(p1[2]*p1[2]*mz*mz));
    // linear interpolation
    r   = r0   +(r1   -r0   )*t;
    p[0]= p0[0]+(p1[0]-p0[0])*t;
    p[1]= p0[1]+(p1[1]-p0[1])*t;
    p[2]=(p0[2]+(p1[2]-p0[2])*t)*mz;
    // correct radius and rescale back
    r/=sqrt((p[0]*p[0])+(p[1]*p[1])+(p[2]*p[2]));
    pt[0]=p[0]*r;
    pt[1]=p[1]*r;
    pt[2]=p[2]*r*_mz;
    }

And preview:

Yellow squares are the used p0,p1 Cartesian coordinates, the White curve is the interpolated path where t=<0.0,1.0> ...