function GeoCorrTIFF(data, rpcfilename, tiffname, useDEM)
%% 读取RPC文件
RPC = Rpc_Reader(rpcfilename);

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

%% 定义角点和位置矩阵
pixel_coords = [1,                  1,                  0;
    rows,               1,                  0;
    1,                  cols,               0;
    rows,               cols,               0;
    round(rows/2),      round(cols/2),      0;
    round(rows/4),      round(cols/4),      0;
    round(3*rows/4),    round(cols/4),      0;
    round(rows/4),      round(3*cols/4),    0;
    round(3*rows/4),    round(3*cols/4),    0;
    round(rows/4),      round(cols/2),      0;
    round(3*rows/4),    round(cols/2),      0;
    round(rows/2),      round(cols/4),      0;
    round(rows/2),      round(3*cols/4),    0;];

% 预分配内存空间
pixel_ll = zeros(13, 2);
srcPoints = pixel_coords(:, 1:2); % 直接提取前两列

%% 计算每个角点的值
for i = 1:13
    pixel_ll(i, :) = RpcBackwardCompute(RPC, pixel_coords(i, :));
    % pixel_ll(i, :) = RS_backward_compute_rpc(rpcfilename, pixel_coords(i, :)); % 速度慢
end

%% 使用DEM进行精校正
if(useDEM)
    % 以提取处的经纬度坐标，读取DEM影像，获取其高程值
    lat_range = [min(pixel_ll(:, 1)), max(pixel_ll(:, 1))];
    lon_range = [min(pixel_ll(:, 2)), max(pixel_ll(:, 2))];
    total_range = [lat_range, lon_range];

    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:13
        % 以下两个是阈值相关的参数
        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
%%
row_resolution = haversine(pixel_ll(1,1), pixel_ll(1,2), pixel_ll(2,1), pixel_ll(2,2)) * 1000 / rows;
col_resolution = haversine(pixel_ll(1,1), pixel_ll(1,2), pixel_ll(3,1), pixel_ll(3,2)) * 1000 / cols;
fprintf("影像原始分辨率：行向：%.2f(m); 列向：%.2f(m)。\n", row_resolution, col_resolution);
srcPoints(:, 1) = srcPoints(:, 1) * row_resolution / 30;
srcPoints(:, 2) = srcPoints(:, 2) * col_resolution / 30;

%%
[dstPoints, epsg_code] = LatLonToXY(pixel_ll); % 转换为投影坐标系的坐标(m)
% fprintf("投影坐标系：%d\n", epsg_code);
%% 设置 x 和 y 的世界坐标限制
xWorldLimits = [min(dstPoints(:, 2)), max(dstPoints(:, 2))];
yWorldLimits = [min(dstPoints(:, 1)), max(dstPoints(:, 1))];

% 地理坐标系转换及缩放
dstPoints(:, 1) = max(dstPoints(:, 1)) - dstPoints(:, 1);
dstPoints(:, 2) = dstPoints(:, 2) - min(dstPoints(:, 2));
dstPoints = dstPoints / 30; % 分辨率为 30m

%% 重投影
tform = fitgeotrans(srcPoints, dstPoints, 'projective');
% 重采样并创建栅格信息
outputImage = imwarp(data, tform, 'InterpolationMethod', 'bilinear');
R = maprefcells(xWorldLimits, yWorldLimits, size(outputImage));

%% 导出为投影坐标系
% 创建 GeoTIFF 标签
GeoKeyDirectoryTag = struct('GTModelTypeGeoKey', 1, ...
    'GTRasterTypeGeoKey', 1, ...
    'ProjectedCSTypeGeoKey', epsg_code);

% 保存 GeoTIFF 文件
geotiffwrite(tiffname, outputImage, R, ...
    'CoordRefSysCode', epsg_code, ...
    'GeoKeyDirectoryTag', GeoKeyDirectoryTag);

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

end
