tomov2

# Tu umieść swój skrypt

#Załadowanie niezbędnych bibiotek
import imageio
import numpy as np
from scipy import ndimage, signal
import matplotlib.pyplot as plt

# Plot the filter pulse response
def plot_response(fs, w, h, title):
    "Utility function to plot response functions"
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.plot(0.5*fs*w/np.pi, 20*np.log10(np.abs(h)))
    ax.set_ylim(-40, 5)
    ax.set_xlim(0, 0.5*fs)
    ax.grid(True)
    ax.set_xlabel('Frequency (Hz)')
    ax.set_ylabel('Gain (dB)')
    ax.set_title(title)


# Band-pass filter (BPF)
def butter_bandpass(lowcut, highcut, fs, order=5):
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq
    b,a = signal.butter(order, [low, high], btype='band')
    return b, a

# Linear filtering through vector
def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
    b, a = butter_bandpass(lowcut, highcut, fs, order=order)
    y = signal.lfilter(b, a, data)
    return y

# Forward-backward filter
def butter_linear(data, lowcut, highcut, fs, order=5):
    b, a = butter_bandpass(lowcut, highcut, fs, order=order)
    y = signal.filtfilt(b, a, data)
    return y


######################################
#         Sekcja konfiguracji        #
######################################

#Nazwa pliku z sinogramem do wczytania
sFileName = "Sinogram1.png"
#Nazwa pliku z przekrojem do zapisania
sOutFileName = "Tomo1.png"

Img = imageio.imread(sFileName)
(AngleSteps, ImgWidth) = Img.shape

RecImg = np.zeros((ImgWidth, ImgWidth))
c = np.zeros((ImgWidth, ImgWidth))

#Tu wstaw swój kod
for AngleId in range(0, AngleSteps):
    print("{}/{}".format(AngleId, AngleSteps))

# High-pass filter (HPF)
    # Nyquist rate.
    nyq_rate = 48000 / 2
    # Width of the roll-off region.
    width = 150 / nyq_rate
    # Attenuation in the stop band.
    ripple_db = 15.0
    num_of_taps, beta = signal.kaiserord(ripple_db, width)
    if num_of_taps % 2 == 0:
        num_of_taps = num_of_taps + 1
    # Cut-off frequency.
    cutoff_hz = 350.0
    # Estimate the filter coefficients.
    taps = signal.firwin(num_of_taps, cutoff_hz / nyq_rate, window=('kaiser', beta), pass_zero=False)

    line = Img[AngleId, :]
    #filtered = np.convolve(line, taps, mode='same')
    filtered = butter_linear(line, 400.0, 8000.0, 48000.0, order=2)

    for i in range(0, 200):
        c[i, :] = filtered
    RotatedImage = ndimage.rotate(c, -360.0 * AngleId / AngleSteps, reshape=False)
    RecImg += RotatedImage

n = sum(sum(RecImg))
plt.imshow(RecImg/n, cmap='gray')
plt.show()

#Zapis wyniku do pliku
imageio.imwrite(sOutFileName, RecImg)

# Wniosek - metodą prób i błędów należy dobrać parametry filtrów,
# jego częstotliwości odcięcia, natomiast wciąż, nawet przy filtrach
# wyższego rzędu nie jesteśmy w stanie ostatnie pozbyć się artefaktów rekonstrukcji
#
# Podobne rezultaty otrzymamy za równo dla filtra górno-przepustowego
# jak i pasmowo-przepustowego.