Skylighty

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.  
  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

Adblocker detected! Please consider disabling it...

We've detected AdBlock Plus or some other adblocking software preventing Pastebin.com from fully loading.

We don't have any obnoxious sound, or popup ads, we actively block these annoying types of ads!

Please add Pastebin.com to your ad blocker whitelist or disable your adblocking software.

×