Skip to content

Instantly share code, notes, and snippets.

@tknopp
Last active March 8, 2024 07:34
Show Gist options
  • Save tknopp/7872297 to your computer and use it in GitHub Desktop.
Save tknopp/7872297 to your computer and use it in GitHub Desktop.
Shepp Logan Phantom
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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment