줄리아, 매트랩, 파이썬에서 Shepp-Logan 팬텀 쓰는 법

줄리아, 매트랩, 파이썬에서 Shepp-Logan 팬텀 쓰는 법

매트랩

phantom으로 Shepp-Logan 팬텀을 생성할 수 있다. 아무런 변수도 입력하지 않으면 $256 \times 256$의 배열이 만들어진다. 첫번째 변수에는 팬텀의 종류, 두번째 변수에는 크기를 입력하여 임의의 팬텀을 생성할 수 있다. 아래는 예제코드와 실행결과이다.

%팬텀 생성
p1 = phantom();
p2 = phantom('Modified Shepp-Logan',2^10);

%그림 그리기
figure()

%첫번째 그림 p1
subplot(1,2,1)
imagesc(p1)
colorbar
title('p1')

%두번째 그림 p2
subplot(1,2,2)
imagesc(p2)
colorbar
title('p2')

환경

Matlab.png

파이썬

사이킷이미지 패키지 skimagedata 모듈에서 shepp_logan_phantom함수를 쓰면 된다. 변수를 입력하지 않으면 $400 \times 400$의 팬텀을 생성한다. 아래는 예제코드와 실행결과이다.

import matplotlib.pyplot as plt

from skimage.data import shepp_logan_phantom

#Phantom 생성
P = shepp_logan_phantom() #default size=(400,400)

#그림 출력
plt.imshow(P)
plt.show()

Python.png

환경

줄리아

지원하지 않는다. 직접 코드를 짜서 실행해보면 다음과 같이 작동한다. 가독성을 위해 최하단에 phantom의 전체 코드를 올려두었다.

julia> p=phantom(100,2)
100×100 Array{Float64,2}:

julia> using Plots

julia> fig=heatmap(p, yflip=true)

julia> savefig(fig,"Julia.png")

Julia.png

전체코드

타원에 대한 값은 위키피디아와 매트랩의 phantom 코드를 뜯어봐서 참고하였다.

function phantom(N::Int64, n::Int64) #N:이미지 크기, n:팬텀의 타입
    #Shepp-Logan phantom n=1
    #            A     a       b        x0     y0        phi     
    #            k1    k2      k3       k4     k5        k6          
    # ---------------------------------------------------------
    sheplog = [ 1.0   0.69    0.92     0.0    0.0       0.0   ;   
               -0.98  0.6624  0.8740   0.0   -0.0184    0.0   ;
                0.01  0.2100  0.2500   0.0    0.35      0.0   ;
                0.01  0.0460  0.0460   0.0    0.1       0.0   ;
                0.01  0.0460  0.0460   0.0   -0.1       0.0   ;
                0.01  0.0460  0.0230  -0.08  -0.605     0.0   ;
                0.01  0.0230  0.0230   0.0   -0.606     0.0   ;
                0.01  0.0230  0.0460   0.06  -0.605     0.0   ;
               -0.02  0.1100  0.3100   0.22   0.0      -0.314 ;
               -0.02  0.1600  0.4100  -0.22   0.0       0.314 ]
  
    #modified-Shepp_Logan phantom n=2
    #              A    a       b        x0     y0       phi     
    #              k1   k2      k3       k4     k5       k6          
    #----------------------------------------------------------
    modishep = [  1.0  0.69    0.92     0.0    0.0      0.0   ;   
                 -0.8  0.6624  0.8740   0.0   -0.0184   0.0   ;
                  0.1  0.2100  0.2500   0.0    0.35     0.0   ;
                  0.1  0.0460  0.0460   0.0    0.1      0.0   ;
                  0.1  0.0460  0.0460   0.0   -0.1      0.0   ;
                  0.1  0.0460  0.0230  -0.08  -0.605    0.0   ;
                  0.1  0.0230  0.0230   0.0   -0.606    0.0   ;
                  0.1  0.0230  0.0460   0.06  -0.605    0.0   ;
                 -0.2  0.1100  0.3100   0.22   0.0     -0.314 ;
                 -0.2  0.1600  0.4100  -0.22   0.0      0.314 ]
  
    if n==1
      sl=sheplog
    elseif n==2
      sl=modishep
    end
  
    P = zeros(Float64, N, N)
    x = LinRange(-1,1,N)
    y = - x
  
    for k=1:10
  
      if k < 9 #회전이 없는 타원
        for i=1:N, j=1:N #i, j 헷갈리지 않게 주의
          if ( (x[j]-sl[k,4])/sl[k,2] )^2 + ( (y[i]-sl[k,5])/sl[k,3] )^2 < 1
            P[i, j]+= sl[k,1]
          end
        end
  
      else #회전이 있는 타원
        for i=1:N, j=1:N
          if ( (cos(sl[k,6])*(x[j]-sl[k,4])+sin(sl[k,6])*(y[i]-sl[k,5]))/sl[k,2] )^2+( (sin(sl[k,6])*(x[j]-sl[k,4])-cos(sl[k,6])*(y[i]-sl[k,5]))/sl[k,3] )^2< 1
            P[i, j]+= sl[k,1]
          end
        end
      end
  
    end
  
    return P
  
 end

환경

댓글