###### tags: `影像處理`
# *2022/04/19 影像處理 HW07*
## Coding Part
### 實戰(2)
#### 題目:
#### 程式碼：

#實戰演練(2)
import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft

x=np.array([1,1,4,2])
fft_X=fft(x)
print("Original x=",x)
print("fft_X=",fft_X)

dft_x=np.fft.fft(x,16)
plt.plot(dft_x)
plt.ylabel('Frequency')
plt.show()

#### 輸出結果：
Original x= [1 1 4 2]
fft_X= [ 8.-0.j -3.+1.j 2.-0.j -3.-1.j]

### 實戰(3)
#### 題目:
#### 理想濾波器程式碼：

#理想濾波器
import numpy as np
import cv2
from numpy.fft import fft2,ifft2,fftshift
from matplotlib import pyplot as plt
from google.colab import drive
from google.colab.patches import cv2_imshow
drive.mount('/content/drive')
img = cv2.imread("/content/drive/My Drive/Colab Notebooks/imgs/Barbara.bmp",0)

def spectrum(f):
    F=fft2(f)
    Fshift=fftshift(F)
    mag=20*np.log(np.abs(Fshift)+1)
    mag=mag/mag.max()*255.0
    g=np.uint8(mag)
    return g

def phase_spectrum(f):
    F=fft2(f)
    phase=np.angle(F,deg=True)
    nr,nc=phase.shape[:2]
    for x in range(nr):
        for y in range(nc):
            if phase[x,y]<0:
                phase[x,y]=phase[x,y]+360
            phase[x,y]=int(phase[x,y]*255/360)
        g=np.uint8(np.clip(phase,0,255))
        return g

def frequency_filtering(f,filter,D0):
    nr,nc = f.shape[:2]
    fp = np.zeros([nr,nc])
    for x in range(nr):
        for y in range(nc):
            fp[x,y] = pow(-1,x+y)*f[x,y]
    F = fft2(fp)
    G = F.copy()
    if filter ==1:
        for u in range(nr):
            for v in range(nc):
                dist = np.sqrt((u-nr/2)*(u-nr/2))+((v-nc/2)*(v-nc/2))
                if dist > D0:
                    G[u,v]=0
    elif filter==2:
        for u in range (nr):
            for v in range(nc):
                dist = np.sqrt((u-nr/2)*(u-nr/2))+((v-nc/2)*(v-nc/2))
                if dist <= D0:
                    G[u,v]=0
    gp = ifft2(G)
    gp2 = np.zeros([nr,nc])
    for x in range(nr):
        for y in range(nc):
            gp2[x,y]=round(pow(-1,x+y)*np.real(gp[x,y]),0)
    g = np.uint8(np.clip(gp2,0,255))
    return g

Magnitude_img = spectrum(img)
Phase=phase_spectrum(img)

img1_lpft = frequency_filtering(img,1,50)
Magnitude_img_1=spectrum(img1_lpft)
Phase_1=phase_spectrum(img1_lpft)

img2_lpft = frequency_filtering(img,2,50)
Magnitude_img_2=spectrum(img2_lpft)
Phase_2=phase_spectrum(img2_lpft)

titles = ["Original","Original Magnitude","Original Phase","ideal LowPass","ideal LowPass_Magnitude","ideal LowPass_Phase","ideal HighPass","ideal HighPass_Magnitude","ideal HighPass_Phase"]
images = [img,Magnitude_img,Phase,img1_lpft,Magnitude_img_1,Phase_1,img2_lpft,Magnitude_img_2,Phase_2]
plt.figure(figsize=(20,15))

for i in range(9):
    plt.subplot(3,3,i+1)
    plt.imshow(images[i],'gray')
    plt.title(titles[i],fontsize=15,color='r')

plt.tight_layout()
plt.show()

#### 輸出結果：

#### 高斯濾波器程式碼：

#高斯濾波器
import numpy as np
import cv2
from numpy.fft import fft2,ifft2,fftshift
from matplotlib import pyplot as plt
from google.colab import drive
from google.colab.patches import cv2_imshow
drive.mount('/content/drive')
img = cv2.imread("/content/drive/My Drive/Colab Notebooks/imgs/Barbara.bmp",0)

def spectrum(f):
    F=fft2(f)
    Fshift=fftshift(F)
    mag=20*np.log(np.abs(Fshift)+1)
    mag=mag/mag.max()*255.0
    g=np.uint8(mag)
    return g

def phase_spectrum(f):
    F=fft2(f)
    phase=np.angle(F,deg=True)
    nr,nc=phase.shape[:2]
    for x in range(nr):
        for y in range(nc):
            if phase[x,y]<0:
                phase[x,y]=phase[x,y]+360
            phase[x,y]=int(phase[x,y]*255/360)
        g=np.uint8(np.clip(phase,0,255))
        return g

def frequency_filtering(f,filter,D0):
    nr,nc = f.shape[:2]
    fp = np.zeros([nr,nc])
    for x in range(nr):
        for y in range(nc):
            fp[x,y] = pow(-1,x+y)*f[x,y]
    F = fft2(fp)
    G = F.copy()
    if filter ==3:
        for u in range(nr):
            for v in range(nc):
                dist = np.sqrt((u-nr/2)*(u-nr/2))+((v-nc/2)*(v-nc/2))
                H=np.exp(-(dist*dist)/(2*D0*D0))
                G[u,v]*=H
    elif filter==4:
        for u in range (nr):
            for v in range(nc):
                dist = np.sqrt((u-nr/2)*(u-nr/2))+((v-nc/2)*(v-nc/2))
                H=1-np.exp(-(dist*dist)/(2*D0*D0))
                G[u,v]*=H

    gp = ifft2(G)
    gp2 = np.zeros([nr,nc])
    for x in range(nr):
        for y in range(nc):
            gp2[x,y]=round(pow(-1,x+y)*np.real(gp[x,y]),0)
    g = np.uint8(np.clip(gp2,0,255))
    return g

Magnitude_img = spectrum(img)
Phase=phase_spectrum(img)

img1_glp = frequency_filtering(img,3,50)
Magnitude_img_glp_1=spectrum(img1_glp)
Phase_glp_1=phase_spectrum(img1_glp)

img2_glp = frequency_filtering(img,4,50)
Magnitude_img_glp_2=spectrum(img2_glp)
Phase_glp_2=phase_spectrum(img2_glp)

titles = ["Original","Original Magnitude","Original Phase","Gaussian LowPass","Gaussian LowPass_Magnitude","Gaussian LowPass_Phase","Gaussian HighPass","Gaussian HighPass_Magnitude","Gaussian HighPass_Phase"]
images = [img,Magnitude_img,Phase,img1_glp,Magnitude_img_glp_1,Phase_glp_1,img2_glp,Magnitude_img_glp_2,Phase_glp_2]
plt.figure(figsize=(20,15))

for i in range(9):
    plt.subplot(3,3,i+1)
    plt.imshow(images[i],'gray')
    plt.title(titles[i],fontsize=15,color='g')

plt.tight_layout()
plt.show()

#### 輸出結果:

#### 巴特沃斯濾波器程式碼:

#巴斯特濾波器
import numpy as np
import cv2
from numpy.fft import fft2,ifft2,fftshift
from matplotlib import pyplot as plt
from google.colab import drive
from google.colab.patches import cv2_imshow
drive.mount('/content/drive')
img = cv2.imread("/content/drive/My Drive/Colab Notebooks/imgs/Barbara.bmp",0)

def spectrum(f):
    F=fft2(f)
    Fshift=fftshift(F)
    mag=20*np.log(np.abs(Fshift)+1)
    mag=mag/mag.max()*255.0
    g=np.uint8(mag)
    return g

def phase_spectrum(f):
    F=fft2(f)
    phase=np.angle(F,deg=True)
    nr,nc=phase.shape[:2]
    for x in range(nr):
        for y in range(nc):
            if phase[x,y]<0:
                phase[x,y]=phase[x,y]+360
            phase[x,y]=int(phase[x,y]*255/360)
        g=np.uint8(np.clip(phase,0,255))
        return g

def frequency_filtering(f,filter,D0):
    nr,nc = f.shape[:2]
    fp = np.zeros([nr,nc])
    for x in range(nr):
        for y in range(nc):
            fp[x,y] = pow(-1,x+y)*f[x,y]
    F = fft2(fp)
    G = F.copy()
    if filter ==5:
        for u in range(nr):
            for v in range(nc):
                dist = np.sqrt((u-nr/2)*(u-nr/2))+((v-nc/2)*(v-nc/2))
                H=1.0/(1.0+pow(dist/D0,2))
                G[u,v]*=H
    elif filter==6:
        for u in range (nr):
            for v in range(nc):
                dist = np.sqrt((u-nr/2)*(u-nr/2))+((v-nc/2)*(v-nc/2))
                H=1.0-1.0/(1.0+pow(dist/D0,2))
                G[u,v]*=H

    gp = ifft2(G)
    gp2 = np.zeros([nr,nc])
    for x in range(nr):
        for y in range(nc):
            gp2[x,y]=round(pow(-1,x+y)*np.real(gp[x,y]),0)
    g = np.uint8(np.clip(gp2,0,255))
    return g

Magnitude_img = spectrum(img)
Phase=phase_spectrum(img)

img1_blp = frequency_filtering(img,5,50)
Magnitude_img_blp_1=spectrum(img1_blp)
Phase_blp_1=phase_spectrum(img1_blp)

img2_blp = frequency_filtering(img,5,50)
Magnitude_img_blp_2=spectrum(img2_blp)
Phase_blp_2=phase_spectrum(img2_blp)

titles = ["Original","Original Magnitude","Original Phase","Butter_LowPass","Butter_LowPass_Magnitude","Butter_LowPass_Phase","Butter_HighPass","Butter_HighPass_Magnitude","Butter_HighPass_Phase"]
images = [img,Magnitude_img,Phase,img1_blp,Magnitude_img_blp_1,Phase_blp_1,img2_blp,Magnitude_img_blp_2,Phase_blp_2]
plt.figure(figsize=(20,15))

for i in range(9):
    plt.subplot(3,3,i+1)
    plt.imshow(images[i],'gray')
    plt.title(titles[i],fontsize=15,color='b')

plt.tight_layout()
plt.show()

#### 輸出結果:

## 簡答題：
1. 將能量集中便於使用率波器與觀察，將低頻放在中間，高頻放在四個角落
2. 在空間域內影像對濾波器作捲積運算就等效於在頻率域內影像對濾波器作乘積運算