基于双线性内插的DEM高程点内插实现

针对tif格式DEM数据，由于其数据组成为等间距横向和纵向分布的离散高程点构成，可以看做为普通影像的像素值，点与点之间间距为一个像素，那么对DEM数据的内插，就可以当做在每个像素内均匀增加点，每个点的高程值则由周围四个点通过双线性内插计算结果来确定。

假设当点P(x,y)位于P1(x1,y1)、P2(x2,y2)、P3(x3,y3)和P4(x4,y4)四个点构成的矩阵内，如图所示：

设定四个点高程值为p1、p2、p3、p4，中间P点高程为p，P点高程通过P1-P4四个点高程进行拟合。这里采用距离加权的方法对内插点高程进行赋值，即通过计算P点坐标与其他几个点距离来确定每个值所占权重，假设P点刚好位于矩阵中心，那么p刚好等于p1-p4的均值。

问题的核心在于如何计算P点与周围四个点距离，由于P点在矩阵内，假设一个矩阵边长为1，x和y与周围点坐标的差值都小于等于1，那么可以通过如下公式来计算p点高程值。

通过以上公式，当P点位于矩阵中心时，x-x1=0.5,y-y1=0.5,P1点权重为0.5*0.5=0.25，P2到P4点权重也一样是0.25，如果P点刚好位于P1-P4点任意一点上，通过公式计算结果也是p1-p4。

核心代码：

double result = p1 * ((1 - Math.abs(x - x1)) * (1 - Math.abs(y - y1)))
				+ p2 * ((1 - Math.abs(x - x2)) * (1 - Math.abs(y - y2)))
				+ p3 * ((1 - Math.abs(x - x3)) * (1 - Math.abs(y - y3))) 
				+ p4 * ((1 - Math.abs(x - x4)) * (1 - Math.abs(y - y4)));

假设P1 = 100.25, p2= 120.23, p3 = 103.49, p4 = 105.66，P点坐标为(0.2,0.3)，计算结果p=112.217，四个点均值为107.4075，P点坐标最接近P2点，P2点最高，因此P点高程大于均值也更符合预期。

具体实现：

public void interpolation(File inputDEM,File outputDEM,int interpolation) {
		DemData data = new DemData();
		try {
			gdal.AllRegister();
			Dataset dataset = gdal.Open(inputDEM.toString(), gdalconstConstants.GA_ReadOnly);
			Band band = dataset.GetRasterBand(1);
			int width = dataset.GetRasterXSize();
			int height = dataset.GetRasterYSize();
			data.setWidth(width);
			data.setHeight(height);		
			double[] values = new double[width * height];
			double[] buf = new double[width];
			for (int i = 0; i < height; i++) {
				// 一行一行读取
				band.ReadRaster(0, i, width, 1, buf);
				for (int j = 0; j < width; j++) {
					int index = i * width + j;
					values[index] = buf[j];
				}
			}
			data.setValues(values);
			
			Driver driver = gdal.GetDriverByName("GTiff");
			int destWidth = (width - 1) * interpolation;
			int destHeight = (height - 1) * interpolation;
			Dataset demDataset = driver.Create(outputDEM.toString(), destWidth, destHeight, 1);
			Band demBand = demDataset.GetRasterBand(1);
			demBand.SetMetadata(band.GetMetadata_Dict());
			double[] resultLine = new double[destWidth];
			for (int i = 0; i < destHeight; i++) {
				for (int j = 0; j < destWidth; j++) {
					double x = (double)j / interpolation;
					double y = (double)i / interpolation;
					int x1 = (int)x;
					int x2 = x1;
					int x3 = x1 + 1;
					int x4 = x3;
					int y1 = (int)y + 1;
					int y2 = y1 - 1;
					int y3 = y2;
					int y4 = y1;
					double p1 = data.getValue(y1, x1);
					double p2 = data.getValue(y2, x2);
					double p3 = data.getValue(y3, x3);
					double p4 = data.getValue(y4, x4);
					double result = caculateElevation(x, y, x1, x2, x3, x4, y1, y2, y3, y4, p1, p2, p3, p4);
					resultLine[j] = result;
				}
				demBand.WriteRaster(0, i, destWidth, 1, resultLine);
			}
			demDataset.SetMetadata(dataset.GetMetadata_Dict());
			dataset.delete();
			demDataset.delete();
			} catch (Exception e) {
			e.printStackTrace();
		}
	}
	
	public double caculateElevation(double x, double y, int x1, int x2, int x3, int x4, int y1, int y2, int y3, int y4,
			double p1, double p2, double p3, double p4) {
		double result = p1 * ((1 - Math.abs(x - x1)) * (1 - Math.abs(y - y1)))
				+ p2 * ((1 - Math.abs(x - x2)) * (1 - Math.abs(y - y2)))
				+ p3 * ((1 - Math.abs(x - x3)) * (1 - Math.abs(y - y3))) 
				+ p4 * ((1 - Math.abs(x - x4)) * (1 - Math.abs(y - y4)));
		return result;
	}