pagrangeInterpolation for curve fitting - 使用拉格朗日插值平滑曲线

使用拉格朗日插值平滑曲线

1 平滑前后效果对比：

平滑前

平滑后

2 拉格朗日插值代码:(参考cesium源码)

结合拉格朗日原理：https://zhuanlan.zhihu.com/p/135229305

源码如下：https://github.com/CesiumGS/cesium/blob/b30/Source/Core/LagrangePolynomialApproximation.js#L39

/**
 * Interpolates values using Lagrange Polynomial Approximation.
 *
 * @param {Number} x The independent variable for which the dependent variables will be interpolated.
 * @param {Number[]} xTable The array of independent variables to use to interpolate.  The values
 * in this array must be in increasing order and the same value must not occur twice in the array.
 * @param {Number[]} yTable The array of dependent variables to use to interpolate.  For a set of three
 * dependent values (p,q,w) at time 1 and time 2 this should be as follows: {p1, q1, w1, p2, q2, w2}.
 * @param {Number} yStride The number of dependent variable values in yTable corresponding to
 * each independent variable value in xTable.
 * @param {Number[]} [result] An existing array into which to store the result.
 * @returns {Number[]} The array of interpolated values, or the result parameter if one was provided.
 */
LagrangePolynomialApproximation.interpolateOrderZero = function(x, xTable, yTable, yStride, result) {
    if (!defined(result)) {
        result = new Array(yStride);
    }

    var i;
    var j;
    var length = xTable.length;

    for (i = 0; i < yStride; i++) {
        result[i] = 0;
    }

    for (i = 0; i < length; i++) {
        var coefficient = 1;

        for (j = 0; j < length; j++) {
            if (j !== i) {
                var diffX = xTable[i] - xTable[j];
                coefficient *= (x - xTable[j]) / diffX;
            }
        }

        for (j = 0; j < yStride; j++) {
            result[j] += coefficient * yTable[i * yStride + j];
        }
    }

    return result;
};

3 cesium中实现代码

放入https://sandcastle.cesium.com/index.html?src=Callback%20Property.html运行即可

// This example illustrates a Callback Property, a property whose
// value is lazily evaluated by a callback function.
// Use a CallbackProperty when your data can't be pre-computed
// or needs to be derived from other properties at runtime.
var viewer = new Cesium.Viewer("cesiumContainer");
viewer.clock.shouldAnimate = true;

var startLatitude = 35;
var startLongitude = -120;
var endLongitude;
var startTime = Cesium.JulianDate.now();
var stTime = startTime;


// Add a polyline to the scene. Positions are dynamic.
var isConstant = false;
var redLine = viewer.entities.add({
  polyline: {
    // This callback updates positions each frame.
    positions: new Cesium.CallbackProperty(function (time, result) {
      endLongitude =
        startLongitude +
        0.001 * Cesium.JulianDate.secondsDifference(time, startTime);
      
      // 平滑前效果
      // return computeCirclularFlight(-112.110693, 36.0994841, 0.03, time);
      // 平滑后效果
      return computeCirclularFlightSmoothly(-112.110693, 36.0994841, 0.03, 500, stTime, time);
      /*
      return Cesium.Cartesian3.fromDegreesArray(
        [startLongitude, startLatitude, endLongitude, startLatitude],
        Cesium.Ellipsoid.WGS84,
        result
      );
      */
    }, isConstant),
    width: 5,
    material: Cesium.Color.RED,
  },
});


//Generate a random circular pattern with varying heights.
function computeCirclularFlight(lon, lat, radius, curTime) {

        // let a = JulianDate.secondsDifference(JulianDate.now(), stTime);
         let positions = [];
        // Generate a random circular pattern with letying heights.
        let property = new Cesium.SampledPositionProperty();
        for (let i = 0; i <= 360; i += 45) {
            let radians = Cesium.Math.toRadians(i);
            let time = Cesium.JulianDate.addSeconds(
                stTime,
                i,
                new Cesium.JulianDate()
            );
        
            let position = Cesium.Cartesian3.fromDegrees(
                lon + radius * 1.5 * Math.cos(radians),
                lat + radius * Math.sin(radians),
                Cesium.Math.nextRandomNumber() * 500 + 
                    Math.random() * 20 * Cesium.JulianDate.secondsDifference(curTime, stTime)
            );
        
            positions.push(position);
            property.addSample(time, position);
      }

      return positions;
}

    // Generate a random circular pattern with varying heights.
    function computeCirclularFlightSmoothly(lon, lat, radius, levelHeight, stTime, curTime) {

        // let a = JulianDate.secondsDifference(JulianDate.now(), stTime);
        let positions = [];
        let llhPositions = [];
        // let property = new PositionProperty();
        // Generate a random circular pattern with letying heights.
        for (let i = 0; i <= 360; i += 45) {
            let radians = Cesium.Math.toRadians(i);
            let time = Cesium.JulianDate.addSeconds(
                stTime,
                i,
                new Cesium.JulianDate()
            );

            let tmpPoint = new Cesium.Cartographic(
                lon + radius * 1.5 * Math.cos(radians),
                lat + radius * Math.sin(radians),
                Cesium.Math.nextRandomNumber() * 0.1 * levelHeight + levelHeight +
                    Cesium.Math.random() * 20 * Cesium.JulianDate.secondsDifference(curTime, stTime)
            );
            llhPositions.push(tmpPoint);
            positions.push(
                Cesium.Cartesian3.fromDegrees(
                    tmpPoint.longitude,
                    tmpPoint.latitude,
                    tmpPoint.height
            ));
        }
 
        // let iptPositions = Array<Cartographic>();
        let iptPositions = [];
        let xRes = [];
        let xTable = [];        // interpolate at xTable[i]
        for(let ix = 0; ix < positions.length; ++ix) {
            xTable.push(ix * 100);
        }
        let yTable = [];        // used to interpolate, in {la1, lon1, h1, l2, lon2, h2}
        for(let iy = 0; iy < positions.length; ++iy) {
            yTable.push(llhPositions[iy].longitude);
            yTable.push(llhPositions[iy].latitude);
            yTable.push(llhPositions[iy].height);
        }
        let yStride = 3;        // 3 dependent vaule in yTable is viewed as one item
        
        for (let ix = 0; ix < xTable[positions.length - 1]; ++ix) {
            let iptRes = [];
            Cesium.LagrangePolynomialApproximation.interpolateOrderZero(
                ix,
                xTable,
                yTable,
                yStride,
                iptRes
            );

            iptPositions.push(
                Cesium.Cartesian3.fromDegrees(
                iptRes[0],
                iptRes[1],
                iptRes[2]
            ));
        }

        return  iptPositions;
    }

// Keep the view centered.
viewer.trackedEntity = redLine;