# 第四堂社課 | 電研X天文聯合社課
## 卷積&去卷積&天文影像處理

---

## 時（空）域、頻域（頻譜圖）
:::info
一張照片的兩個面向
:::

我們平常看到的像素照片，每一個像素都記錄了RGB的值 (相機記錄到的光能)
將RGB值作為Y軸畫成折線圖後，即是時（空）域，他具有這樣一個特徵：
* 時（空）域的訊號是雜亂的

而為了使電腦可以更好的理解這個無規律的訊號，我們需要對照片進行一些處理，把它表示成 "多個有規律的訊號的疊加"，此過程即是傅立葉轉換。

* 傅立葉轉換(Fourier Transfrom)（公式如圖一）
任何一個有規律的函數，都可以表示成其他有規律的函數的疊加 (圖二 ａb)
也就是說透過把一堆有規律的函數相疊加
就能得出一個更複雜且有規律的函數

$$f(x)=\frac{a_0}{2}+\sum^{\infty}_{n=1}(a_ncos{\frac{2\pi nx}{t}}+b_nsin{\frac{2\pi nx}{t}})$$

$$a_0=\frac{1}{\frac{t}{2}}\int^{\frac{t}{2}}_{-\frac{t}{2}}f(x)dx$$

$$a_0=\frac{1}{\frac{t}{2}}\int^{\frac{t}{2}}_{-\frac{t}{2}}\cos\frac{2\pi nx}{t}dx$$

$$a_0=\frac{1}{\frac{t}{2}}\int^{\frac{t}{2}}_{-\frac{t}{2}}sin\frac{2\pi nx}{t}dx$$

$$（圖一：傅立葉級數的公式）$$

(圖二：將一個雜亂的函數表示成很多個有規律的sin/cos函數的疊加)

* 那沒規律的函數呢?
把一個無規律的函數圖形看成是有規律的，只是它的週期是無限大
然後把這個 "週期無限大" 的函數拆解成很多個$\sin$和$\cos$函數的疊加
$\sin$和$\cos$都是有規律的，所以電腦可以比較輕鬆的理解這段函數

* 頻域(頻譜圖)
圖一 b的左邊可以看到紅色的標記
此紅色的標記就是各個$\sin$ or $\cos$函數的振幅
一個$\sin$函數會有一種頻率，同時會對應到一個振幅
如果我們把各個$\sin$函數的頻率作為x軸，振幅作為y軸
這樣畫出來的圖，就是頻譜圖(頻域) (圖二 c)

$$(圖三：時域與頻域的對照圖) (source：改編自二十四橋KG)$$

* 光譜學上的用處
光譜即是一張時（空）域，我們對其進行傅立葉轉換後
便可得到許多不同頻率的$\sin$和$\cos$波形
如此一來，我們便可以輕鬆將一些特定頻率的訊號(例如：光害的訊號)去掉
然後再把剩下的訊號（天體發出的訊號）重新疊合在一起
即可得到乾淨的天體訊號（沒有大氣影響或光害的）

---

## 卷積(Convolution)
### 卷積的定義
:::info
卷積是將訊號疊加
在影像處理上類似於套用濾鏡
:::
有一間工廠生產麵包
遵循以下兩條規則

$$
\begin{aligned}
&f(x)=隨時間麵包的產量\\
&g(x)=隨時間麵包未腐壞率
\end{aligned}
$$

那如果想知道在第$t$個時間單位共有多少麵包呢?

我們可以先思考一下
在第0個時間單位生產的麵包數量為

$$f(0)g(t)$$

::: success
因為第0秒有$f(0)$個麵包被生產出來
又在$t$時間後
第0秒生產的麵包的未腐壞率為$g(t)$
:::
所以在第$x$秒生產的麵包數量為
$$f(x)g(t-x)$$
如果我們要求$t$秒內共有多少麵包
那就是
$$\sum^t_{x=0}f(x)g(t-x)$$
如果把時間當作連續
則變成
$$\int_0^t f(x)g(t-x)dx$$
這就是我們說的
**卷積是在做訊號疊加**
我們有兩個函數(訊號)$f(x)、g(x)$
就能藉由卷積將兩個函數疊加在一起

$$(f\ast g)(t)=\int_0^t f(\tau)g(t-\tau)d\tau$$

另一個例子
我們有原訊號及系統響應(也就是訊號隨時間會不斷衰減)


所以我們如果想求出實際訊號值
就必須將輸入訊號與系統響應疊加
那在時間$t$後
第0秒的訊號就會變成$f(0)g(t)$
又因為訊號是不停地被輸入
所以第$t$秒時的總訊號就是在$t$秒內的所有乘積和
即
$$\int_0^t f(x)g(t-x)dx$$
這樣應該就能理解為什麼算卷積時要將其中一個函數翻轉了


### 練習：用陣列算卷積

我們可以把陣列$[4,5,6]當作[f(1),f(2),f(3)]$
這樣就能用離散的方式練習卷積了

現在要算出
$[1,2,3]\ast[4,5,6]$

**算法如下：**

先把其中一個陣列反轉
然後依序移動做對應項的乘積和

目前答案: [**1*4**, , , , ]

目前答案: [4, **1\*5+2\*4**, , , ]

目前答案: [4, 13, **1\*6+2\*5+3\*4**, , ]

目前答案: [4, 13, 28, **2\*6+3\*5**, ]

目前答案: [4, 13, 28, 27, **3\*6**]

所以最後得出$[1,2,3]\ast[4,5,6]=[4, 13, 28, 27, 18]$

## 影像卷積（二維卷積）

### 圖的數值化


在圖像中會以8個bit儲存顏色資訊
也就是說每個像素可以調整出256種不同的顏色強度

彩色圖片則可以分為RGB三個channel
一樣以8bit分別儲存顏色資訊


因為圖被數值化了，所以我們可以對圖做二維卷積。

### 二維卷積

我們先前講過卷積其實是在計算輸入訊號受到系統響應的影響
所以圖片同理
我們要決定每個像素值對最終結果的影響大小
也可以說是要決定每個像素的權重

這些權重稱作卷積核(kernel，又稱為filter weight)
其大小和值是由使用者決定的
不同卷積核對圖片又會有不同影響
類似於圖片的濾鏡(下面有範例)

接著將kernel放在圖上計算對應格子的乘積
再像sliding window那樣將kernel滑過去，再計算
重複至整張圖都被跑過為止


## 影像卷積效果展示

>原圖
:::info
還要更快
:::

>```python=
>kernel = np.array([0])
>```
>:::info
>不意外的整張變黑了，因為每個像素都被*0
>:::

>```python=
>kernel = np.array([[1/81 for _ in range(9)] for __ in range(9)])
>```
>:::info
>使用9*9的kernel，每格的值都是1/81，所以整張圖被模糊了9倍
>:::

>```python=
>kernel = np.array(
>[[0,1,1,2,2,2,1,1,0],
>[1,2,4,5,5,5,4,2,1],
>[1,4,5,3,0,3,5,4,1],
>[2,5,3,-12,-24,-12,3,5,2],
>[2,5,0,-24,-40,-24,0,5,2],
>[2,5,3,-12,-24,-12,3,5,2],
>[1,4,5,3,0,3,4,4,1],
>[1,2,4,5,5,5,4,2,1],
>[0,1,1,2,2,2,1,1,0]])
>```
>:::info
>高斯拉普拉斯算子，可以偵測物體邊緣
>:::

>```python=
>kernel = np.array([[-1,-2,-1],[-2,13,-2],[-1,-2,-1]])
>```
>:::info
>銳利化
>:::

## 遊樂場
環境自己架
gui.py是主程式

https://github.com/ChiuDeYuan/conv_playground

## 去卷積（Deconvolution）
:::info
把圖片還原回來
:::
既然卷積是如此運算
可以把一張圖片變模糊
$$P\ast K=L$$

那能不能這樣做
就能恢復成清楚的照片了
$$P=L*K^{-1}$$
:::info
這裡卷積核右上方的-1表示這個卷積核是做反操作
:::

這個把圖片變清楚的過程就稱為**去卷積(Deconvolution)**
如果去卷積可以被計算出來
就能夠修復模糊的圖片了

這在天文照片的處理有很大的用處
因為天文照片時常受到大氣、光害等干擾
使用去卷積就能修復出較真實的圖像

:::warning
:warning: 還有其他因素可能會使照片失真
* 相機抖動
* 鏡頭缺陷(如球面像差)
* 光的繞射(Diffraction)
* 物體移動
* 景深限制
* 後期處理的干擾
:::

### 如何執行去卷積?
$$\begin{aligned}我們這樣&表示卷積\\x\;\ast\ &c=b\end{aligned}$$
$$\begin{aligned}空域中的卷積相當&於頻域中的乘法\\F(x)\;\cdot\ F(&c)=F(b)\end{aligned}$$
$$\begin{aligned}除法就相當&於去卷積\\F(x_{deconv})=&\;F(b)/F(c)\end{aligned}$$
$$\begin{aligned}再對結果作反傅立葉轉&換就能得到原圖了\\x_{deconv}=\;F^{-1}&(F(b)/F(c))\end{aligned}$$

<br>

### 點擴散函數(Point Spread Function, PSF)
::: info
用來描述光線的擴散、模糊程度
:::
:::success
PSF與物體的卷積就會是我們所看到的影像
:::
由於一些擾動(如鏡頭缺陷、光的繞射)，
導致影像與理想值有差距，
所以我們用PSF來描述影像的擴展斑點(影像光線的特徵)


你也可以說PSF就是$x_{deconv}=F^{-1}(F(b)/F(c))中的c$
只要知道一張模糊圖像的PSF
就能以去卷積還原圖像

---

## 非盲去卷積(Non-blind Deconvolution)

非盲卷積是指知道PSF的情況下的去卷積

### 維納濾波

由於圖像中可能會有噪聲
$$b=x\ast c+n$$
所以沒辦法直接用$x_{deconv}=F^{-1}(F(b)/F(c))$算出來
因為算出來只會是一片噪聲

所以Norbert Wiener在1942年提出了一個算法
公式如下
$$x_{deconv}=F^{-1}(\frac{\vert F(c)\vert^2}{\vert F(c)\vert^2+1/SNR(\omega)}{\cdot}\frac{F(b)}{F(c)})$$
:::info
$SNR(\omega)$是信噪比，當圖片失真越小，信噪比越大
:::

這個公式是通過最小化期望誤差而來
$$min_H E[\Vert X-H\tilde{B}\Vert^2]$$
:::info
$X=F(x)$
$H$是一個執行去卷積的函數
算出來也就是維納濾波的公式
最後公式是
$$x_{deconv}=F^{-1}(HB)=F^{-1}(\frac{\vert F(c)\vert^2}{\vert F(c)\vert^2+1/SNR(\omega)}{\cdot}\frac{F(b)}{F(c)})$$
:::

### Richardson-Lucy方法
這是用迭代方法算反卷積
$$\hat{f}_{k+1}(x)=\hat{f}_k(x,y)[h(-x,-y)\ast\frac{g(x,y)}{h(x,y)\ast \hat{f}_k(x,y)}]$$

---

## 盲去卷積(Blind Deconvolution)
>圖片來源&參考
>https://zhuanlan.zhihu.com/p/105500403

盲去卷積是在不知道PSF的情況下做去卷積
但其實就是先算出PSF後再用非盲去卷積

### 什麼是我們要的結果?
一般照片像素的顏色梯度分布是這樣的
因為邊界很清楚
所以會有許多梯度較大的區域
這種梯度分布叫做**Heavy-Tail**

:::info
梯度：你可以想像是顏色的變化量
:::
而比較模糊的照片的梯度分布可能會長這樣
因為色塊都糊在一起了(顏色邊界不清楚)
所以會有很多小梯度


### 計算PSF

現在我們可以構造估計PSF的式子
$$argmax_K\;p(K|P)$$
:::info
K是模糊核
P是模糊圖像
:::
再用 **變分貝葉斯(Variational Bayesian)** 來求取近似條件機率

### 整體流程
1. 預處理圖像：因為變分貝葉斯計算量大，所以要先選取一小塊圖像
2. 估計模糊核：使用變分貝葉斯
3. 重建圖像：知道PSF後，用非盲去卷積算法

---

## 天文照片處理

* 預處理(第一步)
1. 下載Siril (https://siril.org/)

1. 從USB中下載各個檔案資料夾，並把資料夾名稱中的-IC434去掉
(你的資料夾應該會變成lights、flats、darks、biases)
1. 打開Siril，點擊左上角藍色的房子圖案
1. 然後選取你下載的資料夾所在的位置(假設你下載在D槽，那就選D槽)
1. 點擊上方中間偏左的"Scripts"，然後點擊"OSC_Preprocessing"
(這一步驟叫做預處理"，通常需要一點時間)
1. 點擊左上角的"Open"，然後找到你剛剛預處理完的檔案，打開它
(通常檔案名會叫 "result+一串數字" ，副檔名是.fit檔)

* 線性處理(第二步)
1. 點擊下方中間偏左的"Linear"，然後改成"Auto Stretch"
1. 將鼠標移至畫面中左的照片上，按住左鍵，框選照片
1. 選取好你要的照片範圍後，在照片上按右鍵，然後點擊"crop"
(此步驟可以裁切照片，去除掉照片邊緣一些暗角之類的雜訊)
1. 對照片中看起來是背景星空的地方點選左鍵，生成紅色的點
(注意不要點到目標天體(中間的馬頭星雲和火焰星雲))
1. 按"compute background"，:warning: 然後按"apply"，再關閉"Background Extraction"的介面
1. 一樣是點擊上方的"Image Processing"，然後選取"Color Calibration"
1. 一樣按住右鍵，框選照片中看起來最暗的部分(最黑)，然後點擊上面的"use current section"(框選的部分不用很大)
1. 重複上面的步驟，但改成選取照片中最亮的部分，然後點擊下方的"use current section"，然後按"apply"

* 去卷積(第2.5步)
1. 一樣點擊上方的"Image Processing"，然後選取"Noise Reduction"
1. 將"Secondary denoising algorithm"改成"Anscombe VST"，並將下方的"Modulatuion"調到0.5~0.6之間，然後按"apply"
1. 接著，將照片的右下角框起來，然後像之前一樣，點右鍵，然後按"Crop"
1. 按右上角一個像"三"的符號，然後點"image information"，再按最下面的"Dynamic PSF"
1. 將左下角的"Profile Type"改為"Moffat"，然後按右下角那一排最左邊，三顆星星的符號
(此步驟會幫你偵測受大氣擾動後變模糊的星點，並生成它們的平均PSF)

1. 接著，點選上方的"Image Processing"，然後選取"Deconvolution"
1. 點擊上半部分的"PSF from stars"，看看PSF review是不是有出現一個點擴散函數的圖片
1. 將下半部分(non-blind deconvolution)的模式改為"Richardson-Lucy methods"，然後將"Gradient Descent"改成"multiplicative"
1. 更改"Iteration"，試試看星點去卷積後有什麼不一樣
(試完一次後，可以按左上角的Redo，再重新試一次)
1. 當你覺得去卷積後的星點看起來足夠像圓形後，你可以按Redo，把照片回復成大張的(未Crop成一小張之前的樣子)
1. 然後用你剛剛覺得滿意的去卷積參數對整張照片去卷積

* 非線性處理(第三步)
1. 將下方的"Auto Stretch"改回"Linear"
1. 點擊上方的"Image Processing"，然後選取"Asinh Transformation"
1. 介面中，調高"Stretch factor"會讓照片變亮，調高"Black point"會讓照片變暗
1. 你們要讓照片中馬頭星雲稍微出現，然後背景不能出現反白(過曝)的感覺，:warning:然後記得按"apply"

1. 一樣點擊上方的"Image Processing"，然後選取"Histogram Transformation"
1. 然後點擊介面中間偏右的齒輪，:warning:然後一樣記得按"apply"
1. 一樣點擊上方的"Image Processing"，然後選取"Remove Green Noise"，然後直接按apply

1. 一樣點擊上方的"Image Processing"，然後選取"Color Saturation"
1. 然後調整成自己喜歡你 yeah的顏色

### 就這樣！

## 參考
https://zhuanlan.zhihu.com/p/104340592
https://zhuanlan.zhihu.com/p/105500403