TrendSurfaceInterpolation.java

package heyingzhe;

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

public class TrendSurfaceInterpolation {

    /**
     * 根据趋势面插值填充二维数组。
     *
     * @param inputArray 包含观测值和0的二维数组
     * @return 插值后的二维数组
     */
    public static double[][] trendSurfaceInterpolation(double[][] inputArray) {
        int height = inputArray.length;
        int width = inputArray[0].length;

        // 收集所有非零点作为观测值
        List<Point> points = new ArrayList<>();
        for (int i = 0; i < height; i++) {
            for (int j = 0; j < width; j++) {
                if (inputArray[i][j] != 0) {
                    points.add(new Point(j, i, inputArray[i][j]));
                }
            }
        }

        // 如果没有足够的点进行拟合，直接返回
        if (points.size() < 3) {
            throw new IllegalArgumentException("至少需要三个非零点进行插值计算！");
        }

        // 根据点拟合趋势面: z = a + b*x + c*y
        double[] coefficients = fitTrendSurface(points);

        // 使用拟合模型填充所有零点
        double[][] resultArray = new double[height][width];
        for (int i = 0; i < height; i++) {
            for (int j = 0; j < width; j++) {
                if (inputArray[i][j] == 0) {
                    resultArray[i][j] = coefficients[0] + coefficients[1] * j + coefficients[2] * i;
                } else {
                    resultArray[i][j] = inputArray[i][j];
                }
            }
        }

        return resultArray;
    }

    /**
     * 拟合趋势面 z = a + b*x + c*y。
     *
     * @param points 包含所有观测值的点集合
     * @return 拟合的趋势面系数 [a, b, c]
     */
    private static double[] fitTrendSurface(List<Point> points) {
        int n = points.size();
        double sumX = 0, sumY = 0, sumZ = 0;
        double sumX2 = 0, sumY2 = 0, sumXY = 0, sumXZ = 0, sumYZ = 0;

        for (Point p : points) {
            double x = p.x;
            double y = p.y;
            double z = p.z;

            sumX += x;
            sumY += y;
            sumZ += z;
            sumX2 += x * x;
            sumY2 += y * y;
            sumXY += x * y;
            sumXZ += x * z;
            sumYZ += y * z;
        }

        // 构造矩阵和向量
        double[][] A = {
                {n, sumX, sumY},
                {sumX, sumX2, sumXY},
                {sumY, sumXY, sumY2}
        };
        double[] B = {sumZ, sumXZ, sumYZ};

        // 解线性方程组 A * [a, b, c] = B
        return solveLinearSystem(A, B);
    }

    /**
     * 解线性方程组 A * X = B。
     *
     * @param A 系数矩阵
     * @param B 常数向量
     * @return 解向量 X
     */
    private static double[] solveLinearSystem(double[][] A, double[] B) {
        int n = B.length;
        double[] X = new double[n];

        // 高斯消元法
        for (int i = 0; i < n; i++) {
            // 寻找主元
            int max = i;
            for (int j = i + 1; j < n; j++) {
                if (Math.abs(A[j][i]) > Math.abs(A[max][i])) {
                    max = j;
                }
            }

            // 交换行
            double[] temp = A[i];
            A[i] = A[max];
            A[max] = temp;

            double t = B[i];
            B[i] = B[max];
            B[max] = t;

            // 归一化当前行
            for (int j = i + 1; j < n; j++) {
                double factor = A[j][i] / A[i][i];
                B[j] -= factor * B[i];
                for (int k = i; k < n; k++) {
                    A[j][k] -= factor * A[i][k];
                }
            }
        }

        // 回代求解
        for (int i = n - 1; i >= 0; i--) {
            double sum = 0;
            for (int j = i + 1; j < n; j++) {
                sum += A[i][j] * X[j];
            }
            X[i] = (B[i] - sum) / A[i][i];
        }

        return X;
    }

    /**
     * 点类，存储 x, y, z 坐标。
     */
    private static class Point {
        double x, y, z;

        Point(double x, double y, double z) {
            this.x = x;
            this.y = y;
            this.z = z;
        }
    }

    public static void main(String[] args) {
        // 示例输入数组
        double[][] inputArray = {
                {0, 0, 0, 0, 0},
                {0, 10, 0, 20, 0},
                {0, 0, 0, 0, 0},
                {0, 30, 0, 40, 0},
                {0, 0, 0, 0, 0}
        };

        // 调用插值函数
        double[][] resultArray = trendSurfaceInterpolation(inputArray);

        // 输出结果
        for (double[] row : resultArray) {
            for (double val : row) {
                System.out.printf("%.2f ", val);
            }
            System.out.println();
        }
    }
}