I have implemented the Gain Compensation algorithm as mentioned in the paper Automatic Panoramic Image Stitching using Invariant Features. The derivation looks fine, here is the attached one for reference:

The gains values for RGB are for each images:

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

When I compensate for the above gains, I get a dark-banded panoramic image as shown below:

Below is the minimal code to reproduce the output image. Here is the link to the image files. Any help to fix this is appreciated!

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): %fn', 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): %fn', 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