clear
clc

level = 255;

%%INPUTS:
%%Input directory
dirInput  = '..../'; 
%%Output directory
dirOutput = '..../';
%%Input volume name
ImageFile1 = '...'; %%SEGMENTED IMAGE
ImageFile2 = '...'; %%GROUND TRUTH IMAGE

mhdFile1 = [dirInput ImageFile1 '.mhd'];
rawFile1 = [dirInput ImageFile1 '.raw'];
mhdFile2 = [dirInput ImageFile2 '.mhd'];
rawFile2 = [dirInput ImageFile2 '.raw'];
mhdText = fileread(mhdFile1);
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(rawFile1,'rb');
ima1 = fread(fp,ROW*COL*SLI,'uint8');
fclose(fp);
ima1 = reshape(ima1, [ROW COL SLI]);

fp  = fopen(rawFile2,'rb');
ima2 = fread(fp,ROW*COL*SLI,'uint8');
fclose(fp);
ima2 = reshape(ima2, [ROW COL SLI]);

cont1 = 0;
cont2 = 0;
cont3 = 0;

for k=1:SLI,
    for i=1:ROW,
        for j=1:COL,
            a = ima1(i,j,k);
            b = ima2(i,j,k);
            out1 = 0;
            out2 = 0;
            if(a ~= 0)
                out1 = 1;
                cont1 = cont1 + 1;
            else
                out1 = 0;
            end
            if(b ~= 0)
                out2 = 1;
                cont2 = cont2 + 1;
            else
                out2 = 0;
            end
            if((out1==1)  & (out2==1))
                cont3 = cont3 + 1;
            end
        end
    end
end


Dice_ = (2.0*cont3)/abs(cont1 + cont2);
fprintf('The Dice Coefficient is: %6.4f or %5.2f%%\n',Dice_,100*Dice_);


        