function retinex( filename, outname )
% 说明：该算法仅对RGB色度空间有效，即输入图像波段数应该为3
% filename - 输入图像文件路径，格式为ENVI标准格式
% outname - 输出图像文件路径， 格式为ENVI标准格式

clc
tic

origImg = readENVIImg( filename );
figure; imshow( origImg );
title('Original image');

%% 参数设置
k1=225;
k2=3000; 
k3=80000;
k4=10000;
r=81;
c = 0.01;
nn=floor((r+1)/2);   % 计算中心

Img1 = double( origImg(:,:,1) );
Img2 = double( origImg(:,:,2) );
Img3 = double( origImg(:,:,3) );
[m, n]=size( Img1 );
b = zeros(r,r);
bb = zeros(r,r);
bbb = zeros(r,r);
%% 高斯函数 共有3个
for i=1:r
    for j=1:r
        b(i,j) =exp(-((i-nn)^2+(j-nn)^2)/k1)/k4; 
   end
end
tiaojie1 = 1/sum(sum(b));
b = tiaojie1.*b;

for i=1:r
    for j=1:r
        bb(i,j) =exp(-((i-nn)^2+(j-nn)^2)/k2)/k4; 
   end
end
tiaojie2 = 1/sum(sum(bb));
bb = tiaojie2.*bb;
for i=1:r
    for j=1:r
        bbb(i,j) =exp(-((i-nn)^2+(j-nn)^2)/k3)/k4; 
   end
end
tiaojie3 = 1/sum(sum(bbb));
bbb = tiaojie3.*bbb;

%% 对R分量的处理
K=imfilter(Img1,b);
KK=imfilter(Img1,bb);
KKK=imfilter(Img1,bbb);

G = zeros(m,n);
for i=1:m
    for j=1:n      
       G(i,j)=1/3*(log(Img1(i,j)+c)-log(K(i,j)+c));
        G(i,j)=1/3*(log(Img1(i,j)+c)-log(KK(i,j)+c))+G(i,j);
         G(i,j)=1/3*(log(Img1(i,j)+c)-log(KKK(i,j)+c))+G(i,j);
    end
end
G = 3*G+50;
mi=min(min(G));
ma=max(max(G));
L=(G-mi)*255/(ma-mi);

%% 对G分量的处理
K=imfilter(Img2,b);
KK=imfilter(Img2,bb);
KKK=imfilter(Img2,bbb);

for i=1:m
    for j=1:n      
       G(i,j)=1/3*(log(Img2(i,j)+c)-log(K(i,j)+c));
        G(i,j)=1/3*(log(Img2(i,j)+c)-log(KK(i,j)+c))+G(i,j);
         G(i,j)=1/3*(log(Img2(i,j)+c)-log(KKK(i,j)+c))+G(i,j);
    end
end
G = 3*G+50;
mi=min(min(G));
ma=max(max(G));
LL=(G-mi)*255/(ma-mi);
%% 对B分量的处理
K=imfilter(Img3,b);
KK=imfilter(Img3,bb);
KKK=imfilter(Img3,bbb);

for i=1:m
    for j=1:n      
       G(i,j)=1/3*(log(Img3(i,j)+c)-log(K(i,j)+c));
        G(i,j)=1/3*(log(Img3(i,j)+c)-log(KK(i,j)+c))+G(i,j);
         G(i,j)=1/3*(log(Img3(i,j)+c)-log(KKK(i,j)+c))+G(i,j);
    end
end
G = 3*G+50;
mi=min(min(G));
ma=max(max(G));
LLL=(G-mi)*255/(ma-mi);
%% 对彩色图像的综合处理 
msr = cat(3, L, LL, LLL);

writeENVIImg(msr, outname, 'float');
toc
end

% ---------- sub functions ----------
% Author : Jacory Gao
function writeENVIImg(filename, img, format)
% write envi standard format images, only support bsq format and can't
% specify header offset.
% pamameters:
%   filename  :  output file path.
%   img       :  image variable.
%   format    :  format string

expectedFormats = { 'uint8', 'int16','int32', 'float32'};
parser = inputParser;
addRequired(parser, 'filename');
addRequired(parser, 'img');
addRequired(parser, 'format', @(x) any(validatestring(x, expectedFormats)));
parse( parser, filename, img, format);

basename = strsplit(filename, '.');
if size( basename ) == 1
    hdrpath = strjoin( [basename, '.hdr'], '');
else
    hdrpath = strjoin( [basename(1:end-1), '.hdr'], '');
end
hdrfile = fopen(hdrpath, 'w');


[lines, samples, bandCount] = size(img);

switch format
    case 'uint8'
        dataType = 1;
    case 'int16'
        dataType = 2;
    case 'int32'
        dataType = 3;
    case 'float32'
        dataType = 4;
    otherwise
        dataType = 4;
end

fprintf(hdrfile, '%s\n', 'ENVI description = { matlab writeENVIImg function }');
fprintf(hdrfile, '%s\n', ['samples = ', int2str(samples)]);
fprintf(hdrfile, '%s\n', ['lines = ', int2str(lines)]);
fprintf(hdrfile, '%s\n', ['bands = ', int2str(bandCount)]);
fprintf(hdrfile, '%s\n', 'header offset = 0');
fprintf(hdrfile, '%s\n', 'file type = ENVI Standard');
fprintf(hdrfile, '%s\n', ['data type = ', int2str(dataType)]);
fprintf(hdrfile, '%s\n', 'interleave = bsq');

fclose(hdrfile);

outfile = fopen(filename, 'w');
img = permute(img, [2,1,3]);
fwrite(outfile, img, format);
fclose(outfile);

end

% Author : Jacory Gao
function [Image_Cube, lines, samples, bandCount, datatype] = readENVIImg( filepath );
    % read ENVI standard image file into matlab 
    % filepath: image file path
    % return：Image_Cube - 3 dimensional array of image cube
	%				 lines - image height
	%				 samples - image width
	%				 bandCount - band count of image file
	%				 datatype - data type of image file
    
    %% read header file
    basename = strsplit(filepath, '.');
    if size( basename ) == 1
        hdrpath = strjoin( [basename, '.hdr'], '');
    else
        hdrpath = strjoin( [basename(1:end-1), '.hdr'], '');
    end
    fid = fopen( hdrpath, 'r' );
    
    linestr = '';
    
    samples = 0;
    lines = 0;
    bandCount = 0;
    datatype = 0;
    
    while ~feof(fid)
        linestr = fgetl(fid);
        linearr = strsplit(linestr, '=');
        if  strcmp(strtrim(linearr(1)), 'samples')
            samples = str2double( cell2mat( strtrim( linearr(2) ) ) );
        elseif strcmp(strtrim(linearr(1)), 'lines')
            lines = str2double( cell2mat( strtrim( linearr(2) ) ) );
        elseif strcmp(strtrim(linearr(1)), 'bands')
            bandCount = str2double( cell2mat( strtrim( linearr(2) ) ) );
        elseif strcmp(strtrim(linearr(1)), 'data type')
            datatype = str2double( cell2mat( strtrim( linearr(2) ) ) );
        end
    end
    fclose(fid);
    
    if samples == 0 || lines == 0 || bandCount == 0 || datatype == 0
        error('read hdr file error.');
    end
    
    %% read image file
    fid = fopen(filepath,'r');
    switch datatype
        case 1                       % 8 bit
            PixData = fread(fid, 'uint8');
            len = length( PixData );
            if len ~= samples * lines * bandCount
                error('image input size does not match the metadata description!');
            end
            Image_Cube = reshape(PixData,[samples,lines,bandCount]);
        case 2                       % 16 bit
            PixData = fread(fid,'int16'); 
            len = length( PixData );
            if len ~= samples * lines * bandCount
                error('image input size does not match the metadata description!');
            end
            Image_Cube = reshape(PixData,[samples,lines,bandCount]);
        case 3                       % 32 bit
            PixData = fread(fid,'int32');
            len = length( PixData );
            if len ~= samples * lines * bandCount
                error('image input size does not match the metadata description!');
            end
            Image_Cube = reshape(PixData,[samples,lines,bandCount]);
        case 4                       % 32 bit floating point
            PixData = fread(fid,'float32');
            len = length( PixData );
            if len ~= samples * lines * bandCount
                error('image input size does not match the metadata description!');
            end
            Image_Cube = reshape(PixData,[samples,lines,bandCount]);
        otherwise
            error('data type is wrong.');
    end
    Image_Cube = permute(Image_Cube,[2 1 3]);
end