# mado: spline.c
[src/spline.c](https://github.com/sysprog21/mado/blob/main/src/spline.c)

[mado 發音](http://ilrdc.tw/grammar/sound/13/A2-10-8.mp3)
## Flow

## functions and types
### fixed-point types
* [Fixed-point operation with MCTS by vax-r](https://hackmd.io/@vax-r/linux2024-homework1#Fixed-point-operation-with-MCTS)
* Fixed-points standard ISO/IEC TR 18037
* [Fixed-point types in GCC](https://gcc.gnu.org/onlinedocs/gcc/Fixed-Point.html)

根據 ISO/IEC TR 18037:2008(E) 規格書，其中第 4.2 節 Detailed changes to ISO/IEC 9899:1999 給出基於 C99 擴增的 C 語言新規範。

* Clause 5.2.4.2.3 - Characteristics of fixed-point types
一定點數 $x$ 定義為：
$$
x= s\times n \times b^f
$$
其中，$s$ 為正負號(-1 或 1)；
$p$ 為精度；
$n$ 為一非負整數且小於 $b^p$；
$b$ 為一基數，在二進位系統中通常為 $2$；
$f$ 為一係數，用於表示小數部分的位數。

舉例來說，以 16 位元的有號整數來表示一定點數 $x$，其中 4 個位元用於表示分數部分，11 個位元用於表示整數部分，1 個位元用於表示正負號。$f=-4$、$p=11+4=15$、$0\leq n <2^{15}$。假如有一定點數為 $549$ 其二進位表示為 `0000 0010 0010 0101` 則此定點數 $x$ 實際數值為
$$
x = 549 \times 2^{-4} = 34.3125 = 2^5 + 2^1 + 2^{-2}+2^{-1}
$$

* Clause 6.2.6.3 - Fixed-point types
* 對有號定點數，可以將位元分成三個部分：
1. 正負號位元(sign bit)
2. 數值位元(value bits)
3. 填充位元(padding bits)

對於數值位元又可以分成：
1. 整數位元(integral bits)
2. 分數位元(fractional bits)

其中，假如數值位元為 $N$ 個且整數位元為 $L$ 個，則分數位元部分為 $N-L$ 個。對於用於表示純分數的定點數，該整數位元為 $L=0$ 個，則每一數值位元用於表示在 $2^{-1}$ 到 $2^{-N}$ 之間不同 $2$ 的冪，所以該定點數的實際數值範圍為 $-1$ 至 $1-2^N$；
對於用於表示帶分數的定點數，則每一數值位元用於表示在 $2^{L-1}$ 到 $2^{N-L}$ 之間不同 $2$ 的冪，所以該定點數的實際數值範圍為 $-2^L$ 至 $2^L-2^{N-L}$。
* 對無號定點數，可以將位元分成兩個部分：
1. 數值位元(value bits)
2. 填充位元(padding bits)

對於數值位元又可以分成：
1. 整數位元(integral bits)
2. 分數位元(fractional bits)

其中，假如數值位元為 $N$ 個且整數位元為 $L$ 個，則分數位元部分為 $N-L$ 個。對於用於表示純分數的定點數，該整數位元為 $L=0$ 個，則每一數值位元用於表示在 $2^{-1}$ 到 $2^{-N}$ 之間不同 $2$ 的冪，所以該定點數的實際數值範圍為 $0$ 至 $1-2^{-N}$；
對於用於表示帶分數的定點數，則每一數值位元用於表示在 $2^{L-1}$ 到 $2^{N-L}$ 之間不同 $2$ 的冪，所以該定點數的實際數值範圍為 $0$ 至 $2^L-2^{N-L}$。

typedef int16_t twin_sfixed_t; /* 12.4 format */
typedef int32_t twin_dfixed_t; /* 24.8 format */

在 `mado` 中採用的是

/* Geometrical objects */
typedef struct _twin_spoint {
    twin_sfixed_t x, y;
} twin_spoint_t;

### _twin_distance_to_line_squared

一條樣到一直線的距離。輸入有三個 `twin_spoint_t` 輸出為一 `twin_dfixed_t`。

回顧[點到直線距離](https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line)，直線方程的一般式為：
$$
Ax+By+C = 0
$$

假設有兩點 $(x_1,y_1)$ 和 $(x_2,y_2)$ 形成的直線寫成直線方程為
$$
(y-y_1)=\dfrac{y_2-y_1}{x_2-x_1}(x-x_1) \text{ and }
(y-y_2)=\dfrac{y_2-y_1}{x_2-x_1}(x-x_2)
$$
整理得到
$$
x(y_2-y_1) - y(x_2-x_1) - x_1(y_2-y_1) + y_1(x_2-x_1) = 0
$$

其中
$$
\begin{align}
A&=(y_2-y_1) \\
B&=(x_1-x_2) \\
C&=(y_1x_2-x_1y_2)
\end{align}
$$

則一直線方程外的點 $(x_0,y_0)$ 到該直線方程的距離則為：
$$
\dfrac{\lvert Ax_0+By_0+C \rvert}{\sqrt{A^2+B^2}}
$$

將其平方
$$
\dfrac{(Ax_0+By_0+C)^2}{A^2+B^2}
$$

### `_lerp`

為一線性插值函數，lerp 為 liner interploation 的縮寫

在 mado 中使用的曲線插值方式是貝茲曲線，兩點 $P_1$、$P_2$ 間的凸組合構成的插值方式
$$
B(t) = P_1 +(P_2-P_1)t = (1-t)P_1 + tP_2
$$
其中 $t\in[0,1]$。

把 $t$ 寫成 $\frac{1}{2}$ 的冪，則第二個等式可以表示為
$$
P_1 +(P_2-P_1)2^{-k}
$$
其中 $k$ 為一非負整數。

由於為 $(P_2-P_1)$ 除以 $2$ 的冪，所以改成用位元運算可以寫成
$$
P_1 +((P_2-P_1) >> k)
$$

但在此 $t$ 的選擇就被縮限在 $t=1,2^{-1},\dots,2^{-k},\dots$

而

static void _lerp(twin_spoint_t *a,
                  twin_spoint_t *b,
                  int shift,
                  twin_spoint_t *result)
{
    result->x = a->x + ((b->x - a->x) >> shift);
    result->y = a->y + ((b->y - a->y) >> shift);
}

### `_de_casteljau`
[Linear-time geometric algorithm for evaluating Bézier curves Filip Chudy, Paweł Woźny](https://arxiv.org/abs/1803.06843)

在 mado 中的貝茲曲線是三次貝茲曲線，這裡的三次是基於變數 $t$ 的三次。

根據一開始四個點 $a,b,c,d$ 依照順序兩兩做一次線性插值則可以得到
$$
\begin{align}
ab(t) &= (1-t)a + tb \\
bc(t) &= (1-t)b + tc \\
cd(t) &= (1-t)c + td
\end{align}
$$

其中 $ab(t),bc(t),cd(t)$ 為不同線性插植中的點，根據一次線性插植中的點做二次插值，
$$
\begin{align}
abbc(t) &= (1-t)ab(t) + tbc(t) \\
&= (1-t)((1-t)a + tb) + t((1-t)b + tc)\\
&= (1-t)^2a+2t(1-t)b+t^2c \\
bccd(t) &= (1-t)bc(t) + tcd(t) \\
&= (1-t)((1-t)b + tc) + t((1-t)c + td)\\
&= (1-t)^2b+2t(1-t)c+t^2d \\
\end{align}
$$

三次插值為
$$
\begin{align}
\text{final}(t)
&= (1-t)abbc(t) + tbccd(t) \\
&= (1-t)((1-t)^2a+2t(1-t)b+t^2c)\\
& + t((1-t)^2b+2t(1-t)c+t^2d) \\
&= (1-t)^3+ 3t(1-t)^2b + 3t^2(1-t)c + t^3d
\end{align}
$$

De Casteljau 的演算法則是藉由，每次迭代算出一線性插植點，再用於下一次迭代中。

* 藍色為初始點。

* 令插值函數 $t=2^-1$。紅色為第一次線性插值出來的點，黃色點為以紅色點做線性插值出來的點，黑色點為以黃色點做線性插值出來的點。

### `_twin_spline_decompose`

[include/twin_private.h](https://github.com/sysprog21/mado/blob/main/include/twin_private.h)

#define TWIN_SFIXED_ONE (0x10)
#define TWIN_SFIXED_TOLERANCE (TWIN_SFIXED_ONE >> 2)

`TWIN_SFIXED_ONE` in binary format is `0000 0000 0001 0000` where four bits for representing the fractional part.

And the default tolerance `TWIN_SFIXED_TOLERANCE` is `0000 0000 0000 0100` which is $\frac{1}{2^2}=0.25$ in deciaml in the 12.4 fixed-point representation.

`TWIN_SFIXED_TOLERANCE * TWIN_SFIXED_TOLERANCE`

基於容忍量 tolerance 選取 shift 的策略。

數學上的插值函數是建立在一連續定義域跟一對應域中，在電腦中是以二進位數值為定義域跟定義域，更具體的來說是一有限域。詳細的內容可以參考[解讀計算機編碼](https://hackmd.io/@sysprog/binary-representation)。

每執行一次三次插值法，就會多 $5$ 個切割線段

## flatness
https://blend2d.com/research/simplify_and_offset_bezier_curves.pdf

##

$$
d(p_0, L)^2 \leq \epsilon^2
$$

自己選的 $\epsilon$

$$
g(t)=d(p_0, L)^2 = \dfrac{((y_2-y_1)p_{x,0}+(x_1-x_2)p_{y,0}+(y_1x_2-x_1y_2))^2}{(y_2-y_1)^2+(x_1-x_2)^2} \leq \epsilon^2
$$

$$
\begin{align}
A&=(y_2-y_1) \\
B&=(x_1-x_2) \\
C&=(y_1x_2-x_1y_2)
\end{align}
$$

## 論文 GPU-friendly Stroke Expansion, Google 2024
2024 RAPH LEVIEN and ARMAN UGURAY

fig 3. $s$ 為分段圓形弧長，該分段圓形弧長的兩端點形成一線段，弦長 $d$ 為該分段圓形弧長至該線段的最遠距離
$$
d = (1-\cos{\theta})\frac{s}{2\theta}
$$
見
https://en.wikipedia.org/wiki/Circular_segment
https://en.wikipedia.org/wiki/Sagitta_(geometry)

對一曲線 $\hat{s}$ 曲線上任意點曲率 $\kappa(s)$ 其原文使用積分表示誤差
$$
d\approx \frac{1}{8} \left(\int^{\hat{s}}_0 \sqrt{|\kappa(s)|} ds\right)^2
$$