Quick way to generate the Shepp-Logan phantom

I have some very old code for generating a Shepp-Logan phantom:

SheppLoganData = {
{{0, 0}, {0.92, 0.69}, 90°, 2}, {{0, -0.0184}, {0.874, 0.6624}, 90°, -0.9},
{{0.22, 0}, {0.31, 0.11}, 72°, -0.1}, {{-0.22, 0}, {0.41, 0.16}, 108°, -0.1},
{{0, 0.35}, {0.25, 0.21}, 90°, 0.3}, {{0, 0.1}, {0.046, 0.046}, 0, 0.3},
{{0, -0.1}, {0.046, 0.046}, 0, 0.3}, {{-0.08, -0.605}, {0.046, 0.023}, 0, 0.3},
{{0.06, -0.605}, {0.046, 0.023}, 90°, 0.3}, {{0, -0.605}, {0.023, 0.023}, 0, 0.3}
};

DensityPlot[Total[#4 Boole[Norm[({x, y} - #1).RotationMatrix[#3]/#2]^2 < 1] & @@@
SheppLoganData], {x, -1, 1}, {y, -1, 1}, ColorFunction -> GrayLevel,
Frame -> False, Mesh -> False, PlotPoints -> 500]


(The previous version used RotationMatrix2D from the Rotations package, but I have updated it to use the current built-in.)

This is very slow, however. Image is apparently what is used today to represent images in Mathematica, so is there any way to adapt my old code to produce an Image of the Shepp-Logan phantom for a specified size (e.g. a 1024×1024 image)? Thanks for any ideas.

• What do the numbers under "Gray Level" mean in the table in Wikipedia? – Szabolcs Jun 30 '16 at 10:08
• @Szabolcs, I had a look at the original paper; it seems those values get added up whenever a given point is within any of those ten ellipses. – J. M.'s torpor Jun 30 '16 at 10:54
• You could also do the following: take the symbolic expression you generate for this: Total[#4 Boole[norm[({x, y} - #1).RotationMatrix[#3]/#2]^2 < 1] & @@@ SheppLoganData]. It may be better to define norm = Sqrt[#.#]& instead of Norm to save some unneeded Abs calls. Take the resulting expression, and copy-paste it inside of cf = Compile[{x,y}, ...]. Then use this function as Array[cf, {512, 512}, {{-1, 1}, {-1, 1}}] // Transpose // Reverse // Image // ImageAdjust. It's faster than my answer, and fairly straightforward too. – Szabolcs Jun 30 '16 at 11:13
• If instead of compiling, you could take the time to rewrite this in vectorizable form, it would be faster still. – Szabolcs Jun 30 '16 at 11:14
• @Szabolcs, many thanks for your contribution. Is there any way Melania Black original code could run faster while retaining "DensityPlot"? – Dean May 1 '20 at 0:06

Using Graphics

Using the data supplied by yourself, you could do something along the lines of:

grs = Graphics[{White, Rotate[Disk[#1, #2], #3]}, Background -> Black,
PlotRange -> {{-1, 1}, {-1, 1}}] & @@@ SheppLoganData;

#2 ImageData[
ColorConvert[
Rasterize[Style[#, Antialiasing -> False], "Image",
ImageSize -> 512],
"Grayscale"],
"Real"] &,
{
grs,
SheppLoganData[[All, -1]]
}
]


It's not fast, but neither is it unbearably slow.

Alternative method

Take the expression you wrote yourself, and use it in a compiled function. I will make a small change and replace Norm by norm = Sqrt[#.#]& to get rid of some Abs calls.

Total[#4 Boole[norm[({x, y} - #1).RotationMatrix[#3]/#2]^2 < 1] & @@@
SheppLoganData] // Simplify

0.3 Boole[22.6757 x^2 + 16. (-0.35 + y)^2 < 1] +
0.3 Boole[472.59 x^2 + 472.59 (-0.1 + y)^2 < 1] -
0.1 Boole[
0.0351981 + 1. x^2 + 0.123325 y + 0.228447 y^2 <
0.44 x + 0.560566 x y] -
0.1 Boole[
0.0205452 + 0.44 x + 1. x^2 + 0.119275 y + 0.542158 x y +
0.253783 y^2 < 0] + 2 Boole[2.1004 x^2 + 1.18147 y^2 < 1] -
0.9 Boole[2.27908 x^2 + 1.30911 (0.0184 + y)^2 < 1] +
0.3 Boole[472.59 x^2 + 472.59 (0.1 + y)^2 < 1] +
0.3 Boole[1890.36 (-0.06 + x)^2 + 472.59 (0.605 + y)^2 < 1] +
0.3 Boole[1890.36 x^2 + 1890.36 (0.605 + y)^2 < 1] +
0.3 Boole[472.59 (0.08 + x)^2 + 1890.36 (0.605 + y)^2 < 1]


In order to avoid having to deal with the subtleties of code generation, I just copy this and paste it into a compiled function:

cf = Compile[{x, y},
0.3 Boole[22.675736961451246 x^2 + 16. (-0.35 + y)^2 < 1] +
0.3 Boole[
472.5897920604915 x^2 + 472.5897920604915 (-0.1 + y)^2 <
1] - 0.1 Boole[
0.03519805880091542 + 1. x^2 + 0.1233245922675553 y +
0.22844664024299033 y^2 <
0.44 x + 0.5605663284888875 x y] -
0.1 Boole[
0.02054517798250305 + 0.4400000000000001 x + 1. x^2 +
0.11927486660565205 y + 0.5421584845711455 x y +
0.2537828638062259 y^2 < 0] +
2 Boole[2.100399075824407 x^2 + 1.1814744801512287 y^2 < 1] -
0.9 Boole[
2.2790788583164137 x^2 + 1.3091129973974835 (0.0184 + y)^2 <
1] + 0.3 Boole[
472.5897920604915 x^2 + 472.5897920604915 (0.1 + y)^2 < 1] +
0.3 Boole[
1890.359168241966 (-0.06 + x)^2 +
472.5897920604915 (0.605 + y)^2 < 1] +
0.3 Boole[
1890.359168241966 x^2 + 1890.359168241966 (0.605 + y)^2 <
1] + 0.3 Boole[
472.5897920604915 (0.08 + x)^2 +
1890.359168241966 (0.605 + y)^2 < 1]
];


We can verify that it will make no callbacks to the main evaluator.

Now instead of using ListDensityPlot, we use it with Array:

Array[cf, {512, 512}, {{-1, 1}, {-1, 1}}] // Transpose // Reverse //