function GeoCorrTIFF4(data, rpcfilename, tiffname, varargin)
if(nargin==3)
    useDEM = false;
elseif(nargin==4)
    useDEM = varargin{1};
end
%% 读取RPC文件
RPC = Rpc_Reader(rpcfilename);

%% 获取影像的行列数并翻转数据
[rows, cols, wavs] = size(data);

%% 定义角点和位置矩阵
nn = 4;
row_ori = linspace(1, rows, nn);
col_ori = linspace(1, cols, nn);
[mycol, myrow] = meshgrid(col_ori, row_ori);

pixel_coords = [myrow(:), mycol(:), zeros(nn^2, 1)];

points_num = size(pixel_coords, 1);

% 预分配内存空间
pixel_ll = zeros(points_num, 2);

%% 计算每个角点的值
for i = 1 : points_num
    pixel_ll(i, :) = RpcBackwardCompute(RPC, pixel_coords(i, :));
end
pixel_lat = pixel_ll(:, 1);
pixel_lon = pixel_ll(:, 2);

% 以提取处的经纬度坐标，读取DEM影像，获取其高程值
lat_range = [min(pixel_lat), max(pixel_lat)];
lon_range = [min(pixel_lon), max(pixel_lon)];

%% 使用DEM进行精校正
% 精校正需要SRTM数据，下载时开全局代理
% url = 'https://e4ftl01.cr.usgs.gov/MEASURES/SRTMGL1.003/2000.02.11';
total_range = [lat_range, lon_range];
if(useDEM)
    key = input('[请输入NASA账户密码]\nUsename  :whu_why2021\nLogin Key:', 's');
    DEM = readhgt(total_range, 'srtm1',...
        'login', 'whu_why2021', key, 'interp',...
        'outdir', 'E:/SRTM_data'); % 'plot' 用于绘制

    for i = 1 : points_num
        % 以下两个是阈值相关的参数
        old_ll = [0, 0];
        end_yuzhi = 100;
        while( end_yuzhi > 1.5e-4)
            value_bijin=[pixel_ll(i, 1),pixel_ll(i, 2)];
            pixel_ll(i, :) = RpcBackwardCompute(RPC, pixel_coords(i, :), value_bijin);
            dem_lat = DEM.lat;
            dem_lon = DEM.lon;
            dem_z = DEM.z;
            [~, row] = min(abs(pixel_ll(i, 1) - dem_lat));
            [~, col] = min(abs(pixel_ll(i, 2) - dem_lon));
            value_dem = double(dem_z(row, col));
            pixel_coords(i, 3) = value_dem;
            end_yuzhi = sum(abs(pixel_ll - old_ll));
            old_ll = pixel_ll;
        end
    end
end
pixel_lat = pixel_ll(:, 1);
pixel_lon = pixel_ll(:, 2);
%% 获取图像上每个点的经纬度坐标
[COLS, ROWS] = meshgrid(1:cols, 1:rows);
% 使用 griddata 进行插值 (这里选择 'linear' 方法)
pixel_lat_all = griddata(mycol(:), myrow(:), pixel_lat, COLS, ROWS, 'linear');
pixel_lon_all = griddata(mycol(:), myrow(:), pixel_lon, COLS, ROWS, 'linear');

%% 设置输出分辨率
res = 3e-4; % 约为30m
[lonq, latq] = meshgrid(lon_range(1):res:lon_range(2), lat_range(1):res:lat_range(2));
% 使用 griddata 进行插值 (这里选择 'linear' 方法)
data_ori = reshape(data, rows*cols, wavs);
grid_data = zeros([size(lonq), wavs]);
% 并行计算
parfor j = 1:wavs
    grid_data(:,:,j) = griddata(pixel_lon_all(:), pixel_lat_all(:), double(data_ori(:, j)), lonq, latq, 'linear');
end
data_new = cast(grid_data, class(data));

%% 导出为地理坐标系 EPSG：4326
% 创建新的空间参考，单位为地理坐标（经纬度）
R_geo = georefcells(lat_range, lon_range, size(data_new));
% 创建 GeoTIFF 标签
GeoKeyDirectoryTag = struct('GTModelTypeGeoKey', 2, ... % 地理坐标系类型
    'GTRasterTypeGeoKey', 1); % 栅格像素的类型
geotiffwrite(tiffname, data_new, R_geo, 'CoordRefSysCode', 4326,'GeoKeyDirectoryTag', GeoKeyDirectoryTag);

end