function tonemap()
  compute('atriumnight');
  compute('church');
  compute('doll_doll');
  compute('lips');
  compute('rend07_o85C');
  compute('rosette');
end

function compute(imageFile)
  RawImage = getpfmraw(strcat(imageFile, '.pfm'));
  [HOriginal, SOriginal, VOriginal] = rgb2hsv(RawImage);
  
  
  VLinear = linearmap(VOriginal);
  LinearImage = hsv2rgb(HOriginal, SOriginal, VLinear);
  imwrite(LinearImage, strcat(imageFile, '_linear.bmp'));
  
  VSimple = simplemap(VOriginal, 0.18);
  SimpleImage = hsv2rgb(HOriginal, SOriginal, VSimple);
  imwrite(SimpleImage, strcat(imageFile, '_key_0.18.bmp'));
   
  
  %VSimple = simplemap(VOriginal, 0.18);
  
  VGlobal = globalmap(VSimple);
  GlobalImage = hsv2rgb(HOriginal, SOriginal, VGlobal);  
  GlobalImage = imadjust(GlobalImage, [], [], 1.0 / 2.2);
  imwrite(GlobalImage, strcat(imageFile, '_global.bmp'));  
  
  VDodgeBurn = dodgeburn(VSimple, 0.18);
  DodgeBurnImage = hsv2rgb(HOriginal, SOriginal, VDodgeBurn);
  DodgeBurnImage = imadjust(DodgeBurnImage, [], [], 1.0 / 2.2);
  imwrite(DodgeBurnImage, strcat(imageFile, '_local.bmp'));
end

function VGlobal = globalmap(VOriginal)
  VGlobal = VOriginal;
  for i = 1 : size(VOriginal, 1) * size(VOriginal, 2)
    VGlobal(i) = VOriginal(i) / (1.0 + VOriginal(i));
  end
end

function VLinear = linearmap(VOriginal)
  maxV = max(max(VOriginal));  
  VLinear = VOriginal / maxV * 255.0;   
end

function VSimple = simplemap(VOriginal, a)
  epsilon = 0.00001;
  lw = exp(mean2(log(epsilon + VOriginal)));
  VSimple = VOriginal * (a / lw);
end

function VDodgeBurn = dodgeburn(VOriginal, a)
  numLevel = 12;
  low = 1;  
  threshold = 0.05;
  phi = 8.0;
  alpha = [1.0, 1.2];
  r = ones(1, numLevel + 1);  
  for i = 1 : numLevel + 1
    r(1, i) = floor(low * alpha(2) ^ i);
  end   
  V1 = zeros(size(VOriginal, 1), size(VOriginal, 2), numLevel + 1);
  V2 = zeros(size(VOriginal, 1), size(VOriginal, 2), numLevel + 1);
  for i = 1 : numLevel
    R1 = fspecial('gaussian', r(1, i) * 2.0 + 1.0, r(1, i) / sqrt(2.0));    
    R2 = fspecial('gaussian', r(1, i + 1) * 2.0 + 1.0, r(1, i + 1) / sqrt(2.0));
    V1(:, :, i) = imfilter(VOriginal, R1);
    V2(:, :, i) = imfilter(VOriginal, R2);
  end
  VDodgeBurn = VOriginal;
  for i = 1 : size(VOriginal, 1)
    for j = 1 : size(VOriginal, 2)
      v1 = VOriginal(i, j);
      for k = numLevel : -1 : 1
        v = abs(V1(i, j, k) - V2(i, j, k)) / (V1(i, j, k) + 2.0 ^ phi * a / r(1, k) ^ 2.0);   
        if v < threshold && abs(V1(i, j, k) - VOriginal(i, j)) / (min(V1(i, j, k), VOriginal(i, j)) + 0.0001) < 0.8
          v1 = V1(i, j, k);          
          break;
        end
      end
      VDodgeBurn(i, j) = VOriginal(i, j) / (1.0 + v1);
    end
  end  
end