function BasicHE(path)

input = imread(path);
[y,x,z] = size(input);

figure;
subplot(2,3,1); imshow(input); title('Before');
subplot(2,3,2); imhist(input(:)); title('Histogram before');

lambda = 1;
gamma = 1024;

if z==1
    output = HEmodule(input,x,y,lambda,gamma);
else
    % Y=T(Y)   Cb=Cb   Cr=Cr
    output = YUVconvert(input,x,y);
    output(:,:,1) = HEmodule(uint8(output(:,:,1)),x,y,lambda,gamma);
    output = RGBconvert(output,x,y);
end

subplot(2,3,4); imshow(output); title('After');
subplot(2,3,5); imhist(output(:)); title('Histogram after');

end

function output = YUVconvert(input,x,y)

output = zeros(y,x,3);

for i=1:x
    for j=1:y
        output(j,i,1) = 0.299*input(j,i,1) + 0.587*input(j,i,2) + 0.114*input(j,i,3);
        output(j,i,2) = 128 - 0.168736*input(j,i,1) - 0.331264*input(j,i,2) + 0.5*input(j,i,3);
        output(j,i,3) = 128 + 0.5*input(j,i,1) - 0.418688*input(j,i,2) - 0.081312*input(j,i,3);
    end
end

end

function output = RGBconvert(input,x,y)

output = zeros(y,x,3,'uint8');

for i=1:x
    for j=1:y
        output(j,i,1) = input(j,i,1) + 1.402*(input(j,i,3)-128);
        output(j,i,2) = input(j,i,1) - 0.344136*(input(j,i,2)-128) - 0.714136*(input(j,i,3)-128);
        output(j,i,3) = input(j,i,1) + 1.772*(input(j,i,2)-128);
    end
end

end

function output = HEmodule(input,x,y,lambda,gamma)

H = imhist(input);

D = zeros(255,256);
for i=1:256
    for j=1:255
        if i==j
            D(j,i) = -1;
            D(j,i+1) = 1;
        end
    end
end
I = eye(256);
S = (1+lambda)*I + gamma*(D'*D);
S = inv(S);

%HWS_print(S);

Hi = zeros(1,256,'uint8');
for i=1:256
    Hi(i) = sum(H(1:i))*255/(x*y);
end
subplot(2,3,3); plot(Hi); title('Standard transformation function');

u = 0:255;
%C = (double(Hi) + lambda*u)/(1+lambda);
C = S*(double(Hi) + lambda*u)';
C = uint8(C);
subplot(2,3,6); plot(C); title('Smooth transformation function');

output = zeros(y,x,'uint8');
for i=1:x
    for j=1:y
        f = double(input(j,i))+1;
        output(j,i) = C(f);
    end
end

end

function HWS_print(S)

HWS = S*1024;
HWS = uint8(HWS);

fname = 'HWS.txt';
fid = fopen(fname,'w');
if fid ~= -1
    for i=1:256
        for j=1:256
            fprintf(fid,'\t%d\t:\t',(i-1)*256+(j-1));
            temp = de2bi(HWS(i,j),5);
            for k=5:-1:1
                fprintf(fid,'%d',temp(k));
            end
            fprintf(fid,';\n');
        end
    end
    fclose(fid);
end

end
