Skylighty

tomov2

Oct 28th, 2020
792
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. # Tu umieść swój skrypt
  2.  
  3. #Załadowanie niezbędnych bibiotek
  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()
  13.     ax = fig.add_subplot(111)
  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.  
  53. Img = imageio.imread(sFileName)
  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