Help understanding the algorithms used in this code.

Hello! I came across this volumetric imaging code and I'm trying to better understand how the algorithms work. Could someone help explain what each algorithm is doing?

clear;
vol=zeros(98,98,111);
fp = fopen('pigSmall.img','r');
for z=1:110;
z
for i=1:97
for j=1:97
rho=fread(fp,1,'char');
vol(i,j,z)=rho;
end
end
end

img=zeros(98,98);
minth=20;
P=eye(4);
P(3,3)=0;
vox=zeros(4,1);
pos=zeros(4,1);
for i=1:97
for j=1:97
for z=1:110
vox=transpose([i,j,z,1]);
pos=P*vox;
rho=vol(i,j,z);
if pos(1,1) > 0 & pos(1,1) < 98 & pos(2,1) > 0 & pos(2,1) < 98
if rho > img(pos(1,1),pos(2,1)) & rho > minth
img(pos(1,1),pos(2,1)) = rho;
end
end
end
end
end
minint=min(min(img));
img=img-minint;
maxint = max(max(img));
img=img/maxint*64;
colormap(gray)
image(img)

Since mathworks wont let me attach the pigSmall.img file this is the output of the code:

Answer 1

[jonas]

I've tried to decipher this rather cryptic code. It appears to essentially read a binary file of size 98 x 98 x 111 and store the maximum value along the depth dimension in a 98x98 matrix. The code is extremely inefficient, to the point where it seems as if it's almost intentionally trying to mislead a reader.

clear;

% Preallocate a volume matrix size 98x98x111
vol=zeros(98,98,111);

% Open img file, not sure what the content is
fp = fopen('pigSmall.img','r');

% Reads each element of the binary img file and stores it in the vol matrix...
% which then becomes a binary matrix with 111 layers
for z=1:110;
    z
    for i=1:97
        for j=1:97
            rho=fread(fp,1,'char');
            vol(i,j,z)=rho;
        end
    end
end

% Pre allocate a single layer called img
img=zeros(98,98);

% Probably a "min threshold" value
minth=20;

% identity matrix 
P=eye(4);

% Set id 3;3 to zero for some reason
P(3,3)=0;

% preallocate some vectors
vox=zeros(4,1);
pos=zeros(4,1);

% Again, loop over all elements of vol
for i=1:97
    for j=1:97
        for z=1:110
            % Same as vox=[i;j;z;1]
            vox=transpose([i,j,z,1]);
            % Next line outputs pos=[i;j;0;1]
            pos=P*vox;
            % Grab each value in vol iteratively
            rho=vol(i,j,z);
            % I believe this condition is always true and thus useless line
            % of code
            if pos(1,1) > 0 & pos(1,1) < 98 & pos(2,1) > 0 & pos(2,1) < 98
                % Store value (i,j,z) in img if it is larger than current
                % value and the threshold value
                if rho > img(pos(1,1),pos(2,1)) & rho > minth
                    img(pos(1,1),pos(2,1)) = rho;
                end
            end
        end
    end
end
% minimum intensity image
minint=min(min(img));
% subtract all values by min intensity
img=img-minint;
% Max intensity
maxint = max(max(img));
% probably normalizing the image
img=img/maxint*64;
colormap(gray)
image(img)

[Westin Messer]

Thanks so much Jonas!