Shepp Logan Phantom

[tknopp]

Hi Tim,

I am not sure whether this package is the best fit for this, but in the following gist their is an implementation of the shepp logan phantom:

https://gist.github.com/tknopp/7872297

I have used this code in Matlab (without imaging toolbox...) and Python when working in the medical image reconstruction field (MRI, CT,...).

[timholy]

That would be a fine thing to have. Do you want to turn it into a PR? It's your own code, right, no licensing concerns?

[tknopp]

I can make a PR of that. Wanted to first get some feedback if you think this is the right place for this.

I am not entirly sure where the original Matlab code is from as it was used by different people at the institute I did my PHD (not what you wanted to here from me...). But the code is actually only about drawing 10 ellipses and it is currently written uggly anywhere. I will therefore pick up the ellipses table from the original paper (http://stat.wharton.upenn.edu/~shepp/publications/33.pdf) rewrite the code to go rowwise through the ellipse-table and then this should be safe from my perspective.

[timholy]

Doing it from the paper is definitely safe. Thanks!

[JeffFessler]

This URL is the first hit for a google search for "julia shepp logan phantom", so in case it is useful to others who end up here, I want to mention that there are now more versions of the Shepp Logan phantom with methods for making images and sinograms and spectra, with Documenter/Literate examples here: https://github.com/JuliaImageRecon/ImagePhantoms.jl

[tknopp]

This is very nice Jeff, Shepp Logan moved into TestImages.jl what IMHO does not fit so well since it now does not output Float64 anymore but Gray. In the long run probably the best would be to remove it from Images.jl and TestImages.jl an cross link ImagePhantoms.jl in both docus.

[johnnychen94]

Shepp Logan moved into TestImages.jl what IMHO does not fit so well since it now does not output Float64 anymore but Gray

This can be supported very easily by introducing shepp_logan(Float64, sz...) usage. https://github.com/JuliaImages/TestImages.jl/issues/138

[tknopp]

its not a matter of supporting it but more a matter of what people that actually use the phantom expect. In the area of tomographic image reconstruction (where most people using this phantom come from) using a floating point number is much more natural.

[johnnychen94]

The whole philosophy of JuliaImages is to bring people from the raw numerical world into a colorant world. If the algorithm was designed to be generic enough, there should be no difference at all whether you're processing Array{Float64, 2} or Array{Gray{Float64}, 2}. The thing that really matters is to design generic functions.

It's understandable that people from other image processing fields and frameworks enjoy their good old raw numerical world. I think this is mainly because JuliaImages documentation doesn't do a good job educating this workflow and best practice and people haven't yet enjoyed this type system. I started some best-practices draft https://github.com/JuliaImages/juliaimages.github.io/pull/214 so if you're interested please check and/or leave some comments there.

[timholy]

Essentially Gray{Float64} is just a "self-documenting type": it tells you that this particular Float64 is a grayscale intensity rather than, e.g., the red channel of an RGB value. Other than that it should be virtually identical to Float64. The people in the world of tomographic imaging should appreciate that. What is an Array{Float64} of size (200, 250, 3): is it a 2d array of RGB values or a 3d grayscale array that happens to have only 3 slices in it? Likewise, should a Matrix{Float64} be displayed as a numeric array or as an image? The world is full of mysteries.

But in Julia:

[JeffFessler]

This discussion is very helpful for me. It reminds me that I have attempted to design ImagePhantoms in a way that supports Unitful arrays, and I should make that more apparent in the Literate examples. The physically appropriate units for Shepp Logan in the CT setting are something like 1/cm for attenuation coefficients, and I should illustrate (and probably encourage) that. Colorants make sense for image processing, whereas in MRI a (200,250,3) array might be (spin density, T1, T2) each with its own physical units. AFAIK the Colorant type system is not designed to support tuples having different units. I'm slowly but systematically trying to update all my code to support units because having units also is a step in the direction of a "self-documenting type".

[timholy]

Interesting! There's an intent to relax the requirement on the T in Gray{T}, and that might accommodate unitful values, if that would help you.

[JeffFessler]

relax the requirement on the T in Gray{T}

That sounds like it could be a useful extension. But to be honest, the term Gray feels a bit unnatural for an image consisting of physical units (like density in g/cm^3) that might be eventually displayed using any possible colormap. It is nice that RGB and HSV have a self-documenting description of their contents in the type name, whereas Gray feels less descriptive, i.e., just about as general as Number.

Gray{Unitful{m}} could be helpful, because now we know it involves length units, but more descriptive would be Gray{Unitful{m}, :depth} where the tag :depth tells us it is a depth map. I'm not really proposing such an extension because probably the array name will be depth and that is a good clue to the reader of the code what it means. But by that logic an array of Float32 named depth is nearly equally descriptive and we're back to numbers :)

I have great admiration for all the work in ImageProcessing.jl to carefully distinguish Gray from colors, but I'll guess that it is probably unlikely that the image formation community will adapt Gray as the generic type, and not just because we are used to using Floats :)

[timholy]

Agreed, it doesn't seem to be a big help in the case of unitful values.

[tknopp]

I agree with Jeff here.

My take here is: The mapping from a physical quantity to a unitless color representation is an important step that should be represented in the types being used. Making everything a color does not help in this respect. But as I said, In the Image processing package family, using a colored representation is perfect (did not want to argue against that). My point was that the Shepp Logan phantom would suite more into the package where no color is put around it. And apparently this goes in line how ImagePhantoms is designed.