function imOut = reduceDR(imIn,alpha,beta,s,filename);

colorlayers = size(imIn,3);

if colorlayers == 3
    Lin = 0.30*imIn(:,:,1) + 0.59*imIn(:,:,2) + 0.11*imIn(:,:,3);
else
    Lin = imIn;
end
    
Lin(Lin == 0) = inf;
Lin(isinf(Lin)) = min2(Lin);
H = log(Lin);

if exist('filename')
    Hmod = H-min2(H);
    Hmod = Hmod/max2(Hmod);
    disp('Writing logarithmic input luminance images to JPEG.');
    imwrite(Hmod,                [filename '_H_l.jpg'],'JPEG','quality',100);
    imwrite(reduce(Hmod),        [filename '_H_m.jpg'],'JPEG','quality',100);
    imwrite(reduce(reduce(Hmod)),[filename '_H_s.jpg'],'JPEG','quality',100);
    figure(1); imshow(Hmod);
end    

disp(' ');
disp('Generating gradient attenuation map...');
phi = gradientAttenuator(H,alpha,beta); 
disp('Gradient attenuation map complete.');
disp(' ');

if exist('filename')
    phimod = phi-min2(phi);
    phimod = phimod*0.5/mean2(phimod);
    phimod(phimod > 1) = 1;
    disp('Writing gradient attenuation functions to JPEG.');
    imwrite(phimod,                [filename '_phi_l.jpg'],'JPEG','quality',100);
    imwrite(reduce(phimod),        [filename '_phi_m.jpg'],'JPEG','quality',100);
    imwrite(reduce(reduce(phimod)),[filename '_phi_s.jpg'],'JPEG','quality',100);
    figure(2); imshow(phimod);
end    

delHx = padarray([H(:,2:end)-H(:,1:(end-1))],[0 1],0,'post');
delHy = padarray([H(2:end,:)-H(1:(end-1),:)],[1 0],0,'post');

Gx = delHx.*phi;
Gy = delHy.*phi;

divGx = padarray([Gx(:,2:end)-Gx(:,1:(end-1))],[0 1],0,'pre');
divGy = padarray([Gy(2:end,:)-Gy(1:(end-1),:)],[1 0],0,'pre');
divG  = divGx + divGy;

if exist('filename')
    divGmod = divG-min2(divG);
    divGmod = divGmod/max2(divGmod);
    disp('Writing logarithmic output luminance images to JPEG.');
    imwrite(divGmod,                [filename '_I_l.jpg'],'JPEG','quality',100);
    imwrite(reduce(divGmod),        [filename '_I_m.jpg'],'JPEG','quality',100);
    imwrite(reduce(reduce(divGmod)),[filename '_I_s.jpg'],'JPEG','quality',100);
    figure(3); imshow(divGmod);
end    

% divdelHx = padarray([delHx(:,2:end)-delHx(:,1:(end-1))],[0 1],0,'pre');
% divdelHy = padarray([delHy(2:end,:)-delHy(1:(end-1),:)],[1 0],0,'pre');
% divdelH  = divdelHx + divdelHy;

disp(' ');
disp('Solving Poisson''s Equation...');
I = poissonSolver2D(divG);
disp('Optimal luminance image found.');
disp(' ');

if exist('filename')
    Imod = I-min2(I);
    Imod = Imod/max2(Imod);
    disp('Writing logarithmic output luminance images to JPEG.');
    imwrite(Imod,                [filename '_I_l.jpg'],'JPEG','quality',100);
    imwrite(reduce(Imod),        [filename '_I_m.jpg'],'JPEG','quality',100);
    imwrite(reduce(reduce(Imod)),[filename '_I_s.jpg'],'JPEG','quality',100);
    figure(4); imshow(Imod);
end    

Lout = exp(I-max2(I));

if exist('filename')
    Loutmod = Lout-min2(Lout);
    Loutmod = Loutmod/max2(Loutmod);
    disp('Writing output luminance images to JPEG.');
    imwrite(Loutmod,                [filename '_Lout_l.jpg'],'JPEG','quality',100);
    imwrite(reduce(Loutmod),        [filename '_Lout_m.jpg'],'JPEG','quality',100);
    imwrite(reduce(reduce(Loutmod)),[filename '_Lout_s.jpg'],'JPEG','quality',100);
    figure(5); imshow(Lout);
end    

if colorlayers > 1
    disp('Generating colored image...');
    imOut = repmat(Lout,[1 1 3]).*(imIn./repmat(Lin,[1 1 3])).^s;
    imOut = imOut/max(max2(imOut));

    if exist('filename')
        imOutmod = imOut-min(min2(imOut));
        imOutmod = imOutmod/max(max2(imOutmod));
        disp('Writing output color images to JPEG.');
        imwrite(imOutmod,                          [filename '_ldr_l.jpg'],'JPEG','quality',100);
        imwrite(multireduce(imOutmod),             [filename '_ldr_m.jpg'],'JPEG','quality',100);
        imwrite(multireduce(multireduce(imOutmod)),[filename '_ldr_s.jpg'],'JPEG','quality',100);
        figure(6); imshow(imOutmod);
    end    
    
else
    imOut = Lout;
end

disp(['The initial image had a dynamic range of ' ...
            num2str((max2(H)-min2(H))/log(10)) ' log units.']);
disp(['The final image has a dynamic range of '   ...
            num2str((max2(I)-min2(I))/log(10)) ' log units.']);
