增益补偿算法实现后全景图出现暗带问题求助

我已实现论文《Automatic Panoramic Image Stitching using Invariant Features》中的Gain Compensation算法，推导过程无误（推导过程参考附图）。各图像的RGB增益值如下：

gainRGB =

0.6426    0.5983    0.6804
0.5293    0.5094    0.5329
0.4528    0.4353    0.4595
1.0402    0.9467    1.1177
0.6170    0.5957    0.6045
0.8882    0.8748    0.8587
1.1474    1.1198    1.2398
0.7990    0.7952    0.7872
1.3044    1.3110    1.2684
0.8528    0.8714    0.8841
1.6154    1.6276    1.7521
0.6150    0.5949    0.6002
0.8667    0.8438    0.8930
0.6304    0.6162    0.6599
1.0281    1.0020    1.0150
0.9782    0.9291    1.0337
0.6405    0.6108    0.6798
0.5141    0.4800    0.5499

应用上述增益补偿后，生成的全景图出现暗带（效果参考附图）。以下是复现该问题的最小化代码，恳请协助解决：

clc; clear;

% Initialze
load matlab.mat

n = length(warpedImages);
sigmaN = 10;
sigmag = 0.1;

Amat = cell(n);
Bvec = zeros(n,n);
Bvec_temp = zeros(n,1);
IuppeIdx = nonzeros(triu(reshape(1:numel(Amat), size(Amat))));
Amat_temp = cell(1,length(IuppeIdx));    
matSize = size(Amat);

% Get the Ibarijs and Nijs
tic
parfor i = 1:length(IuppeIdx)

    % Index to subscripts
    [ii,jj] = ind2sub(matSize, IuppeIdx(i));                 

    if ii == jj
        diag_val_1 = 0;
        diag_val_2 = 0;
       
        Z = 1:n;
        Z(Z==ii) = [];
        for d = Z
            [Ibarij, Ibarji, Nij] = getIbarNij(warpedImages{ii}, warpedImages{d});
            diag_val_1 = diag_val_1 + ((Nij * Ibarij.^2 + Nij * Ibarij.^2) / sigmaN^2);
            diag_val_2 = diag_val_2 + (Nij / sigmag^2);
        end
        
        Amat_temp{i} = diag_val_1 + diag_val_2;
        Bvec_temp(i) = diag_val_2;
    end       

    if ii ~= jj
        [Ibarij,Ibarji,Nij] = getIbarNij(warpedImages{ii}, warpedImages{jj});
        Amat_temp{i} = - ((Nij+Nij) * (Ibarij .* Ibarji)) /sigmaN^2 ;
    end

end
fprintf('Gain Compensation (par-for): %f\n', toc);

% --------------------------------------------------------------------------------------------------------------
% Form matrices
updiagIdx = nonzeros(triu(reshape(1:numel(Amat), size(Amat)),1));
lowdiagIdx = nonzeros(tril(reshape(1:numel(Amat), size(Amat)),-1));

% Populate A matrix
Amat(IuppeIdx) = Amat_temp;
Amat(lowdiagIdx) = Amat(updiagIdx);
Bvec(IuppeIdx) = Bvec_temp; 
Bvec = diag(Bvec);
Amat = cell2mat(Amat);    

% A matrix gain values
gainmatR = Amat(:,1:3:size(Amat,2));
gainmatG = Amat(:,2:3:size(Amat,2));
gainmatB = Amat(:,3:3:size(Amat,2));

% --------------------------------------------------------------------------------------------------------------
% AX = b --> X = A \ b
gR = gainmatR \ Bvec;
gG = gainmatG \ Bvec;
gB = gainmatB \ Bvec;

% Concatenate RGB gains
gainRGB = [gR, gG, gB];

gainRGB

% --------------------------------------------------------------------------------------------------------------
% Compensate gains for images        
gains = num2cell(gainRGB,2);
gains = cellfun(@(x) reshape(x, [1,1,3]), gains, 'UniformOutput',false);
gainImages = cellfun(@(x,y) uint8(double(x).*y), warpedImages, gains,'UniformOutput',false);
gainImagesMask = cellfun(@(x) repmat(imfill(imbinarize(rgb2gray(255 * x)), 'holes'), 1, 1, size(warpedImages{1},3)), ...
                        gainImages, 'UniformOutput',false);

% --------------------------------------------------------------------------------------------------------------
% Construct gain comepensated panorama    
gainpanorama = zeros(size(warpedImages{1}), 'uint8');
tic
for i = 1:n
    gainpanorama(gainImagesMask{i}) = gainImages{i}(gainImagesMask{i});
end
fprintf('Gain gain comepensated panorama (for loop): %f\n', toc);

figure; imshow(gainpanorama)


%--------------------------------------------------------------------------------------------------------
% Auxillary functions
%--------------------------------------------------------------------------------------------------------

%--------------------------------------------------------------------------------------------------------
% [Ibarij,Ibarji,Nij] = getIbarNij(Imij, Imji)
%--------------------------------------------------------------------------------------------------------
%
function [Ibarij,Ibarji,Nij] = getIbarNij(Imij, Imji)

% Overlay the warpedImage onto the panorama.
maski = imbinarize(rgb2gray(255 * Imij));  
maskj = imbinarize(rgb2gray(255 * Imji));

% Find the overlap mask
Nij_im  = maski & maskj;

% Get tje Ibarij and Nijs
if sum(sum(Nij_im)) ~= 0
    Nij_im  = imfill(Nij_im, 'holes');    
    Nijidx  = repmat(Nij_im, 1, 1, size(Imij,3));

    % Get the overlapping region RGB values for two images
    Ibarij = double(reshape(Imij(Nijidx),1,[],3));
    Ibarji = double(reshape(Imji(Nijidx),1,[],3));

    % Nij
    Nij = sum(sum(Nij_im));

    % Ibar ijs
    Ibarij = reshape(sum(Ibarij) ./ Nij, 1, 3);
    Ibarji = reshape(sum(Ibarji) ./ Nij, 1, 3);
else
    % Nij
    Nij = 0;

    % Ibar ijs
    Ibarij = zeros(1, 3);
    Ibarji = zeros(1, 3);
end


end

内容的提问来源于stack exchange，提问作者PManjunatha