clear
clc

%%Parameter n of requires to images enhancement is chosen as 20%
pp = 20.0;
%%Percentage of gray level
pplevel = 98.0;
%%Gray levels
level = 4095.0;
%%Output value of enhancement
out = 1700.0;

%%INPUTS:
%%Input directory
dirInput  = '...../MHD/'; 
%%Output directory
dirOutput = '...../Enhancement/';
%%Input volume name
ImageFile = '....';

mhdFile = [dirInput ImageFile '.mhd'];
rawFile = [dirInput ImageFile '.raw'];
mhdText = fileread(mhdFile);
b1 = strfind(mhdText, 'DimSize = ');
b2 = strfind(mhdText, 'AnatomicalOrientation = ');
imageDims = mhdText(b1+9:b2-1);
dims = regexp(imageDims,'\d*','Match');
ROW = str2double(cell2mat(dims(1)));
COL = str2double(cell2mat(dims(2)));
SLI = str2double(cell2mat(dims(3)));

fp  = fopen(rawFile,'rb');
ima = fread(fp,ROW*COL*SLI,'int16');
fclose(fp);
ima = reshape(ima, [ROW COL SLI]);
imax = max(max(max(ima)));
imin = min(min(min(ima)));
ima = level*(ima - imin)/(imax - imin);

sli2 = int16(SLI/2);
n = int16(SLI*pp/100);
k0 = sli2 - int16(n/2);
k1 = sli2 + int16(n/2);
h = zeros(n+1,level);

ii = 1;
for k=k0:k1,
    img = ima(:,:,k);
    for i=1:ROW,
        for j=1:COL,
            kk = int16(img(i,j));
            h(ii,kk+1) = h(ii,kk+1) + 1;
        end
    end
    ii = ii + 1;
end
H = sum(h)/3;

level0 = 1;
level1 = 1300;
level2 = 2100;
level3 = 2150;
level4 = double(int16(level*pplevel/100));

f = [H(level0), H(level0)+H(level1), H(level0)+H(level1)+H(level2), H(level0)+H(level1)+H(level2)+H(level3), H(level0)+H(level1)+H(level2)+H(level3)+H(level4)];
[fmax,ifmax]= max(f);
for i=1:ifmax-1, f(i) = 0.0; end
level0  = level0 -  1;
level1  = level1 -  1;
level2  = level2 -  1;
level3  = level3 -  1;
level4  = level;
y = f;
x = [level0 level1 level2 level3 level4];

for i=1:4,
    s(i)=(y(i)-y(i+1)/(x(i)-x(i+1)));
end

mhdText1 = strrep(mhdText, '.raw', 'Rescale.raw');
mhdFile = [dirInput ImageFile 'Rescale.mhd'];
fp = fopen(mhdFile,'w');
fprintf(fp, '%s',mhdText1);
fclose(fp);
rawFile = [dirInput ImageFile 'Rescale.raw'];
fp  = fopen(rawFile,'wb');
fwrite(fp,int16(ima),'int16');
fclose(fp);

figure(1)
imshow(int16(ima(:,:,60)'),[0 4095])

for k=1:SLI,
    for i=1:ROW,
        for j=1:COL,
            k_ = ima(i,j,k);
            if(k_<level1)
                ienh(i,j,k) = y(1);
            end
            if((k_>=level1) & (k_<level2))
               ienh(i,j,k) = s(2)*(k_-x(2))+y(2);
            end
            if(k_>=level2)
               ienh(i,j,k) = y(4);
            end
        end
    end
end

imax = double(max(max(max(ienh))));
imin = double(min(min(min(ienh))));
ienh_ = (level*(ienh - imin)/(imax - imin));

figure(2)
imshow(ienh_(:,:,60)',[0 4095])

mhdText1 = strrep(mhdText, '.raw', 'Enh.raw');
mhdFile = [dirInput ImageFile 'Enh.mhd'];
fp = fopen(mhdFile,'w');
fprintf(fp, '%s',mhdText1);
fclose(fp);
rawFile = [dirInput ImageFile 'Enh.raw'];
fp  = fopen(rawFile,'wb');
fwrite(fp,int16(ienh_),'int16');
fclose(fp);


        