function GeoCorrTIFF3(data, rpcfilename, tiffname, varargin)
if(nargin==3)
    useDEM = false;
elseif(nargin==4)
    useDEM = varargin{1};
end
%% 读取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;];

points_num = size(pixel_coords, 1);

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

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

% 以提取处的经纬度坐标，读取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];

%% 使用DEM进行精校正
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
%% 计算像元坐标
param_lat = 111.194926644558; % 一纬度的距离，固定(km)
param_lon = haversine(mean(lat_range), lon_range(1), mean(lat_range), lon_range(2))/...
    (lon_range(2) - lon_range(1));

row_resolution1 = haversine(pixel_ll(1,1), pixel_ll(1,2), pixel_ll(2,1), pixel_ll(2,2)) * 1000 / rows;
row_resolution2 = haversine(pixel_ll(4,1), pixel_ll(4,2), pixel_ll(3,1), pixel_ll(3,2)) * 1000 / rows;
col_resolution1 = haversine(pixel_ll(1,1), pixel_ll(1,2), pixel_ll(3,1), pixel_ll(3,2)) * 1000 / cols;
col_resolution2 = haversine(pixel_ll(4,1), pixel_ll(4,2), pixel_ll(2,1), pixel_ll(2,2)) * 1000 / cols;
row_resolution = (row_resolution1 + row_resolution2)/2;
col_resolution = (col_resolution1 + col_resolution2)/2;

fprintf("影像原始分辨率：行向：%.2f(m); 列向：%.2f(m)。\n", row_resolution, col_resolution);

%% 插值统一分辨率
rows_new = round(rows * row_resolution / min([row_resolution col_resolution]));
cols_new = round(cols * col_resolution / min([row_resolution col_resolution]));
[COLS, ROWS] = meshgrid(linspace(1, cols_new, cols), linspace(1, rows_new, rows));
[COLS_NEW, ROWS_NEW] = meshgrid(1:cols_new, 1:rows_new);
data_new = zeros(rows_new, cols_new, size(data,3));
for j = 1:size(data,3)
    data_new(:,:,j) = interp2(COLS, ROWS, double(data(:,:,j)), COLS_NEW, ROWS_NEW);
end
data_new = cast(data_new, class(data));
srcPoints2 = srcPoints.* [row_resolution, col_resolution] / min([row_resolution, col_resolution]);
%% 生成目标点 （单位：格）
weight1 = param_lat * 1e3 / 30; % 单位纬度的格网距离
weight2 = param_lon * 1e3 / 30; % 单位经度的格网距离

dstPoints(:, 1) = round((max(pixel_ll(:, 1)) - pixel_ll(:, 1)) .* weight1) + 1;
dstPoints(:, 2) = round((pixel_ll(:, 2) - min(pixel_ll(:, 2))) .* weight2) + 1;
%%
% 定义填充大小
padSize = 100; % 根据你的权重调整padSize
% 使用黑色填充扩展图像
data_padded = padarray(data_new, [padSize, padSize], 0, 'both'); % 使用0作为填充值（即黑色）
% 调整源点坐标以适应新的填充图像
srcPoints2_padded = srcPoints2 + padSize;
%% 重投影
% 计算变换矩阵
tform = fitgeotrans(srcPoints2_padded, dstPoints, 'projective');
% 应用几何变换
outputImage = imwarp(data_padded, tform, 'InterpolationMethod', 'nearest');
%% 移除纯黑边框
% 自动检测并移除纯黑边框
if size(outputImage, 3) == 3
    imgGray = rgb2gray(outputImage);
else
    imgGray = outputImage;
end

% 检测非黑行和列
rowSums = sum(imgGray, 2);
colSums = sum(imgGray, 1);

% 找到非黑边缘
rowIdx = find(rowSums > 0, 1, 'first'):find(rowSums > 0, 1, 'last');
colIdx = find(colSums > 0, 1, 'first'):find(colSums > 0, 1, 'last');

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


