SplineInterpolation.java

package heyingzhe;

import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;

public class SplineInterpolation {

    public static double[][] splineInterpolation(double[][] data) {
        int rows = data.length;
        int cols = data[0].length;

        // 检查是否整个二维数组全是0
        boolean allZero = true;
        for (double[] row : data) {
            for (double value : row) {
                if (value != 0) {
                    allZero = false;
                    break;
                }
            }
            if (!allZero) break;
        }
        if (allZero) {
            return data.clone(); // 如果没有有效数据，直接返回原始数组
        }

        // 创建结果数组
        double[][] result = new double[rows][cols];

        // 行插值
        for (int i = 0; i < rows; i++) {
            result[i] = interpolate1D(data[i]);
        }

        // 列插值
        for (int j = 0; j < cols; j++) {
            double[] column = new double[rows];
            for (int i = 0; i < rows; i++) {
                column[i] = result[i][j];
            }

            double[] interpolatedColumn = interpolate1D(column);
            for (int i = 0; i < rows; i++) {
                result[i][j] = interpolatedColumn[i];
            }
        }

        return result;
    }

    private static double[] interpolate1D(double[] data) {
        int n = data.length;

        // 收集有效数据点
        int count = 0;
        for (double v : data) {
            if (v != 0) count++;
        }

        // 如果没有有效数据点，或有效点少于 2 个，直接返回原数据
        if (count < 2) return data.clone();

        // 存储有效数据点的索引和值
        double[] x = new double[count];
        double[] y = new double[count];
        int idx = 0;
        for (int i = 0; i < n; i++) {
            if (data[i] != 0) {
                x[idx] = i;
                y[idx] = data[i];
                idx++;
            }
        }

        // 计算样条插值
        return cubicSplineInterpolate(x, y, n);
    }

    private static double[] cubicSplineInterpolate(double[] x, double[] y, int n) {
        // 如果输入数据点不足，直接返回全零数组
        if (x.length < 2 || y.length < 2) {
            return new double[n];
        }

        int m = x.length - 1; // 分段数量

        // Step 1: 计算每一段的间隔
        double[] h = new double[m];
        for (int i = 0; i < m; i++) {
            h[i] = x[i + 1] - x[i];
        }

        // Step 2: 构建三对角矩阵的右侧向量
        double[] alpha = new double[m];
        for (int i = 1; i < m; i++) {
            alpha[i] = (3.0 / h[i]) * (y[i + 1] - y[i]) - (3.0 / h[i - 1]) * (y[i] - y[i - 1]);
        }

        // Step 3: 构建三对角矩阵并求解曲率
        double[] l = new double[m + 1];
        double[] mu = new double[m];
        double[] z = new double[m + 1];
        l[0] = 1.0;
        z[0] = 0.0;
        mu[0] = 0.0;

        for (int i = 1; i < m; i++) {
            l[i] = 2.0 * (x[i + 1] - x[i - 1]) - h[i - 1] * mu[i - 1];
            mu[i] = h[i] / l[i];
            z[i] = (alpha[i] - h[i - 1] * z[i - 1]) / l[i];
        }

        l[m] = 1.0;
        z[m] = 0.0;

        double[] c = new double[m + 1];
        double[] b = new double[m];
        double[] d = new double[m];

        for (int j = m - 1; j >= 0; j--) {
            c[j] = z[j] - mu[j] * c[j + 1];
            b[j] = (y[j + 1] - y[j]) / h[j] - h[j] * (c[j + 1] + 2.0 * c[j]) / 3.0;
            d[j] = (c[j + 1] - c[j]) / (3.0 * h[j]);
        }

        // Step 4: 计算插值点的值
        double[] result = new double[n];
        for (int i = 0; i < n; i++) {
            result[i] = 0; // 初始化无效数据点为0
            for (int j = 0; j < m; j++) {
                if (x[j] <= i && i <= x[j + 1]) {
                    double diff = i - x[j];
                    result[i] = y[j] + b[j] * diff + c[j] * diff * diff + d[j] * diff * diff * diff;
                    break;
                }
            }
        }

        return result;
    }

    public static void main(String[] args) {
        double[][] data = {
                {0, 0, 0.5, 0, 0},
                {0, 1.2, 0, 0.8, 0},
                {0, 0, 0, 0, 0},
                {2.1, 0, 0, 0, 0.7}
        };

        double[][] interpolatedData = splineInterpolation(data);

        for (double[] row : interpolatedData) {
            for (double value : row) {
                System.out.printf("%.2f ", value);
            }
            System.out.println();
        }
    }
}