Quick way to generate the Shepp-Logan phantom

Question

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.

Answer 1

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;

ImageAdjust@Image@Total@MapThread[
    #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 // 
  Image // ImageAdjust

This took half a second on my machine.