Shepp Logan Phantom

shepp_logan

function shepp_logan(N,M)
  P = zeros(N,M)

  x = linspace(-1,1,N)
  y = (-1)*linspace(-1,1,M)'

  s = 0.4

  P += 1.0 * ( (x/0.69).^2 .+ (y/0.92).^2 .< 1)
  P += - s * ( (x/0.6624).^2 .+ ((y+0.0184)/0.874).^2 .< 1 )
  P += - 0.2 * ( ((cos(-18.0/360*2*pi).*(x-0.22).+sin(-18.0/360*2*pi)*y)/0.11).^2 .+
      ( (sin(-18.0/360*2*pi)*(x-0.22).-cos(-18.0/360*2*pi)*y)/0.31).^2 .< 1 )
  P += - 0.2 * ( ( (cos( 18.0/360*2*pi).*(x +0.22).+sin( 18.0/360*2*pi)*y)/0.16).^2 .+
      ( (sin( 18.0/360*2*pi).*(x +0.22).-cos( 18.0/360*2*pi)*y)/0.41).^2 .< 1 )
  P += 0.1 * ( (x/0.21).^2 .+ ((y-0.35)/0.25).^2 .< 1)
  P += 0.1 * ( (x/0.046).^2 .+ ((y-0.1)/0.046).^2 .< 1 )
  P += 0.1 * ( (x/0.046).^2 .+ ((y+0.1)/0.046).^2 .< 1 )
  P += 0.1 * ( ((x+0.08)/0.046).^2 .+ ((y+0.605)/0.023).^2 .< 1 )
  P += 0.1 * ( (x/0.023).^2 .+ ((y+0.606)/0.023).^2 .< 1 )
  P += 0.1 * ( ((x-0.06)/0.023).^2 .+ ((y+0.605)/0.046).^2 .< 1 )

  return P'
end

function shepp_logan(N)
  shepp_logan(N,N)
end