solver3a

function Aest = solver3(imageSize, queryLimit)

submasksizevector = [2 4 6 8 10 12 14 16 Inf];
submasksize = sqrt(imageSize^2 / queryLimit);
index = find(submasksize < submasksizevector, 1, 'first');
submasksize = submasksizevector(index);
submasksize = min(submasksize, 16);

[submask, submasksum] = getsubmask(submasksize);
one_by_integer_lookup_table = 1.0 ./ (1 : submasksum);

nrOfSubmasks = ceil(imageSize / submasksize);
mask = zeros(imageSize);
Aest = zeros(imageSize);
nrOfQueries = 0;
pixelMeanMatrix = zeros(nrOfSubmasks);
for row = 1:nrOfSubmasks
  for col = 1:nrOfSubmasks

    [Aest, mask, pixelMean] = updatepixels(Aest, mask, imageSize, submask, ...
      submasksize, submasksum, row, col, one_by_integer_lookup_table);

    nrOfQueries = nrOfQueries + 1;

    % Remember pixel mean
    pixelMeanMatrix(row, col) = pixelMean;
  end
end

idx0 = 1:nrOfSubmasks;
idx1 = 1:(nrOfSubmasks-1);
idx2 = 2:nrOfSubmasks;
pixelMeanDiff1 = abs(pixelMeanMatrix(idx1, idx0) - pixelMeanMatrix(idx2, idx0));
pixelMeanDiff2 = abs(pixelMeanMatrix(idx0, idx1) - pixelMeanMatrix(idx0, idx2));
[ignore, sortIndex1] = sort(pixelMeanDiff1(:), 'descend');
[ignore, sortIndex2] = sort(pixelMeanDiff2(:), 'descend');

submasksize2  = submasksize  / 2;
nrOfSubmasks2 = ceil(imageSize / submasksize2);
[submask2, submasksum2] = getsubmask(submasksize2);

refined = zeros(nrOfSubmasks);
i1 = 1; i2 = 1;
iMax = nrOfSubmasks*(nrOfSubmasks-1);
while 1
  if (i1 <= iMax) && ( (i2 > iMax) || ...
      (pixelMeanDiff1(sortIndex1(i1)) > pixelMeanDiff2(sortIndex2(i2))))
    [rows(1), cols(1)] = ind2sub(nrOfSubmasks-1, sortIndex1(i1));
    rows(2) = rows(1) + 1; cols(2) = cols(1);
    i1 = i1 + 1;
  elseif (i2 <= iMax)
    [rows(1), cols(1)] = ind2sub(nrOfSubmasks, sortIndex2(i2));
    rows(2) = rows(1); cols(2) = cols(1) + 1;
    i2 = i2 + 1;
  else
    return
  end

  for k = 1:2
    if refined(rows(k), cols(k))
      continue;
    end
    for j = [0 0 1 1; 0 1 0 1]
      subrow = 2*rows(k)-1 + j(1);
      subcol = 2*cols(k)-1 + j(2);

      if (subrow <= nrOfSubmasks2) && (subcol <= nrOfSubmasks2)
        [Aest, mask] = updatepixels(Aest, mask, imageSize, submask2, ...
          submasksize2, submasksum2, subrow, subcol, one_by_integer_lookup_table);

        nrOfQueries = nrOfQueries + 1;
        if nrOfQueries >= queryLimit
          %q = (i1+i2) / (2*iMax)
          return
        end
      end
    end
    refined(rows(k), cols(k)) = 1;
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Aest, mask, pixelMean] = updatepixels(Aest, mask, imageSize, submask, ...
  submasksize, submasksum, row, col, one_by_integer_lookup_table)

offset1 = 1+(row-1)*submasksize;
offset2 = 1+(col-1)*submasksize;
[index1, index2, n1, n2, border] = getsubmaskindex(imageSize, submasksize, offset1, offset2);

% Get pixel sum
[mask, masksum] = setmask(mask, index1, index2, n1, n2, border, submasksum, submask);
pixelSum = queryImage(logical(mask));
mask = resetmask(mask, index1, index2);

% Set pixels
pixelMean = pixelSum * one_by_integer_lookup_table(masksum);
Aest = setpixels(Aest, index1, index2, pixelMean);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [index1, index2, n1, n2, border] = getsubmaskindex(imageSize, submasksize, offset1, offset2)

border = 0;

idx1 = (offset1+submasksize-1);
if idx1 > imageSize
  idx1 = imageSize;
  border = 1;
  n1 = idx1 - offset1 + 1;
  if n1 < 0
    n1 = 0;
  end
else
  n1 = submasksize;
end
index1 = offset1:idx1;

idx2 = (offset2+submasksize-1);
if idx2 > imageSize
  idx2 = imageSize;
  border = 1;
  n2 = idx2 - offset2 + 1;
  if n2 < 0
    n2 = 0;
  end
else
  n2 = submasksize;
end
index2 = offset2:idx2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [mask, masksum] = setmask(mask, index1, index2, n1, n2, border, submasksum, submask)

mask(index1, index2) = submask(1:n1, 1:n2);
if border
  masksum = sum(sum(mask));
  if masksum == 0
    mask(index1(1), index2(1)) = 1;
    masksum = 1;
  end
else
  masksum = submasksum;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function mask = resetmask(mask, index1, index2)
mask(index1, index2) = 0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Aest = setpixels(Aest, index1, index2, pixelValue)
Aest(index1, index2) = pixelValue;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [submask, submasksum] = getsubmask(submasksize)
switch submasksize
  case {1, 2, 3, 4, 5}
    submask    = ones(submasksize);
    submasksum = submasksize^2;

  case {6, 7, 8, 10}
    submask = reshape(1:(submasksize*(submasksize+1)), submasksize+1, submasksize);
    submask = mod(submask(1:submasksize, :), 2) == 1;
    submasksum = sum(submask(:));
    
  case 12
    submask = reshape(1:(submasksize*(submasksize+1)), submasksize+1, submasksize);
    submask = mod(submask(1:submasksize, :), 3) == 1;
    submasksum = sum(submask(:));
    
  case {14, 16}
    submask = reshape(1:(submasksize*(submasksize)), submasksize, submasksize);
    submask = mod(submask(1:submasksize, :), 3) == 1;
    submasksum = sum(submask(:));
end