close all; clc
%% =====================================================================================================
% Parameters to be modified

% read an image with saturated regions
I_1 = imread('picture 1.jpg');

% read an image without saturation
I_2 = imread('picture 2.jpg');

% The pixels of the image with saturated regions whose values are between 255 and 255-level
% are replaced by the corresponding pixels from the non-saturated image.
% Try adjusting the value of "level" to see the effect on the HDR image.
level = 5;




% XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
%
%                             DO NOT MODIFY THE SCRIPT BEYOND THIS POINT
%
% XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX



%% =====================================================================================================
% Process data

% Consdier the Green channel
I1 = I_1(:,:,1);
I2 = I_2(:,:,1);

R = I1;

% find the saturation level of R - adjustment by 5 counts, (max(R) - 5)
m = max(max(R))-level;

% find the indices of the saturated pixels;
i = find(R > m);

% fusion of images by replacing the saturated pixel i in I1 by that from image I2
R(i) = I2(i);

%% =====================================================================================================
% Plot data

figure('color','w','position',[250 250 1200 400])

% show image 1
subplot(1,3,1)
imagesc(I1,[0 255])
colormap(gray)
axis square
title('Image with saturation')
xlabel('Position x (pixel)')
ylabel('Position y (pixel)')

% show image 2
subplot(1,3,2)
imagesc(I2,[0 255])
colormap(gray)
axis square
title('Image without saturation')
xlabel('Position x (pixel)')
ylabel('Position y (pixel)')

% show the combined image
subplot(1,3,3)
imagesc(R,[0 255])
colormap(gray)
axis square
title('HDR image')
xlabel('Position x (pixel)')
ylabel('Position y (pixel)')