# tomov2

Oct 28th, 2020
878
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. # Tu umieść swój skrypt
2.
4. import imageio
5. import numpy as np
6. from scipy import ndimage, signal
7. import matplotlib.pyplot as plt
8.
9. # Plot the filter pulse response
10. def plot_response(fs, w, h, title):
11.     "Utility function to plot response functions"
12.     fig = plt.figure()
14.     ax.plot(0.5*fs*w/np.pi, 20*np.log10(np.abs(h)))
15.     ax.set_ylim(-40, 5)
16.     ax.set_xlim(0, 0.5*fs)
17.     ax.grid(True)
18.     ax.set_xlabel('Frequency (Hz)')
19.     ax.set_ylabel('Gain (dB)')
20.     ax.set_title(title)
21.
22.
23. # Band-pass filter (BPF)
24. def butter_bandpass(lowcut, highcut, fs, order=5):
25.     nyq = 0.5 * fs
26.     low = lowcut / nyq
27.     high = highcut / nyq
28.     b,a = signal.butter(order, [low, high], btype='band')
29.     return b, a
30.
31. # Linear filtering through vector
32. def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
33.     b, a = butter_bandpass(lowcut, highcut, fs, order=order)
34.     y = signal.lfilter(b, a, data)
35.     return y
36.
37. # Forward-backward filter
38. def butter_linear(data, lowcut, highcut, fs, order=5):
39.     b, a = butter_bandpass(lowcut, highcut, fs, order=order)
40.     y = signal.filtfilt(b, a, data)
41.     return y
42.
43.
44. ######################################
45. #         Sekcja konfiguracji        #
46. ######################################
47.
48. #Nazwa pliku z sinogramem do wczytania
49. sFileName = "Sinogram1.png"
50. #Nazwa pliku z przekrojem do zapisania
51. sOutFileName = "Tomo1.png"
52.
54. (AngleSteps, ImgWidth) = Img.shape
55.
56. RecImg = np.zeros((ImgWidth, ImgWidth))
57. c = np.zeros((ImgWidth, ImgWidth))
58.
59. #Tu wstaw swój kod
60. for AngleId in range(0, AngleSteps):
61.     print("{}/{}".format(AngleId, AngleSteps))
62.
63. # High-pass filter (HPF)
64.     # Nyquist rate.
65.     nyq_rate = 48000 / 2
66.     # Width of the roll-off region.
67.     width = 150 / nyq_rate
68.     # Attenuation in the stop band.
69.     ripple_db = 15.0
70.     num_of_taps, beta = signal.kaiserord(ripple_db, width)
71.     if num_of_taps % 2 == 0:
72.         num_of_taps = num_of_taps + 1
73.     # Cut-off frequency.
74.     cutoff_hz = 350.0
75.     # Estimate the filter coefficients.
76.     taps = signal.firwin(num_of_taps, cutoff_hz / nyq_rate, window=('kaiser', beta), pass_zero=False)
77.
78.     line = Img[AngleId, :]
79.     #filtered = np.convolve(line, taps, mode='same')
80.     filtered = butter_linear(line, 400.0, 8000.0, 48000.0, order=2)
81.
82.     for i in range(0, 200):
83.         c[i, :] = filtered
84.     RotatedImage = ndimage.rotate(c, -360.0 * AngleId / AngleSteps, reshape=False)
85.     RecImg += RotatedImage
86.
87. n = sum(sum(RecImg))
88. plt.imshow(RecImg/n, cmap='gray')
89. plt.show()
90.
91. #Zapis wyniku do pliku
92. imageio.imwrite(sOutFileName, RecImg)
93.
94. # Wniosek - metodą prób i błędów należy dobrać parametry filtrów,
95. # jego częstotliwości odcięcia, natomiast wciąż, nawet przy filtrach
96. # wyższego rzędu nie jesteśmy w stanie ostatnie pozbyć się artefaktów rekonstrukcji
97. #
98. # Podobne rezultaty otrzymamy za równo dla filtra górno-przepustowego
99. # jak i pasmowo-przepustowego.
RAW Paste Data