MATLAB code

Water sphere Geometry Writer Code:

%Geometry Writer

%This writer produces a spherical water phantom 160mm in diameter

close all

clear all

Geometry=zeros(200,200,200);

Radius=80;

for k=1:200

for j=1:200

for i=1:200

Distance=((i-100.5)^2+(j-100.5)^2+(k-100.5)^2).^.5;

if Distance <= Radius

Geometry(i,j,k)=1;

end

end

end

end



save 'Geom.mat' Geometry;



% for k=1:200

% imagesc(squeeze(Geometry(:,:,k)))

% pause(0.5)

% end

%Now write a Geometry file

fid = fopen('GEOMETRY.IN','w')

fprintf(fid,'# 1mm Water Spherical Phantom with 1mm thick Bone
shell.\n');

fprintf(fid,'# Total number of material typed\n');

fprintf(fid,'3\n');

fprintf(fid,'# names asigned to each material in the material input
file.\n');

fprintf(fid,'# NB: these must be in upper case.\n');

fprintf(fid,'AIRVOID \n');

fprintf(fid,'WATER \n');

fprintf(fid,'BONE \n');

fprintf(fid,'# integer code asigned to each material\n');

fprintf(fid,'0\n');

fprintf(fid,'1\n');

fprintf(fid,'2\n');

fprintf(fid,'# number of voxels along Y\n');

fprintf(fid,'200\n');

fprintf(fid,'# number of voxels along X\n');

fprintf(fid,'200\n');

fprintf(fid,'# number of voxels along Z\n');

fprintf(fid,'200\n');

fprintf(fid,'# voxel dimension in mm (X,Y,Z)\n');

fprintf(fid,'1.000000\n');

fprintf(fid,'1.000000\n');

fprintf(fid,'1.000000\n');

fprintf(fid,'# anatomical data set, sweep through i, j then k.\n');

for k=1:200

for j=1:200

for i=1:200

fprintf(fid,'%d',Geometry(i,j,k));

end

fprintf(fid,'\n');

end

end

fclose(fid)

We start with the trapezoidal prism.

We need to set the value of the array to "1" for all the points that are inside the prism.
From the drawing we can see that we need to scan for the points that have their z
coordinate between z0 (the height of the prism floor above the skull's bottom) and z0+H,
H being the prism height.

The same holds for the x-coordinate. We need to scan from x0 (the distance of the
rightmost prism wall) to x0+L, L being the prism's length.




So far, pretty easy ­ we have two standard loops: z from z0 to z0+H and x from x0 to x0+L.
Now let's look at side view:




We see that we have to scan a variable length in the y-direction. This length is the
distance from the front face of the prism (the tilted red line) and the back wall, which is
located at y0+W1 and since it is a function of the height of the point from the floor of the
prism we denote the distance d(z-z0).
To find out the explicit expression for d(z-z0) look at the following drawing:




Note that the triangles ADE and ABC are similar (they have the same angles), thus the
ratios of their respective sides is the same.

Also note that AC=KG=H (the height of the prism).

Therefore:
DE CE
=
BA AC


DE =
(BA)(CE )
AC
Now,
CE = z - z 0
BA = W1 - W2
AC = H

Hence:
(W1 - W2 )(z - z 0 )
DE =
H

So the total distance from the tilted face to the back wall is:

d ( z ) = DF = W2 +
(W1 - W2 )(z - z0 )
H

And:
(W - W2 )(z - z 0 ) (z - z 0 )
y i = [ y 0 + W1 ] - d ( z ) = y 0 + W1 - W2 - 1 = y 0 + (W1 - W2 )1 - H
H

And therefore along the y coordinate we scan along the interval from y ( y i , y 0 + W1 )
Note that when z=z0 (the bottom face), we get y i = y 0 + W1 - W2 , and when z=z0+H (the
top face) we get y i = y 0
Also note that when W1=W2, the scan interval is constant from y0 to y0+W1

So the basic algorithm is very simple:

· Set the initial prism parameters: x0,y0,z0 and H,L,W1,W2, make sure that W1>W2

For x=x0 to x0+L

· For z=z0 to z0+H

(z - z 0 )
o Calculate y i = ceil y 0 + (W1 - W2 )1 -
H

o For y=yi to y0+W1

· Geometry(x,y,z)=0
o Next y

· Next z

Next x
Few Notes: since x,y,z represent the indices of the grid point as well as their spatial
locations, make sure that all the initial values are integers.

Because of the same reason as above, you want to use the ceil function to round the initial
point yi to the higher integer (to start at yi =7 rather at yi =6.2).


For the triangular prism (shape a): just note that this is a special case of the trapezoidal
prism with W2=0, and in our case the shape c is flipped:




And we get:
DE BC
=
AE AC

Since:
AC = H
BC = W1 - W2
AE = H - z + z 0
We obtain:

(W1 - W2 )(H - z + z 0 )
DE =
H
(W1 - W2 )(H - z + z 0 )
d ( z ) = W2 + DE = W2 +
H
(W - W2 )(H - z + z 0 ) (H - z + z 0 )
y i = W1 - d ( z ) = W1 - W2 - 1 = (W1 - W2 )1 -
H H

And we can now stick this formula in the algorithm. Since shape c is a special case of a
trapezoidal prism in which W2=0, then when z=z0+H (the final z value along the z axis)
we get:


(H + z 0 - z )
y i = ceil y 0 + (W1 - W2 )1 -
H
(H + z 0 - z 0 - H )
y i = ceil y 0 + W1 1 -
H
0
y i = ceil y 0 + W1 1 - = y 0 + W1
H


But this is also the final value for the scan along the y-axis.