function jointDistribution = jointDistFromImages(image1, image2, powerOfTwo)

  bins = 2 ^ powerOfTwo;
  valuesPerBin = 256 / bins;

  rows = size(image1, 1);
  cols = size(image1, 2);

  jointDistribution = zeros(bins, bins);

  for row = 1 : rows
    for col = 1 : cols
      binCol = ceil(double(image1(row, col) + 1) / double(valuesPerBin));
      binRow = ceil(double(image2(row, col) + 1) / double(valuesPerBin));
      jointDistribution(binRow, binCol) = jointDistribution(binRow, binCol) + 1;
    end
  end

  jointDistribution = jointDistribution * (1.0 / double(rows * cols));

end
