% ImOriginal:原始图像

% type:'add'表示分量是加性的，如雾天图像；'mult'表示分量是乘性的，如对照度的估计

[m,n,z] = size(ImOriginal);

ImOut = uint8(zeros(m,n,z));

for i = 1:z

    if strcmp(type,'add')

        ImChannel = double(ImOriginal(:,:,i))+eps;

    elseif strcmp(type,'mult')

        ImChannel = log(double(ImOriginal(:,:,i))+eps);

    else

        error('type must be ''add'' or ''mult''');

    end

    ImOut(:,:,i) = EnhanceOneChannel(ImChannel);

end

ImOut = max(min(ImOut,255), 0);

end

function ImOut = EnhanceOneChannel(ImChannel)

% 计算计算单个通道的反射分量

% 1.对全图进行照射分量估计

% 2.减去照射分量

% 3.灰度拉伸

ImChannel = ImChannel./max(ImChannel(:));

ImRetinex = round(exp(ImChannel.*5.54));

ImOut = uint8(ImRetinex);

end

function Retinex = retinex_frankle_mccann(L, nIterations)

% RETINEX_FRANKLE_McCANN:

%         Computes the raw Retinex output from an intensity image, based on the

%         original model described in:

%         Frankle, J. and McCann, J., "Method and Apparatus for Lightness Imaging"

%         US Patent #4,384,336, May 17, 1983

%

% INPUT:  L           - logarithmic single-channel intensity image to be processed

%         nIterations - number of Retinex iterations

%

% OUTPUT: Retinex     - raw Retinex output

%

% NOTES:  - The input image is assumed to be logarithmic and in the range [0..1]

%         - To obtain the retinex "sensation" prediction, a look-up-table needs to

%         be applied to the raw retinex output

%         - For colour images, apply the algorithm individually for each channel

%

% AUTHORS: Florian Ciurea, Brian Funt and John McCann.

%          Code developed at Simon Fraser University.

%

% For information about the code see: Brian Funt, Florian Ciurea, and John McCann

% "Retinex in Matlab," by Proceedings of the IS&T/SID Eighth Color Imaging

% Conference: Color Science, Systems and Applications, 2000, pp 112-121.

%

% paper available online at http://www.cs.sfu.ca/~colour/publications/IST-2000/

%

% Copyright 2000. Permission granted to use and copy the code for research and

% educational purposes only.  Sale of the code is not permitted. The code may be

% redistributed so long as the original source and authors are cited.

global RR IP OP NP Maximum

RR = L;

Maximum = max(L(:));                                 % maximum color value in the image

[nrows, ncols] = size(L);

shift = 2^(fix(log2(min(nrows, ncols)))-1);          % initial shift

OP = Maximum*ones(nrows, ncols);                     % initialize Old Product

while (abs(shift) >= 1)

   for i = 1:nIterations

      CompareWith(0, shift);                         % horizontal step

      CompareWith(shift, 0);                         % vertical step

   end

   shift = -shift/2;                                 % update the shift

end

Retinex = NP;

function CompareWith(s_row, s_col)

global RR IP OP NP Maximum

IP = OP;

if (s_row + s_col > 0)

   IP((s_row+1):end, (s_col+1):end) = OP(1:(end-s_row), 1:(end-s_col)) + ...

   RR((s_row+1):end, (s_col+1):end) - RR(1:(end-s_row), 1:(end-s_col));

else

   IP(1:(end+s_row), 1:(end+s_col)) = OP((1-s_row):end, (1-s_col):end) + ...

   RR(1:(end+s_row),1:(end+s_col)) - RR((1-s_row):end, (1-s_col):end);

end

IP(IP > Maximum) = Maximum;                          % The Reset operation

NP = (IP + OP)/2;                                    % average with the previous Old Product

OP = NP;                                             % get ready for the next comparison

function Retinex = retinex_mccann99(L, nIterations)

global OPE RRE Maximum
[nrows ncols] = size(L);                             % get size of the input image
nLayers = ComputeLayers(nrows, ncols);               % compute the number of pyramid layers
nrows = nrows/(2^nLayers);                           % size of image to process for layer 0
ncols = ncols/(2^nLayers);
if (nrows*ncols > 25)                                % not processing images of area > 25
  error('invalid image size.')                       % at first layer
end
Maximum = max(L(:));                                 % maximum color value in the image
OP = Maximum*ones([nrows ncols]);                    % initialize Old Product
for layer = 0:nLayers
   RR = ImageDownResolution(L, 2^(nLayers-layer));   % reduce input to required layer size
  
   OPE = [zeros(nrows,1) OP zeros(nrows,1)];         % pad OP with additional columns
   OPE = [zeros(1,ncols+2); OPE; zeros(1,ncols+2)];  % and rows
   RRE = [RR(:,1) RR RR(:,end)];                     % pad RR with additional columns
   RRE = [RRE(1,:); RRE; RRE(end,:)];                % and rows
  
   for iter = 1:nIterations
     CompareWithNeighbor(-1, 0);                     % North
     CompareWithNeighbor(-1, 1);                     % North-East
     CompareWithNeighbor(0, 1);                      % East
     CompareWithNeighbor(1, 1);                      % South-East
     CompareWithNeighbor(1, 0);                      % South
     CompareWithNeighbor(1, -1);                     % South-West
     CompareWithNeighbor(0, -1);                     % West
     CompareWithNeighbor(-1, -1);                    % North-West
   end
   NP = OPE(2:(end-1), 2:(end-1));
   OP = NP(:, [fix(1:0.5:ncols) ncols]);             %%% these two lines are equivalent with
   OP = OP([fix(1:0.5:nrows) nrows], :);             %%% OP = imresize(NP, 2) if using Image
   nrows = 2*nrows; ncols = 2*ncols;                 % Processing Toolbox in MATLAB
end
Retinex = NP;

function CompareWithNeighbor(dif_row, dif_col)
global OPE RRE Maximum

% Ratio-Product operation
IP = OPE(2+dif_row:(end-1+dif_row), 2+dif_col:(end-1+dif_col)) + ...
     RRE(2:(end-1),2:(end-1)) - RRE(2+dif_row:(end-1+dif_row), 2+dif_col:(end-1+dif_col));
    
IP(IP > Maximum) = Maximum;                          % The Reset step

% ignore the results obtained in the rows or columns for which the neighbors are undefined
if (dif_col == -1) IP(:,1) = OPE(2:(end-1),2); end
if (dif_col == +1) IP(:,end) = OPE(2:(end-1),end-1); end
if (dif_row == -1) IP(1,:) = OPE(2, 2:(end-1)); end
if (dif_row == +1) IP(end,:) = OPE(end-1, 2:(end-1)); end
NP = (OPE(2:(end-1),2:(end-1)) + IP)/2;              % The Averaging operation
OPE(2:(end-1), 2:(end-1)) = NP;

function Layers = ComputeLayers(nrows, ncols)
power = 2^fix(log2(gcd(nrows, ncols)));              % start from the Greatest Common Divisor
while(power > 1 & ((rem(nrows, power) ~= 0) | (rem(ncols, power) ~= 0)))
   power = power/2;                                  % and find the greatest common divisor
end                                                  % that is a power of 2
Layers = log2(power);

function Result = ImageDownResolution(A, blocksize)
[rows, cols] = size(A);                              % the input matrix A is viewed as
result_rows = rows/blocksize;                        % a series of square blocks
result_cols = cols/blocksize;                        % of size = blocksize
Result = zeros([result_rows result_cols]);
for crt_row = 1:result_rows                          % then each pixel is computed as
   for crt_col = 1:result_cols                       % the average of each such block
      Result(crt_row, crt_col) = mean2(A(1+(crt_row-1)*blocksize:crt_row*blocksize, ...
                                       1+(crt_col-1)*blocksize:crt_col*blocksize));
   end
en