// A Practical Introduction to Metropolis Light Transport
// 2.2 Color Images
void makeHistogram(int w, int h, PVector[] F, PVector[] histogram, int mutations) {
  // Create an initial sample point
  DomainLocation X = new DomainLocation();
  X.xloc = randomInteger(0, w - 1);
  X.yloc = randomInteger(0, h - 1);
  PVector colorX = new PVector(F[X.yloc * w + X.xloc].x, F[X.yloc * w + X.xloc].y, F[X.yloc * w + X.xloc].z);
  float Fx = luminance(colorX);
  if (Fx > 0)
    colorX.mult(1.0 / Fx);

  DomainLocation Y = new DomainLocation();

  // Create a histogram of values using Metropolis sampling.
  for (int i = 0; i < mutations; ++i) {
    // choose a tentative next sample according to T.
    mutateAccordingToT(w, h, X, Y);
    float Tyx = T(w, h, Y, X);
    float Txy = T(w, h, X, Y);

    PVector colorY = F[Y.yloc * w + Y.xloc];
    float Fy = luminance(colorY);
    float Axy = min(1.0, (Fy * Txy) / (Fx * Tyx));  // equation 2.
    if (randomReal(0.0, 1.0) < Axy) {
      X.set(Y);
      Fx = Fy;
      colorX.set(colorY.x / Fy, colorY.y / Fy, colorY.z / Fy);
    }
    histogram[X.yloc * w + X.xloc].add(colorX);
  }
}

void mutateAccordingToT(int w, int h, DomainLocation X, DomainLocation Y) {
  Y.xloc = randomInteger(0, w - 1);
  Y.yloc = randomInteger(0, h - 1);
}

float T(int w, int h, DomainLocation X, DomainLocation Y) {
  return 1.0 / (w * h);
}

float luminance(PVector v) {
  return 0.299 * v.x + 0.587 * v.y + 0.114 * v.z;
}