Mado

source code

請參考附件，展示類似 Cairo APIs [1] 風格的向量繪圖函式庫，後續的改進:

• 不依賴 libm，完全用 fixed point 建構
• 整合進 twin: https://github.com/Joouis/EVGUI_Viewer
• 以 perf 一類的工具找出效能瓶頸，並提出改進方案
• 在 Linux framebuffer 展現

[1] https://www.cairographics.org/manual/

TODO: Issue 2

問題說明 :

The usual method to convert a spline into a sequence of line segments is to split the spline in half using de Casteljau's algorithm recursively until the spline can be approximated by a straight line within a defined tolerance. The commit 6d85b42 has already changed the original recursion to an iterative method, avoiding unbounded recursion. Four levels of recursion will consume more than 64 coordinates' worth of space, which can easily overflow the stack of a tiny machine. However, as Keith mentioned in Slightly Better Iterative Spline Decomposition, it is possible to make the overall computational cost lower than a straightforward binary decomposition, and no recursion is required.
See also: Decomposing Splines Without Recursion

討論 :

jserv
The font-edit is not only useful for observing how font rendering works but also serves as a testbed for evaluating the enhancement of iterative spline decomposition, where you can operate with floating point arithmetic instead of fixed-point arithmetic.
jouae
Here are two books that might be helpful:
A Primer on Bézier Curves
This book gives some interactive content can visually aid in understanding the significance of parameters.

Computer Aided Geometric Design
Discussing curve interpolation with rigorous mathematical exposition.
In the second book 10.6 Error Bounds,
We can assure that the approximation error will be less than a specified tolerance by using line segments whose endpoints are evenly spaced in

古老的演算法

The usual method to convert a spline into a sequence of line segments is to split the spline in half using de Casteljau's algorithm recursively until the spline can be approximated by a straight line within a defined tolerance. The commit 6d85b42 has already changed the original recursion to an iterative method, avoiding unbounded recursion.

This is certainly straightforward, but suffers from an obvious flaw — there's unbounded recursion. With two splines in the stack frame, each containing eight coordinates, the stack will grow rapidly; 4 levels of recursion will consume more than 64 coordinates space. This can easily overflow the stack of a tiny machine.

iterative-splines 的實作

commit 6d85b42
先看這個結構體一個 spline 是由四個點所組成的：

• a 曲線端點之一，作為起始點。
• b 第一控制點，用於控制 a 點延伸的方向。
• c 第二控制點，用於控制 b 點延伸的方向。
• d 曲線端點之一，作為終止點。

typedef struct _twin_spline {
    twin_spoint_t a, b, c, d;
} twin_spline_t;

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

在 mado 中使用的曲線插值方式是貝茲曲線，兩點

$P_1$、

$P_2$ 間的凸組合構成的插值方式

$$B(t)=P_1+(P_2-P_1)t=(1-t)P_1+t{P}_2$$
其中

$t\in [0,1]$。

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);
}

把

$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$。

_de_casteljau 在一個 spline 中找到左邊的 spline 與右邊的 spline 。

原本以為 spline 的操作是這樣。

但後來發現當到 shift=3 的時候需要知道它已經平坦了，要確認某一個 spline 是否平坦時，需要去計算，而下面的的算平坦的方式是最簡單的道理。算點與直線的最短距離。
然而同時也發現這樣的操作，當在 shift = 3 的時候確認已經足夠平坦後，到 shift = 4，便可以看到這個趨勢，所計算的點的所落在的位置已經越來越接近點 A 與點 B 了，這相當於失去控制點 B 與 C 原先對於最終產生的曲線的權重。會讓這個曲線越來越平。就如同下圖。

twin_dfixed_t _twin_distance_to_line_squared(twin_spoint_t *p,
                                             twin_spoint_t *p1,
                                             twin_spoint_t *p2)
{
    /*
     * Convert to normal form (AX + BY + C = 0)
     *
     * (X - x1) * (y2 - y1) = (Y - y1) * (x2 - x1)
     *
     * X * (y2 - y1) - Y * (x2 - x1) - x1 * (y2 - y1) + y1 * (x2 - x1) = 0
     * A = (y2 - y1)
     * B = (x1 - x2)
     * C = (y1x2 - x1y2)
     *
     * distance² = (AX + BC + C)² / (A² + B²)
     */
    twin_dfixed_t A = p2->y - p1->y;
    twin_dfixed_t B = p1->x - p2->x;
    twin_dfixed_t C =
        ((twin_dfixed_t) p1->y * p2->x - (twin_dfixed_t) p1->x * p2->y);
    twin_dfixed_t den, num;

    num = A * p->x + B * p->y + C;
    if (num < 0)
        num = -num;
    den = A * A + B * B;
    if (den == 0 || num >= 0x8000)
        return _twin_distance_to_point_squared(p, p1);
    return (num * num) / den;
}

這段程式碼是用來計算一個點 p 到直線(由點 p1 與點 p2 連線而成)的距離的平方。直線的點斜式方程：

兩點

$p1(x1,y1)$ 和

$p2(x2,y2)$ :

$(y-y1)=\frac{y2-y1}{x2-x1}\cdot (x-x1)$

$A\cdot X+B\cdot Y+C=0$

$A=y2-y1$,

$B=x1-x2$,

$C=y1\cdot x2-x1\cdot y2$

Calculate numerator and denominator:
num 代表 Numerator 分子 : num

$=A\cdot x+B\cdot y+C$
den 代表 Denominator 分母 : den

$=A^2+B^2$

點到直線的距離公式：
給定一個點

$p(x,y)$，它到直線

$A\cdot X+B\cdot Y+C=0$ 的距離的平方可以表示為：

${\text{distance}}^2=\frac{(A\cdot x+B\cdot y+C{)}^2}{A^2+B^2}$
這裡的 num 就是

$A\cdot x+B\cdot y+C$，而 den 是

$A^2+B^2$。

Check conditions:
如果分母 ( den ) 等於零，或者 ( num ) 大於等於 0x8000（十進制的32768），則返回到另外一個函數 _twin_distance_to_point_squared(p, p1)。

Compute squared distance:
否則，計算距離的平方：

${\text{distance}}^2=\frac{nu{m}^2}{den}$

最後返回計算出來的距離的平方。

/*
 * Return an upper bound on the error (squared) that could result from
 * approximating a spline with a line segment connecting the two endpoints.
 */
static twin_dfixed_t _twin_spline_error_squared(twin_spline_t *spline)
{
    twin_dfixed_t berr, cerr;

    berr = _twin_distance_to_line_squared(&spline->b, &spline->a, &spline->d);
    cerr = _twin_distance_to_line_squared(&spline->c, &spline->a, &spline->d);

    if (berr > cerr)
        return berr;
    return cerr;
}

看點 b 點 c 有沒有夠接近點 a 點 d 連成一條線。

/*
 * Check if a spline is flat enough by comparing the error against the
 * tolerance.
 */
static bool is_flat(twin_spline_t *spline, twin_dfixed_t tolerance_squared)
{
    return _twin_spline_error_squared(spline) <= tolerance_squared;
}

這段程式碼的意義，因為他是一直更新最左邊的曲線，讓最左邊的曲線越來越平，而如果最左邊的曲線足夠平坦就更新其右側。然而更新最左邊的點，透過

$\frac{1}{2}\to \frac{1}{4}$ 來讓點 AB 與點 ABBC 非常接近點 A 點與點 final 所連成的線，如下圖。

如果夠接近就將右側的視為新的 spline，重新開始用新的 shift 計算，shift = 1 就會變成

$t=\frac{1}{2}$ 。如果不夠平坦再繼續 shift = 2

$t=\frac{1}{4}$，直到左側的曲線為平坦後，就會將右側的曲線視為新的曲線，並從 shift = 1 開始。

在看左側點 AB 與點 ABBC 是否非常接近點 A 點與點 final 所連成的線。因此這樣就不會一直失去原先控制點 B 與 C 的權重。但這件事情還在思考是否正確，所以要證明這個演算法 de Casteljau Algorithm 真的不會失去其權重。

de Casteljau Algorithm

Suppose

$b(t)$ is the degree n Bezier curve for control points

$b_0,\dots ,b_n$. To find

$b(t)$ proceed as follows: Define:

$b_{0,i}=b_i$ For

$i=0,\dots ,n$

$b_{1,i}=(1-t)b_{0,i}+t\times b_{0,i+1}$ For

$i=0,\dots ,n-1$

$b_{2,i}=(1-t)b_{1,i}+t\times b_{1,i+1}$ For

$i=0,\dots ,n-2$

$b_{j,i}=(1-t)b_{j-1,i}+t\times b_{j-1,i+1}$ For

$i=0,\dots ,n-j$

$b_{n,i}=(1-t)b_{n-1,i}+t\times b_{n-1,i+1}$ For

$i=0$

展開這些行

$b_{0,0}=b_0{b}_{0,1}=b_1{b}_{0,2}=b_2\dots b_{0,n-2}=b_{n-2}{b}_{0,n-1}=b_{n-1}{b}_{0,n}=b_n$

$b_{1,0}{b}_{1,1}{b}_{1,2}\dots b_{1,n-2}{b}_{1,n-1}$

$b_{2,0}{b}_{2,1}{b}_{2,2}\dots b_{2,n-2}$

A :

$b_{0,0}=[1,1]$ B :

$b_{0,1}=[2,8]$ C :

$b_{0,2}=[6,0]$ D :

$b_{0,3}=[8,7]$
AB :

$b_{1,0}=[1\times 0.9+2\times 0.1,1\times 0.9+8\times 0.1]=[1.1,1.7]$
BC :

$b_{1,1}=[2\times 0.9+6\times 0.1,8\times 0.9+0\times 0.1]=[2.4,7.2]$
CD :

$b_{1,2}=[8\times 0.9+6\times 0.1,0\times 0.9+7\times 0.1]=[6.2,0.7]$

ABBC :

$b_{2,0}=[1.1\times 0.9+2.4\times 0.1,1.7\times 0.9+7.2\times 0.1]=[1.23,2.25]$
BCCD :

$b_{2,1}=[2.4\times 0.9+6.2\times 0.1,7.2\times 0.9+0.7\times 0.1]=[2.78,6.55]$
final :

$b_{3,0}=[1.23\times 0.9+2.78\times 0.1,2.25\times 0.9+6.55\times 0.1]=[1.385,2.68]$

而前面所說的取左邊相當於現在只要計算

$[0,0.1]$ 的區間 而

$[0.1,1]$ 的區間則忽略計算。而下個定理則來說明這個方法要如何達到。

Theorem 3.0.2.
Suppose

$b(t)$ is the parametrization of the degree n (在此例為 3) Bezier curve with control points

$b_0,\dots ,b_{n=3}$. Suppose

$t_0$ ∈

$[0,1]$ (在此例中為 0.1) is chosen. The curves parametrized by

$b_t$ on

$[0,t_0]$ and

$b_t$ on

$[t_0,1]$ are also degree n (在此例為 3) Bezier curves and moreover we may obtain control points for these two curves according to the following:

對於左邊的使用為 A :

$b_{0,0}$ , AB :

$b_{1,0}$, ABBC :

$b_{2,0}$, final :

$b_{3,0}$
對於右邊的使用為 final :

$b_{0,0}$ , BCCD :

$b_{3,1}$, CD :

$b_{2,2}$, D :

$b_{1,3}$

	0	1	2	3
0	A	B	C	D
1	AB	BC	CD
2	ABBC	BCCD
3	final

AB =

$A\times (1-t)+B\times t$
BC =

$B\times (1-t)+C\times t$
CD =

$C\times (1-t)+D\times t$
ABBC = AB

$\times (1-t)$ + BC

$\times t$
展開
ABBC = [

$A\times (1-t)+B\times t$ ]

$\times (1-t)$ + [

$B\times (1-t)+C\times t$ ]

$\times t$ 化簡
ABBC =

$A\times (1-t{)}^2+2\times B\times t\times (1-t)+C\times t^2$

BCCD = BC

$\times (1-t)$ + CD

$\times t$
展開
BCCD = [

$B\times (1-t)+C\times t$ ]

$\times (1-t)$ + [

$C\times (1-t)+D\times t$ ]

$\times t$ 化簡
BCCD =

$B\times (1-t{)}^2+2\times C\times t\times (1-t)+D\times t^2$

看到這裡就會發現這是一個巴斯卡三角形
final = [

$A\times (1-t{)}^2+2\times B\times t\times (1-t)+C\times t^2$ ]

$\times (1-t)$ + [

$B\times (1-t{)}^2+2\times C\times t\times (1-t)+D\times t^2$ ]

$\times t$
final =

$A\times (1-t{)}^3+3\times B\times t\times (1-t{)}^2+3\times C\times t^2\times (1-t)+D\times t^3$

因此這個演算法是相當於第一次就是畫出最一開始的 A 點。接下來剖成左邊的 spline 與右邊的 spline 。當左邊的 spline 也就是 AB ABBC 夠接近 A 與 final 所做的線時，代表此次計算的這個 final 足夠讓線條夠平滑。

當這個 final 足夠讓線條夠平滑時，則畫這個 A 點。並將點 final 、點 BCCD 、點 CD 、點 D 視為下次要做的 spline 的點 A 、點 B 、點 C 、點 D。因為每次都是畫 A 點，所以這次的結果 final 這個點就會被畫進去，因此可以理解 _twin_spline_decompose() 這個函式在第 26 行要畫 D 點的原因。

static void _twin_spline_decompose(twin_path_t *path,
                                   twin_spline_t *spline,
                                   twin_dfixed_t tolerance_squared)
{
    /* Draw starting point */
    _twin_path_sdraw(path, spline->a.x, spline->a.y);

    while (!is_flat(spline, tolerance_squared)) {
        int shift = 1;
        twin_spline_t left, right;

        /* FIXME: Find the optimal shift value to decompose the spline */
        do {
            _de_casteljau(spline, shift, &left, &right);
            shift++;
        } while (!is_flat(&left, tolerance_squared));

        /* Draw the left segment */
        _twin_path_sdraw(path, left.a.x, left.a.y);

        /* Update spline to the right segment */
        memcpy(spline, &right, sizeof(twin_spline_t));
    }

    /* Draw the ending point */
    _twin_path_sdraw(path, spline->d.x, spline->d.y);
}

然而我還是想要理解為何點 final 、點 BCCD 、點 CD 、點 D 組成的 spline 是否會合原先的 spline 重疊？ 推導：現在的點 A 、點 B 、點 C 、點 D 分別為上一輪點 final 、點 BCCD 、點 CD 、點 D /src/spline.c#L47

    s1->a = spline->a;
    s1->b = ab;
    s1->c = abbc;
    s1->d = final;

    s2->a = final;
    s2->b = bccd;
    s2->c = cd;
    s2->d = spline->d;

證明 de Casteljau Algorithm

A :

$A\times (1-t{)}^3+3\times B\times t\times (1-t{)}^2+3\times C\times t^2\times (1-t)+D\times t^3$
B :

$B\times (1-t{)}^2+2\times C\times t\times (1-t)+D\times t^2$
C :

$C\times (1-t)+D\times t$
D :

$D$
剛剛證明出 final 會是當前那一輪點 A 、點 B 、點 C 、點 D 由下列公式計算出來
final =

$A\times (1-T{)}^3+3\times B\times T\times (1-T{)}^2+3\times C\times T^2\times (1-T)+D\times T^3$
將這一輪的點 A 、點 B 、點 C 、點 D 代入：
final =
[

$A\times (1-t{)}^3+$

$3\times B\times t\times (1-t{)}^2$

$+$

$3\times C\times t^2\times (1-t)$

$+$

$D\times t^3$ ]

$\times (1-T{)}^3+$
3

$\times$[

$B\times (1-t{)}^2$

$+$

$2\times C\times t\times (1-t)$

$+D\times t^2$]

$\times T\times (1-T{)}^2+$
3

$\times$[

$C\times (1-t)$

$+D\times t$]

$\times T^2\times (1-T)+$

$D$

$\times T^3$
展開
final =

$A\times (1-t{)}^3(1-T{)}^3$

$+$

$3\times B\times t\times (1-t{)}^2\times (1-T{)}^3$

$+$

$3\times B\times (1-t{)}^2\times T\times (1-T{)}^2$

$+$

$3\times C\times t^2\times (1-t)\times (1-T{)}^3$

$+$

$6\times C\times t\times (1-t)\times T\times (1-T{)}^2$

$+$

$3\times C\times (1-t)\times T^2\times (1-T)$

$+$

$[D\times t^3\times (1-T{)}^3+3\times D\times t^2\times T\times (1-T{)}^2+3\times D\times t\times T^2\times (1-T)+D]$

$\times T^3$

經過觀察可以發現現在的

$t$ 應該改為

$(1-t)(1-T)$ 另為

$X$ ，

$(1-X)=(T+t-Tt)$
final =

$A\times (X{)}^3$

$+$

$3\times B\times t\times (X{)}^2\times (1-T)+3\times B\times (X{)}^2\times T$

$+$

$3\times C\times t^2\times (X)\times (1-T{)}^2+6\times C\times (X)\times t\times T\times (1-T)+3\times C\times (X)\times T^2$

$+[D\times t^3\times (1-T{)}^3+3\times D\times t^2\times T\times (1-T{)}^2+3\times D\times t\times T^2\times (1-T)+D]$

$\times T^3$

化簡第二列（藍色）的

$3\times B\times t\times (X{)}^2\times (1-T)+3\times B\times (X{)}^2\times T$

$=3\times B\times X^2\times (t-tT+T)=3\times B\times X^2\times (1-X)$

$3\times C\times t^2\times (X)\times (1-T{)}^2+6\times C\times (X)\times t\times T\times (1-T)+3\times C\times (X)\times T^2$

$=3\times C\times X\times (t^2\times (1-T{)}^2+2\times t\times T\times (1-T)+T^2)$

$=3\times C\times X\times (t^2\times (1-2T+T^2)+2tT-2t{T}^2+T^2)$

$=3\times C\times X\times (T^2+2tT+t^2-2t{T}^2-2t^{2}T+t^2{T}^2)$

$=3\times C\times X\times ((T+t{)}^2-2tT(T+t)+t^2{T}^2)$

$=3\times C\times X\times ((T+t)-tT{)}^2=3\times C\times X\times (1-X{)}^2$ 得證

$D\times (t^3\times (1-T{)}^3+3\times t^2\times T\times (1-T{)}^2+3\times t\times T^2\times (1-T)+1)$

$\times T^3$

$=D\times ((t^3-3t^{3}T+3t^3{T}^2-t^3{T}^3)+3t^2\times (T-2T^2+T^3)+3t\times (T^2-T^3)+1)$

$\times T^3$

$=D\times (t^3+3t^{2}T+3t{T}^2-3t{T}^3-6t^2{T}^2-3t^{3}T+3t^2{T}^3+3t^3{T}^2-t^3{T}^3+1)$

$\times T^3$

$=D\times (t^3+3t^{2}T+3t{T}^2-3tT(T^2)-3tT(2tT)-3tT{t}^2+3t^2{T}^2(T)+3t^2(t)T^2-t^3{T}^3+1)$

$\times T^3$

$=D\times (t^3+3t^{2}T+3t{T}^2-3tT(T^2+2tT+t^2)+3t^2{T}^2(T+t)-t^3{T}^3+1)$

$\times T^3$

$=D\times (t^3+3t^{2}T+3t{T}^2-3tT(T+t{)}^2+3t^2{T}^2(T+t)-t^3{T}^3+1)$

$\times T^3$

$=D\times (t^3+3t^{2}T+3t{T}^2+$

$T^3$

$-3tT(T+t{)}^2+3t^2{T}^2(T+t)-t^3{T}^3+1-$

$T^3$

$)$

$\times T^3$

$=D\times ((T+t{)}^3-3tT(T+t{)}^2+3t^2{T}^2(T+t)-t^3{T}^3+1-$

$T^3$

$)$

$\times T^3$

$=D\times ((T+t-tT{)}^3+1-$

$T^3$

$)$

$\times T^3$

$=D\times ((1-X{)}^3+1-$

$T^3$

$)$

$\times T^3$

然而一開始假設為
final =

$A\times (1-T{)}^3+3\times B\times T\times (1-T{)}^2+3\times C\times T^2\times (1-T)+D\times T^3$
因此為了符合上式

$X=(1-T)\to T=1-X\to T^3=(1-X{)}^3$
最後一項為

$D\times (1-X{)}^3$ 得證。

測試一

我原先想直接用權重的方式去計算：

因此實作了以下的版本：

static void _bezier(twin_spline_t *spline, twin_dfixed_t t, twin_spoint_t *p)
{
    twin_dfixed_t t2 = t * t;
    twin_dfixed_t t3 = t2 * t;
    twin_dfixed_t mt = 1 - t;
    twin_dfixed_t mt2 = mt * mt;
    twin_dfixed_t mt3 = mt2 * mt;

    p->x = spline->a.x * mt3 + spline->b.x * 3 * mt2 * t + spline->c.x * 3 * mt * t2 + spline->d.x * t3;
    p->y = spline->a.y * mt3 + spline->b.y * 3 * mt2 * t + spline->c.y * 3 * mt * t2 + spline->d.y * t3;
}

然後 _new_twin_spline_decompose() 會長這樣。

static void _new_twin_spline_decompose( twin_path_t *path,
                                        twin_spline_t *spline,
                                        twin_dfixed_t tolerance_squared)
{
    /* Draw starting point */
    twin_dfixed_t start = TWIN_SFIXED_HALF;
    twin_spoint_t p;
    while (start != TWIN_SFIXED_ONE){
        _bezier(spline, start, &p);
        if(_twin_distance_to_line_squared(&p, &spline->a, &spline->d) <= tolerance_squared){
            _twin_path_sdraw(path, p.x, p.y);
            start << 1;
        }
        else
            start >> 1;
    }
}

然而這樣會有精度的問題，因為定點數如果要處理 t =

$\frac{1}{16}$ 的時候就遇到困難了。這樣 t 就會被乘三次， 要處理 overflow 的成本太大。而 Bezier Curves 的計算也很複雜

Bezier Curves
Even though the formula for a Bezier curve is explicit and closed it is fairly computationally intensive. For any given t-value we must calculate a sum and each term in the sum involves an (

${}_m^n$) function which is itself fairly intensive.

因為這樣的嘗試後，才有意識到 de Casteljau Algorithm 的用意，此外原本這段程式碼只會更新左邊，我在想為何不兩邊都同時跟新，但經過這個推算才知道，它其實是在找下一個 t 應該為哪？
首先會先看 t =

$\frac{1}{2}$ 時左半邊是否平滑，也就相當於 final 到原先的 A 是否很渺小。然而如果左邊不平滑時，可以嘗試 t =

$\frac{1}{4}$ 時左半邊是否平滑。最後假如找到 t =

$\frac{1}{16}$ 這個點，就再算 t =

$\frac{1}{16}$ 到 t = 1 。

版本二： more-iterative-splines

如果要用第二個版本就要放棄 bitwise 的操作。同時原先覺的它只是讓這個

$t$ 不要只是

$\frac{1}{2}\to \frac{1}{4}$。而是可以變化到大於

$\frac{1}{2}$ 的地方。但後來開始看 Bisection method

iterative splines

判斷有沒有 flat

Flattening quadratic Béziers

這篇 blog 主要為 quadratic Bézier curves，這篇blog提及

It is so good that it’s probably best to flatten other curve types (including cubic Béziers) by converting them to quadratics first, then applying the technique of this blog post.

首先它說明

The core insight of this blog post is a closed-form analytic expression for the number of line segments needed in the continuous limit.

第一個概念

$$\delta y\approx \frac{1}{8}\kappa \delta s^2$$
其中

$\kappa$ 代表曲率(curvature)，用於描述一曲線的彎曲程度。

接下來會用

$y=x^2$ 去計算。

// Determine the x values and scaling to map to y=x^2
map_to_basic() {
    const ddx = 2 * this.x1 - this.x0 - this.x2;
    const ddy = 2 * this.y1 - this.y0 - this.y2;
    const u0 = (this.x1 - this.x0) * ddx + (this.y1 - this.y0) * ddy;
    const u2 = (this.x2 - this.x1) * ddx + (this.y2 - this.y1) * ddy;
    const cross = (this.x2 - this.x0) * ddy - (this.y2 - this.y0) * ddx;
    const x0 = u0 / cross;
    const x2 = u2 / cross;
    // There's probably a more elegant formulation of this...
    const scale = Math.abs(cross) / (Math.hypot(ddx, ddy) * Math.abs(x2 - x0));
    return {x0: x0, x2: x2, scale: scale, cross: cross};
}

但這個演算法之所以可以很快，是因為首先它可以很容易的用

$t$ 算出曲線，因為他是 quadradic beziers 上面計算，但是在 Cubic Bezier 上面計算就會有先前要計算

$\text{final}=A\times t^3+3\times B\times t^2(1-t)+3\times C\times t\times (1-t{)}^2+D(1-t{)}^3$ 計算大量的問題。因為也還不知道要如何將 Cubic Bezier 轉為 Quadradic Bezier ，所以先去讀下一篇。

    my_subdiv(tol) {
        // Map quadratic bezier segment to y = x^2 parabola.
        const params = this.map_to_basic();
        // Compute approximate integral for x at both endpoints.
        const a0 = approx_myint(params.x0);
        const a2 = approx_myint(params.x2);
        const count = 0.5 * Math.abs(a2 - a0) * Math.sqrt(params.scale / tol);
        const n = Math.ceil(count);
        const x0 = approx_inv_myint(a0);
        const x2 = approx_inv_myint(a2);
        let result = [0];
        // Subdivide evenly and compute approximate inverse integral.
        for (let i = 1; i < n; i++) {
            const x = approx_inv_myint(a0 + ((a2 - a0) * i) / n);
            // Map x parameter back to t parameter for the original segment.
            const t = (x - x0) / (x2 - x0);
            result.push(t);
        }
        result.push(1);
        return result;
    }

論文閱讀：Precise Flattening of Cubic Bezier Segments

Abstract

這篇論文就是提出一個想法要將 polyline approximation 的 Cubic Bezier 曲線平坦化。透過產生

$\frac{2}{3}$ 的線段。提昇 37% 的效率。

首先先用前面 Cubic Bezier 的參數公式

$Q(t)=(x(t),y(t))$ 為

$x(t)=(1-t{)}^3{x}_1+3t(1-t{)}^2{x}_2+3t^2(1-t)x_3+t^3{x}_4$

$y(t)=(1-t{)}^3{y}_1+3t(1-t{)}^2{y}_2+3t^2(1-t)y_3+t^3{y}_4$
Cubic Bezier 由

$P_1$

$P_2$

$P_3$

$P_3$ 所組成，

$P_1$ 為

$t=0$ 的點也就是貝茲曲線的起點。

接下來定義兩個單位向量

$\hat{r}$ (速度)與

$\hat{s}$ (與

$\hat{r}$ 是正交(垂直的向量) 的

$\hat{r}$ 與

$\hat{s}$ 內積為 0 )。

$$\hat{r}=\frac{P2-P1}{|P2-P1|}=(\frac{x_2-x_1}{\sqrt{(x_2-x_1{)}^2+(y_2-y_1{)}^2}},\frac{y_2-y_1}{\sqrt{(x_2-x_1{)}^2+(y_2-y_1{)}^2}})$$

$$\hat{s}=\frac{P2-P1}{|P2-P1|}=(\frac{y_2-y_1}{\sqrt{(x_2-x_1{)}^2+(y_2-y_1{)}^2}},\frac{-(x_2-x_1)}{\sqrt{(x_2-x_1{)}^2+(y_2-y_1{)}^2}})$$

點

$P(x,y)$ 的

$r=(P-P_1)\cdot \hat{r}=\frac{(x-x_1)(x_2-x_1)+(y-y_1)(y_2-y_1)}{\sqrt{(x_2-x_1{)}^2+(y_2-y_1{)}^2}}$

$s=(P-P_1)\cdot \hat{s}=\frac{(x-x_1)(y_2-y_1)-(y-y_1)(x_2-x_1)}{\sqrt{(x_2-x_1{)}^2+(y_2-y_1{)}^2}}$

接下來定義控制點，用先前的表格

	0	1	2	3
0	A (P1)	B	C	D (P4)
1	AB (P1')	BC	CD (P3')
2	ABBC (P1'')	BCCD (P2'')
3	final (P1''')

接著在這個座標中，控制點分別為 P₁ (r₁,s₁) P₂ (r₂,s₂) P₃ (r₃,s₃) P₄ (r₄,s₄)，接著考慮在這個座標上的貝茲曲線
s(t) =

$(1-t{)}^3{s}_1+3t(1-t{)}^2{s}_2+3t^2(1-t)s_3+t^3{s}_4$ 。接著 r₁ = s₁ = s₂ = 0 。我原本很好奇為何這個座標系統上 r₁ = s₁ = 0 ，但後來發現這其實就只是將原先的x y座標轉換過來而已。而根據上方點

$P(x,y)$ 的公式，

$r=(P-P_1)\cdot \hat{r}$ 將

$P_1$ 代入

$r_1$

$s_1$ 的確為 0 。

$P_2$ 帶入

$s_2$ 也為 0 。
所以 s(t) 變為
s(t) =

$3t^2(1-t)s_3+t^3{s}_4=3s_3{t}^2+(s_4-3s_3)t^3$ 。
在曲線一開始時，

$t\to 0$， s(t)

$=3s_3{t}^2$

For small values of t, and assuming that

$s_3$ is not close to zero (which would occur if the beginning of the curve were near an inflection point) the first term dominates.

這個曲線的速度為

$\frac{s}{r}$ ，進一步假設沿著這個曲線的加速度

$\frac{ds}{dr}$ 相當小，這個假設建立在是線段如果夠短，就像下圖的點 P₁ 與點 P 的連線接近線段時，沿著曲線的切線的變化會很小。但如果在尖端的話這件事情則不會成立。

If we further assume the acceleration along the curve is reasonably small, i.e., the value of the r-coordinate varies reasonably linearly with t (this assumption is met in all places except near a cusp point, which will be handled as a separate case,)

接下來相當於找到一個點其與該線段差距最大，而這個點的特質就是斜率會與這個線段相同。但這一篇論文一堆特殊的情況要處理。
但我想用「最適合的點 t 這個的特質就是斜率會與這個線段」相同的想法。來找到 t 點，這樣就可以減少曲線的 line segement 。

Bisection 的方式

/*
 * Decompose 's' into straight lines which are
 * within SNEK_DRAW_TOLERANCE of the spline
 */
static void
_spline_decompose(void (*draw)(float x, float y), spline_t s)
{
    /* Start at the beginning of the spline. */
    (*draw)(s[0][0], s[0][1]);

    /* Split the spline until it is flat enough */
    while (_flatness(s) > SCALE_FLAT(SNEK_DRAW_TOLERANCE)) {
        spline_t    s1, s2;
        float       hi = 1.0f;
        float       lo = 0.0f;

        /* Search for an initial section of the spline which
         * is flat, but not too flat
         */
        for (;;) {

            /* Average the lo and hi values for our
             * next estimate
             */
            float t = (hi + lo) / 2.0f;

            /* Split the spline at the target location
             */
            _de_casteljau(s, s1, s2, t);

            /* Compute the flatness and see if s1 is flat
             * enough
             */
            float flat = _flatness(s1);

            if (flat <= SCALE_FLAT(SNEK_DRAW_TOLERANCE)) {

                /* Stop looking when s1 is close
                 * enough to the target tolerance
                 */
                if (flat >= SCALE_FLAT(SNEK_DRAW_TOLERANCE - SNEK_FLAT_TOLERANCE))
                    break;

                /* Flat: t is the new lower interval bound */
                lo = t;
            } else {

                /* Not flat: t is the new upper interval bound */
                hi =  t;
            }
        }

        /* Draw to the end of s1 */
        (*draw)(s1[3][0], s1[3][1]);

        /* Replace s with s2 */
        memcpy(&s[0], &s2[0], sizeof (spline_t));
    }

    /* S is now flat enough, so draw to the end */
    (*draw)(s[3][0], s[3][1]);
}

我的想法

首先一開始我先用第一個

$t=0.5$ 找到第一個點。

再來嘗試用

$t=0.4$ 找到找到第一個點。

可以看到如果只將曲線做一次 line segment 用 t = 0.4 的時候會更接近原先的曲線。再進一步思考這樣更接近有什用處，這樣左半邊的曲線就會因此減少逼近的次數。右半邊也是這就是。

而運用剛剛的想法「最適合的點 t 這個的特質就是切線斜率會與這個線段」，但是我要算一個點的切線斜率就代表我要找到跟原先點 A 與點 D 距離最要遠的線段。但是算多遠也就是用 _twin_distance_to_line_squared() 這個函式，

Original

設定

TWIN_SFIXED_ONE 10
TWIN_SFIXED_TOLERANCE 4
TWIN_SFIXED_TOLERANCE^2 10

commit 1f096f3

實驗結果

Experiment 1:
Original - Average number of _de_casteljau calls per point: 1.99
Original - Average points per character: 18.89

Bisection

首先 besection 有 hi 與 lo 兩個值，所以將這兩個值用 twin_sfixed_t 去存。
並且將 lerp 改為以下，可以計算不同的 t 值，將 t 值不只限制在

$\frac{1}{2}$ 與

$\frac{1}{4}$ 與

$\frac{1}{8}$ 中。

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

shift 改為 twin_sfixed_t

static void dcj(twin_spline_t *spline,
                twin_sfixed_t shift, 
                twin_spline_t *s1, 
                twin_spline_t *s2)

bisection 不是一個好方法，找 tolerance 與 flat tolerance 的 threshold 困難，同時很難超越原先的點點數量，跟 iteration 數量，甚至時常會無法收斂，導致計算 tolerance 會不斷徘徊，不斷地收斂。

#define TWIN_SFIXED_TOLERANCE (TWIN_SFIXED_ONE >> 2)
#define TWIN_FLEX_TOLERANCE (TWIN_SFIXED_TOLERANCE - TWIN_SFIXED_TOLERANCE >> 2)

tolerance_squared = TWIN_SFIXED_TOLERANCE*TWIN_SFIXED_TOLERANCE

static void decomp(pts_t *pts, twin_spline_t *spline, twin_dfixed_t tolerance_squared, 
                    int *iter, int *points)
{
    while (!is_flat(spline, tolerance_squared)){

        twin_sfixed_t hi = TWIN_SFIXED_ONE;
        twin_sfixed_t t, lo = 0;
        twin_spline_t left, right;

        while(true){
            t = (hi + lo) >> 1;
            dcj(spline, t, &left, &right);
            /* Calcuate the iteration times */
            (*iter) = (*iter)+1;
            twin_dfixed_t error = _twin_spline_error_squared(&left);
            
            if(error <=  tolerance_square){
                if(error >= TWIN_FLEX_TOLERANCE*TWIN_FLEX_TOLERANCE)
                    break;
                lo = t;
            }
            else
                hi = t;
            
        }
        /* Draw the left segment */
        add_pt(pts, &left.a);
        (*points) =  (*points) + 1;

        /* Update spline to the right segment */
        memcpy(spline, &right, sizeof(twin_spline_t));
    }
}

首先 bisection 在一個單調遞增與單調遞減的函式逼近會比較容易實作。因為在一起取線上很難知道，現在點 t 點相當於落在曲線的哪裡。就像爬一座山，知道自己落在高處 500 公尺沒有用，因為你不知道你是在左側的山上 500 公尺，還是在右側的山上 500 公尺。如果在左側山上就改變 lo 的值會達到山峰，如果在右側山上就改變 hi 的值才會到達山峰。

比較distance 有無被update

彈性 update

接著想法「最適合的點 t 這個的特質就是切線斜率會與這個線段，而這個點與點 A 與點 D 的距離會是最遠的」。像是現在 t 的切線斜率就與 A 點與 D 點相同。因為切線斜率如果相對於 A 點與 D 點連成的線為 0 則代表那個點是距離趨勢變化的點。

The Convex Hull Property of Bezier Curves

An important property of Bezier curves is that they always lie within the convex hull of their control points. The convex hull can be envisioned by pounding an imaginary nail into each control point, stretching an imaginary rubber band so that it surrounds the group of nails, and then collapsing that rubber band around the nails. The polygon created by that imaginary rubber band is the convex hull.

計算 _twin_spline_error_squared() 用這個特性，計算 B 與 C 點到 A 點與 D 點的連線的距離會大於真正的貝茲曲線到 A 點與 D 點的距離。因為貝茲曲線一定位為控制點之間。

就像Computer Aided Geometric Design 所示紅色的點也就是斜率與 P₀ 與 P₃ 點相連的線的斜率相同的點，都是最適合的 t 點

依照這個想法將他們依序再分貝茲曲線，會得到更好接近曲線的線段。以的三張左下角的說明，如果 t == 0.5，用原先的的方式左半邊計算的超過 tolerance，則繼續切割下去。為 t==0.25 此時左半邊的計算小於 tolerance ，則我們將其視為 t 點。

但其實對與整個曲線而言這不是一個好的分割方式，像是下方的 t == 0.33 會是更好的方式，既符合左半邊的不超過 tolerance ，右半邊也更接近 tolerance 。

也因此才有下面彈性 update 的方式，當找到某個點 t 的時候其左半邊的 tolerance 夠小，先繼續先將其存起來。(第 26 行)。並儲存與點 A 與點 D 連成的線的距離，而這個距離必定小與 tolerance_squared (第 24 行)。接下來基於這個 best 的值，繼續尋找其他 t 點與點 A 與點 D 連成的線的距離會更大，這是為了運用上面的特性:「最適合的點 t 這個的特質就是切線斜率會與這個線段，而這個點與點 A 與點 D 的距離會是最遠的」如同Precise Flattening of Cubic Bezier Segments論文所描述最適合的 t 的性質。因此如果沒有更大就不用找了(第 20 行)。而尋找這個方式用 Bisection 的方式，但並非如同more-iterative-splines 用一個 flat tolerance 的 threshold 去計算。

#define TWIN_SFIXED_TOLERANCE (TWIN_SFIXED_ONE >> 2)

static void _twin_spline_decompose(twin_path_t *path,
                                   twin_spline_t *spline,
                                   twin_dfixed_t tolerance_squared)
{
    /* Draw starting point */
    _twin_path_sdraw(path, spline->a.x, spline->a.y);

    while (!is_flat(spline, tolerance_squared)) {
        twin_dfixed_t hi = TWIN_SFIXED_ONE*TWIN_SFIXED_ONE, lo = 0, best =0, t;
        twin_dfixed_t best_error = 0, error;
        twin_spline_t left, right;

        while(true){
            t = (hi + lo) >> 1;
            if(best == t)
                break;

            _de_casteljau(spline, t, &left, &right);

            error = _twin_spline_error_squared(&left);
            if(error < best_error)
                break;

            /* Left is close enough to the bezier curve */
            if(error <= tolerance_squared ){
                best_error = error;
                best = t;

                /* Try to find better t */
                lo = t;
            }
            else{
                /* Find good t already and cannot find better t */
                if(best)
                    break;
                hi = t;
            }
        }

        _de_casteljau(spline, best, &left, &right);

        /* Draw the left segment */
        _twin_path_sdraw(path, left.d.x, left.d.y);

        /* Update spline to the right segment */
        memcpy(spline, &right, sizeof(twin_spline_t));
    }

    /* Draw the ending point */
    _twin_path_sdraw(path, spline->d.x, spline->d.y);
}

點比較少，但iteration 較多，這跟原先想的很不一樣，原先預期的結果會是演算法總是找到適合的 t ，這樣就如同，減少右半邊的的曲線範圍，這樣找到的 t 也會更容易符合 tolerance。這是一個先苦後甘的想法，因為要找的更適合的 t 要花更多的 iteration 次數是一定的，然而期待右邊的計算會因為前面找的 t 更適合而減少 iteration 的次數，但是實作的結果卻不是這樣，即便點能夠比較少，也只少了一點點， t 的值也就是從

$\frac{1}{8}$ 移到

$\frac{3}{8}$ 而已。而右邊也不會因為這樣的變化而受惠使 iteration 次數減少，反而右邊仍在繼續用更多的 iteration 次數尋找更好的 t 值。因此原先這種先苦後甘的想法不會存在，如果再犧牲計算以取得更好的 t 無法在有限步數看到計算量更小，則不用再繼續期待更久之後 iteration 次數會更少了。就像數學歸納法的道理一樣。

This is done by first proving a simple case, then also showing that if we assume the claim is true for a given case, then the next case is also true.

Precise Flattening of Cubic Bezier Segments論文所描述最適合的 t 的性質 的本質在於先計算出 t 為多所少最好，然後直接用這個 t 計算所在的 x 與 y。而上方彈性 update 的方式是透過先嘗試 t 然後計算這個 t 值的 x 與 y 值，再去透過與點 A 與點 D 的連線距離去判斷這個 t 值好不好。因為是先嘗試了，所以最終的結果在測試 font-edit 下 iteration 次數會變成兩倍。
實驗結果：

Experiment 2:
Flexible - Average number of _de_casteljau calls per point: 4.53
Flexible - Average points per character: 16.30

Computer Aided Geometric Design

float b_horner(float *b, int n, float t)
{
    float u, bc, tn, tmp;
    int i;
    u = 1.0 - t;
    bc = 1;
    tn = 1;
    tmp = b[0]*u;
    for(i=1; i<n; i++){
        tn = tn*t;
        bc = bc*(n-i+1)/i;
        tmp = (tmp + tn*bc*b[i])*u;
    }
    return (tmp + tn*t*b[n]);
}

版本三 : accurate distance

原本的方法是用下列的原理
source

/*
 * Return an upper bound on the distance (squared) that could result from
 * approximating a spline with a line segment connecting the two endpoints,
 * which is based on the Convex Hull Property of Bézier Curves: The Bézier Curve
 * lies completely in the convex hull of the given control points. Therefore, we
 * can use control points B and C to approximate the actual spline.
 */
static twin_dfixed_t _twin_spline_distance_squared(twin_spline_t *spline, twin_dfixed_t t)
{
    twin_spoint_t bc,cb;
    twin_dfixed_t bcdist, cbdist;
    _lerp(&spline->b, &spline->c, t, &bc);
    _lerp(&spline->d, &spline->c, t, &cb);
    bcdist  = _twin_distance_to_line_squared(&bc, &spline->a, &spline->d);
    cbdist  = _twin_distance_to_line_squared(&bc, &spline->a, &spline->d);

    if (bcdist > cbdist)
        return bcdist;
    return cbdist;
}

不能只算一個 bc 不然會出現抖動問題。

/*
 * Check if a spline is flat enough by comparing the distance against the
 * tolerance.
 */
static bool is_flat(twin_spline_t *spline, twin_dfixed_t tolerance_squared)
{
    twin_dfixed_t berr, cerr;

    berr = _twin_distance_to_line_squared(&spline->b, &spline->a, &spline->d);
    cerr = _twin_distance_to_line_squared(&spline->c, &spline->a, &spline->d);

    if (berr > cerr)
        return berr <= tolerance_squared ;
    return cerr <= tolerance_squared ;
}

points_comparison

dcj_calls_comparison

實驗結果

original - Average iterations per point: 1.99
original - Average points per character: 18.89
flexibly update - Average iterations per point: 4.53
flexibly update - Average points per character: 16.30
flexibly update w/ accurate distance - Average iterations per point: 4.40
flexibly update w/ accurate distance - Average points per character: 15.86

版本四

hint : 看@
原先 : Character '@' Points 54 Iteration 114
彈性 update : Character '@' Points 47 Iteration 230
accurate distance : Character '@' Points 47 Iteration 224

因為 @ 的數量幾乎是兩倍以上所以將 @ 畫出來看點的位置。

版本三 : accurate distance

想法 :

1. berr cerr consider both
2. 調整t
3. 用 around order 
4. 數值放大

實驗結果 :

original - Average iterations per point: 1.99
original - Average points per character: 18.89
flexibly update - Average iterations per point: 4.53
flexibly update - Average points per character: 16.30
flexibly update w/ accurate distance - Average iterations per point: 4.40
flexibly update w/ accurate distance - Average points per character: 15.86
original w/ accurate distance - Average iterations per point: 1.94
original w/ accurate distance - Average points per character: 17.88

貝茲曲線參考 - vector-graphics

vector-graphics/src/bezier.h

lerp3 lerp2 lerp1 三次貝茲曲線
cut2_base 的 u
u.x 是 a a 。 u.y 是 a b 。 u.z 是 b b 。
cut3_base 的 u 的 x y z w 是用不同的比例線性差值。
u.x 是 a a a 。 u.y 是 a a b 。 u.z 是 a b b 。 u.w 是 b b b 。

原本 lerp 3 會找出 x0 x1 x2 x3 的最終點。也就是 t = a 的方式。現在則是透過

inline real norm2(real x, real y){
    return x*x + y*y;
}

inline vec3 cross3(real x0, real y0, real w0, real x1, real y1, real w1){
    vec3 v = {y0*w1 - y1*w0, w0*x1 - w1*x0, x0*y1 - x1*y0};
    return v;
}

inline vec4 power3(real x0, real x1, real x2, real x3){
    vec4 v = {x0, 3.f*(x1-x0), 3.f*(x2-2.f*x1+x0), x3-x0-3.f*(x2-x1)};
    return v;
}
vec4 crosspmatrix(real x0, real y0, real x1, real y1, real x2, real y2, real x3, real y3){
    vec4 u = power3(x0, x1, x2, x3);
    vec4 v = power3(y0, y1, y2, y3);
    vec4 d;
    d.x = 0.f;
    d.y = det2(u.z, u.w, v.z, v.w);
    d.z = -det2(u.y, u.w, v.y, v.w);
    d.w = det2(u.y, u.z, v.y, v.z);
    return d;
}

分成 3 類

enum Cubic classify3(real x0, real y0, real x1, real y1, real x2, real y2, real x3, real y3){
    vec4 d = crosspmatrix(x0, y0, x1, y1, x2, y2, x3, y3);
    d.x = 3.f*d.z*d.z-4.f*d.y*d.w;
    if (!is_almost_zero(d.y)) {
        if(d.x > 0.f)
            return SERPENTINE;
        else if(d.x < 0.f)
            return LOOP;
        else
            return CUSP_WITH_INFLECTION_AT_INFINITY;
    }
    else if(!is_almost_zero(d.z))
        return CUSP_WITH_CUSP_AT_INFINITY;
    else if(!is_almost_zero(d.w))
        return DEGENERATED_QUADRATIC;
    else
        return LINE_OR_POINT;
}

//Receives the control points from a cubic bezier curve,
//then if the bezier is degenerated to a lower order curve
//modify the points to store only the necessary points
//and return the type of the curve
enum Cubic regenerate(real* pts){
    enum Cubic type = classify3(pts[0], pts[1], pts[2], pts[3], pts[4], pts[5], pts[6], pts[7]);
    if (type == LINE_OR_POINT) {
        pts[2] = pts[6]; pts[3] = pts[7];
    }
    else if (type == DEGENERATED_QUADRATIC) {
        vec3 d0, d1;
        cubic_tan(&d0, &d1, pts[0], pts[1], pts[2], pts[3], pts[4], pts[5], pts[6], pts[7]);
        vec3 p = cross3(d0.x, d0.y, d0.z, d1.x, d1.y, d1.z);
        
        pts[2] = p.x / p.z; pts[3] = p.y / p.z;
        pts[4] = pts[6]; pts[5] = pts[7];
    }
    return type;
}

貝茲曲線參考 - skia

src/base/SkBezierCurves.cpp

Optimal Curve Fitting to Digital Data

The curve technique has used various ideas for curve design. These ideas include end-point interpolation, intermediate point interpolation, detection of characteristic points, and parameterization.

The subdivision process, in this paper, is managed as follows.

Simplifying Bézier paths

---

TODO: Issue 貝茲曲線的頂點會抖動

會抖動的原因是因為線段比較粗。cap 的畫法問題。因為如果是 TwinCapRound 的話，就不會有這個問題出現，抖動的問題發生在 TwinCapButt, TwinCapProjecting。

static twin_dfixed_t _twin_spline_distance_squared(twin_spline_t *spline)
{
    twin_dfixed_t dist;
    twin_spoint_t bc;
    _lerp(&spline->a, &spline->b, (TWIN_SFIXED_ONE * TWIN_SFIXED_ONE) >> 1, &bc);
    return _twin_distance_to_line_squared(&bc, &spline->a, &spline->d);
}

這樣也會抖動，會抖動是因為頂點反轉了，所以兩種方向的頂點都需要被考慮，

---

TODO: Issue TwinCapButt 半圓

在 twin_path_arc() 會把圓弧的點的製作出來，且只在 pen 這個 path 上製作這些點，不會在 stroke 上。 pen 只紀錄一個圓從 0 度到 360 度，因為在 twin_composite_stroke() 函式中呼叫的函式 twin_path_circle(pen, 0, 0, pen_width / 2); 的x跟y都是 0。然後 pen 裡面只會存這一個以 pen_width / 2 為半徑的圓。

void twin_composite_stroke(twin_pixmap_t *dst,
                           twin_operand_t *src,
                           twin_coord_t src_x,
                           twin_coord_t src_y,
                           twin_path_t *stroke,
                           twin_fixed_t pen_width,
                           twin_operator_t operator)
{
    twin_path_t *pen = twin_path_create();
    twin_path_t *path = twin_path_create();
    twin_matrix_t m = twin_path_current_matrix(stroke);

    m.m[2][0] = 0;
    m.m[2][1] = 0;
    twin_path_set_matrix(pen, m);
    twin_path_set_cap_style(path, twin_path_current_cap_style(stroke));
    twin_path_circle(pen, 0, 0, pen_width / 2);
    twin_path_convolve(path, stroke, pen);
    twin_composite_path(dst, src, src_x, src_y, path, operator);
    twin_path_destroy(path);
    twin_path_destroy(pen);
}

Step 1 : _apps_spline_paint()
Step 2 :twin_paint_stroke()
Step 3 : twin_composite_stroke()
Step 3.a twin_path_convolve()
Step 3.b twin_composite_path() -> twin_composite()

Step 2 :twin_paint_stroke()

畫黑色的線

    // This should be TwinCapButt
    twin_paint_stroke(_apps_spline_pixmap(spline), 0xff404040, path, //畫黑色的線
                      twin_int_to_fixed(50)); 
    twin_path_set_cap_style(path, TwinCapButt);
    twin_paint_stroke(_apps_spline_pixmap(spline), 0xffffff00, path, //畫黃色的線
                      twin_int_to_fixed(2));

Step 3.a twin_path_convolve()

convolve.c
twin_path_convolve() 先呼叫 twin_path_convex_hull() 再呼叫 _twin_subpath_convolve

twin_path_convex_hull()

sort 的比較函式：_twin_hull_vertex_compare()

/* Compare two slopes. Slope angles begin at 0 in the direction of the
   positive X axis and increase in the direction of the positive Y
   axis.

   WARNING: This function only gives correct results if the angular
   difference between a and b is less than PI.

   <  0 => a less positive than b
   == 0 => a equal to be
   >  0 => a more positive than b
*/
static int _twin_slope_compare(twin_slope_t *a, twin_slope_t *b)
{
    twin_dfixed_t diff;
    diff = ((twin_dfixed_t) a->dy * (twin_dfixed_t) b->dx -
            (twin_dfixed_t) b->dy * (twin_dfixed_t) a->dx);

a b
[1,1] [-1,-2] -> -1 - (-2) =1
[1,2] [1,1] -> 1

_twin_subpath_convolve()

typedef enum _twin_cap {
    TwinCapRound,
    TwinCapButt,
    TwinCapProjecting,
} twin_cap_t;

line 黑色的是 cap projecting 中間的是 Butt
spline 兩個都是 Butt
source

5655b42

        case TwinCapButt:
            p = ptarget - 1;
            while (p != ptarget) {
                if (++p == np)
                    p = 0;
                _twin_path_sdraw(path, sp[s].x + pp[p].x, sp[s].y + pp[p].y);
            }
            break;
            // 但會出現半個白圓
            /* fall through … */
        case TwinCapRound:
            while (p != ptarget) {
                if (++p == np)
                    p = 0;
                _twin_path_sdraw(path, sp[s].x + pp[p].x, sp[s].y + pp[p].y);
            }
            break;

Graham scan

半圓形不是來自這個，不會建立新的path

有些點的著色被移除了

matrix 用在字旋轉上

將 twin_path_convolve() 一次換一個顏色，包含不同的 subpath 透過 _twin_subpath_convolve() (很少看到subpath數量

$\ge$ 2 )，其中將 draw last 的部份分為
TwinCapRound

draw to starget + draw last	draw last

TwinCapButt

draw to starget + draw last	draw last

TwinCapProjecting

draw to starget + draw last	draw last

TwinCapButt - case 2
origin

draw to starget + draw last	draw last

從這二張圖看出來有問題的應該為畫 draw to starget 的那一段。因為它多出了

$(x,y)=(2000,-1500)$ 的那一個半圓與

$(x,y)=(4500,-1500)$ 的半圓。因為那一段使得藍色的圈圈不會被塗到。

紅色的點對應到原先黃色的線，讓黑色的spline 出現是用橘色的點。但原本spline 透過四個點計算出來就只有紅色的點。再用 twin_path_circle() 把整個線段變粗。

void twin_path_circle(twin_path_t *path,
                      twin_fixed_t x,
                      twin_fixed_t y,
                      twin_fixed_t radius)
{
    twin_path_ellipse(path, x, y, radius, radius);
}
void twin_path_ellipse(twin_path_t *path,
                       twin_fixed_t x,
                       twin_fixed_t y,
                       twin_fixed_t x_radius,
                       twin_fixed_t y_radius)
{
    twin_path_move(path, x + x_radius, y);
    twin_path_arc(path, x, y, x_radius, y_radius, 0, TWIN_ANGLE_360);
    twin_path_close(path);
}

TwinCapProject - case 2

static void _twin_subpath_convolve(twin_path_t *path,
                                   twin_path_t *stroke,
                                   twin_path_t *pen)
{
    twin_spoint_t *sp = stroke->points;
    twin_spoint_t *pp = pen->points;
    int ns = stroke->npoints;
    int np = pen->npoints;
    twin_spoint_t *sp0 = &sp[0];
    twin_spoint_t *sp1 = &sp[1];
    // 統一從最左邊的畫。
    int start = _twin_path_leftpoint(pen, sp0, sp1); 

每次都只畫這個圓 pn 與 pm 為畫 pen 所存的圓當前點 p 的下一個點與前一個點。如果 pn 也就是下一個點等於圓的盡頭也就是 np - 1 則代表畫完了。如果所有點 pn 與 pm 與 stroke 也就是畫筆的骨幹斜率一樣。則代表該畫筆的骨幹位置的圓已經畫完，則移動到下一個畫筆的骨幹位置 s = sn (第17行)

        /*
         * Convolve the edges
         */
        do {
            int sn = s + inc;
            int pn = (p == np - 1) ? 0 : p + 1;
            int pm = (p == 0) ? np - 1 : p - 1;

            /*
             * step along pen (forwards or backwards) or stroke as appropriate
             */
            if (_around_order(&sp[s], &sp[sn], &pp[p], &pp[pn]) > 0) {
                p = pn;
            } else if (_around_order(&sp[s], &sp[sn], &pp[pm], &pp[p]) < 0) {
                p = pm;
            } else {
                s = sn;
            }
            // printf("%d\t%d\n",sp[s].x + pp[p].x, sp[s].y + pp[p].y);
            _twin_path_sdraw(path, sp[s].x + pp[p].x, sp[s].y + pp[p].y);
        } while (s != starget);

解決辦法:找出butt 最後的水平線，把線段另一側的點都刪掉

m = (sp1->y - sp0->y) / (sp1->x - sp0->x)
-1/m = (sp0->x - sp1->x) / (sp1->y - sp0->y)
Y - sp0->y = ((sp0->x - sp1->x) / (sp1->y - sp0->y)) (X - sp0->x)
(sp1->y - sp0->y)(Y - sp0->y) = (sp0->x - sp1->x)(X - sp0->x)
0 = (sp0->x - sp1->x)X - (sp0->x - sp1->x)(sp0->x) - (sp1->y - sp0->y)Y + (sp1->y - sp0->y)(sp0->y)
0 = (sp0->x - sp1->x)X - (sp1->y - sp0->y)Y - (sp0->x - sp1->x)(sp0->x) + (sp1->y - sp0->y)(sp0->y)
0 = (sp0->x - sp1->x)X + (sp0->y - sp1->y)Y + (sp1->x - sp0->x)(sp0->x) + (sp1->y - sp0->y)(sp0->y)

    twin_dfixed_t A1 = (sp0->x - sp1->x);
    twin_dfixed_t B1 = (sp0->y - sp1->y);
    twin_dfixed_t C1 = (twin_dfixed_t) (sp1->x - sp0->x) * (sp0->x) - 
                        (twin_dfixed_t) (sp1->y - sp0->y) * (sp0->y);

    twin_sfixed_t A1 = (path->spline->a.x - path->spline->b.x);
    twin_sfixed_t B1 = (path->spline->a.y - path->spline->b.y);
    twin_sfixed_t C1 = twin_dfixed_mul((path->spline->b.x - path->spline->a.x),path->spline->a.x) - 
                        twin_dfixed_mul((path->spline->b.y - path->spline->a.y),path->spline->a.y);

cairo-path-stroke
_cairo_stroker_init
_cairo_stroker_limit -> _cairo_boxes_get_extents() 每個邊界加上 stroke 的寬度
_tessellate_fan

目前mado 沒有這個功能

To test whether we need to join two segments of a spline using a round-join or a bevel-join

pen 跟 mado 的作法一樣先畫一個圓。

	_cairo_slope_init (&v->slope_cw, &prev->point, &v->point);
	_cairo_slope_init (&v->slope_ccw, &v->point, &next->point);

兩線段相連處會看是順時針還逆時針

CapButt
cairo-polygon.c
_cairo_polygon_add_external_edge()

cairo_status_t
_cairo_polygon_add_external_edge (void *polygon,
				  const cairo_point_t *p1,
				  const cairo_point_t *p2)
{
    _cairo_polygon_add_edge (polygon, p1, p2, 1);
    return _cairo_polygon_status (polygon);
}

在 cairo-path-stroke.c 傳進來函式的是 closure

cairo_status_t (*add_external_edge) (void *closure,
					 const cairo_point_t *p1,
					 const cairo_point_t *p2);

解決辦法:判斷鈍角三角形

bool triangle(twin_spoint_t *A,
                twin_spoint_t *B,
                twin_sfixed_t cx, twin_sfixed_t cy){
    twin_sfixed_t ax = A->x;
    twin_sfixed_t ay = A->y;
    twin_sfixed_t bx = B->x;
    twin_sfixed_t by = B->y;
    if( ax == bx && ay == by)
        return false;
    twin_dfixed_t AB = twin_dfixed_mul(ax - bx, ax - bx) + twin_dfixed_mul(ay - by, ay - by);
    twin_dfixed_t BC = twin_dfixed_mul(bx - cx, bx - cx) + twin_dfixed_mul(by - cy, by - cy);
    twin_dfixed_t CA = twin_dfixed_mul(cx - ax, cx - ax) + twin_dfixed_mul(cy - ay, cy - ay);
    if (AB + BC <= CA || BC + CA <= AB || CA + AB <= BC)
        return true;
    return false;
}

如果是鈍角三角形就不在路徑內

poly-isect 解決路徑自我重疊。

水彩

            if (TwinCapButt == path->state.cap_style){
                twin_sfixed_t dx = -(sp[sn].y -sp[s].y);
                twin_sfixed_t dy = sp[sn].x -sp[s].x;
                twin_sfixed_t ddx = twin_dfixed_div(twin_dfixed_mul(pen_half_width,dx),_twin_sfixed_sqrt(twin_square_add(dx,dy)));
                twin_sfixed_t ddy = twin_dfixed_div(twin_dfixed_mul(pen_half_width,dy),_twin_sfixed_sqrt(twin_square_add(dx,dy)));
                printf("%d\t%d\n",sp[s].x + dx, sp[s].y + dy);  
                _twin_path_sdraw(path, sp[s].x + dx, sp[s].y + dy);
                if (s == 0){
                    s = sn;
                    continue;
                }
                s = sn;
            } 

鋸齒狀

            if (TwinCapButt == path->state.cap_style){
                twin_sfixed_t up_y = sp[sn].x -sp[s].x;
                twin_sfixed_t up_x = -(sp[sn].y -sp[s].y);
                twin_sfixed_t down_y = -(sp[sn].x -sp[s].x);
                twin_sfixed_t down_x = sp[sn].y -sp[s].y;
                
                _twin_path_sdraw(path, sp[s].x + up_x, sp[s].y + up_y);
                if (s == 0){
                    s = sn;
                    continue;
                }
                up_y = sp[s].x -sp[s-1].x;
                up_x = -(sp[s].y -sp[s-1].y);
                down_y = -(sp[s].x -sp[s-1].x);
                down_x = sp[s].y -sp[s-1].y;
                printf("%d\t%d\n",sp[s].x + up_x, sp[s].y + pp[p].y);  
                _twin_path_sdraw(path, sp[s].x + up_x, sp[s].y + up_y);
                s = sn;

半徑有問題

            if (TwinCapButt == path->state.cap_style){
                twin_sfixed_t dx = -(sp[sn].y -sp[s].y);
                twin_sfixed_t dy = sp[sn].x -sp[s].x;
                twin_sfixed_t c = _twin_sfixed_sqrt(twin_square_add(dx,dy));
                if(c == 0){
                    if (dx == 0){
                        printf("%d\t%d\n",sp[s].x, sp[s].y + pen_half_width);  
                        _twin_path_sdraw(path, sp[s].x, sp[s].y + pen_half_width);
                    }
                    else{
                        printf("%d\t%d\n",sp[s].x + pen_half_width, sp[s].y);  
                        _twin_path_sdraw(path, sp[s].x + pen_half_width, sp[s].y);
                    }
                    s = sn;
                    continue;
                }
                twin_sfixed_t ddx = twin_dfixed_div(twin_dfixed_mul(pen_half_width,dx),_twin_sfixed_sqrt(twin_square_add(dx,dy)));
                twin_sfixed_t ddy = twin_dfixed_div(twin_dfixed_mul(pen_half_width,dy),_twin_sfixed_sqrt(twin_square_add(dx,dy)));
                _twin_path_sdraw(path, sp[s].x + ddx, sp[s].y + ddy);
                s = sn;
            }

offset

座標 (2000, -1800~-1400) 的地方 offset 還是會畫出來

如果 _around_order(&sp[s], &sp[sn], &pp[p], &pp[pn]) < 0 且 _around_order(&sp[s], &sp[sn], &pp[pm], &pp[p]) > 0 就是用原先 pen 做好的圓上兩個相鄰點的向量相減就是平行於路徑。如果用路徑的斜率 m 的垂直 (-1/m) ，計算出來的半徑與 twin_path_circle() 計算出的不一樣。

#define twin_dfixed_mul(a, b) \
    ((twin_sfixed_t) (((twin_dfixed_t) (a) * (twin_dfixed_t) (b)) >> 8))
#define twin_square_add(a ,b) \ 
    ((twin_dfixed_mul(a,a)) + twin_dfixed_mul(b,b))

    twin_sfixed_t pen_half_width = _twin_sfixed_sqrt(twin_square_add(pp[0].x,pp[0].y));

            if (path->state.cap_style == TwinCapButt){
                if ( _around_order(&sp[s], &sp[sn], &pp[p], &pp[pn]) > 0 ) {
                    p = pn;
                } else if ( _around_order(&sp[s], &sp[sn], &pp[pm], &pp[p]) < 0 ) {
                    p = pm;
                } else {
                    s = sn;
                    _twin_path_sdraw(path, sp[s].x + pp[p].x, sp[s].y + pp[p].y);
                }

Cramer's rule

後來我發現 offset 的特性就是他的移動路徑應該要跟原先的曲線一樣。所以如果移動路徑不相同則代表有缺口。

因此透過計算有沒有交點可以消掉缺口

            if (TwinCapButt == path->state.cap_style){
                if ( _around_order(&sp[s], &sp[sn], &pp[p], &pp[pn]) > 0 ) {
                    p = pn;
                } else if ( _around_order(&sp[s], &sp[sn], &pp[pm], &pp[p]) < 0 ) {
                    p = pm;
                } else {
                    s = sn;
                    twin_sfixed_t ax = sp[prev_s].x;
                    twin_sfixed_t ay = sp[prev_s].y;
                    twin_sfixed_t bx = sp[prev_s].x + pp[prev_p].x;
                    twin_sfixed_t by = sp[prev_s].y + pp[prev_p].y;
                    twin_sfixed_t cx = sp[s].x;
                    twin_sfixed_t cy = sp[s].y;
                    twin_sfixed_t dx = sp[s].x + pp[p].x;
                    twin_sfixed_t dy = sp[s].y + pp[p].y;
  
                    //Cramer's rule
                    twin_dfixed_t t_num = (cx-ax)*(dy-cy) - (cy-ay)*(dx-cx);
                    twin_dfixed_t t_den = (bx-ax)*(dy-cy) - (by-ay)*(dx-cx);
                    twin_dfixed_t u_num = (cx-ax)*(by-ay) - (cy-ay)*(bx-ax);
                    twin_dfixed_t u_den = t_den;
                    if(!((  ((t_num < 0 && t_den < 0 )  && (-t_num < -t_den))||
                            ((t_num > 0 && t_den > 0 )  && (t_num < t_den))) &&
                         (  ((u_num < 0 && u_den < 0 )  && (-u_num < -u_den))||
                            ((u_num > 0 && u_den > 0 )  && (u_num < u_den))))){
                        if(firstdraw)
                            _twin_path_smove(path, sp[s].x + pp[p].x, sp[s].y + pp[p].y);
                        else{
                            _twin_path_sdraw(path, sp[s].x + pp[p].x, sp[s].y + pp[p].y);
                            firstdraw = false;
                        }
                    }
                    prev_s = s;
                    prev_p = p;   
                    firstdraw = false;
                }
            }

剩下開頭無法解決，把後面的 butt 判斷式都拿掉

        case TwinCapButt:
            p = ptarget - 1;
+           break;
            /* fall through … */
        case TwinCapRound:
            while (p != ptarget) {
                if (++p == np)
                    p = 0;
                _twin_path_sdraw(path, sp[s].x + pp[p].x, sp[s].y + pp[p].y);
            }
            break;
        }

convolve.c#L178
路徑折返

        /* reach the end of the path?  Go back the other way now */
        inc = -1;
        ptarget = start;
        starget = 0;

因為會路徑折返，所以判斷 firstdraw 初始化要寫在迴圈最外面。

沒拿掉就會從新設一個起點，將折反點視為路徑的起點。

            if (path->state.cap_style == TwinCapButt){
                if ( _around_order(&sp[s], &sp[sn], &pp[p], &pp[pn]) > 0 ) {
                    p = pn;
                } else if ( _around_order(&sp[s], &sp[sn], &pp[pm], &pp[p]) < 0 ) {
                    p = pm;
                } else {
                    s = sn;          
+                    if(prev_s != -1){   
                        twin_sfixed_t ax = sp[prev_s].x;
                        twin_sfixed_t ay = sp[prev_s].y;
                        twin_sfixed_t bx = sp[prev_s].x + pp[prev_p].x;
                        twin_sfixed_t by = sp[prev_s].y + pp[prev_p].y;
                        twin_sfixed_t cx = sp[s].x;
                        twin_sfixed_t cy = sp[s].y;
                        twin_sfixed_t dx = sp[s].x + pp[p].x;
                        twin_sfixed_t dy = sp[s].y + pp[p].y;
                                
                        //Cramer's rule
                        twin_dfixed_t t_num = (cx-ax)*(dy-cy) - (cy-ay)*(dx-cx);
                        twin_dfixed_t t_den = (bx-ax)*(dy-cy) - (by-ay)*(dx-cx);
                        twin_dfixed_t u_num = (cx-ax)*(by-ay) - (cy-ay)*(bx-ax);
                        twin_dfixed_t u_den = t_den;
                        if(!((  ((t_num < 0 && t_den < 0 )  && (-t_num < -t_den))||
                                ((t_num > 0 && t_den > 0 )  && (t_num < t_den))) &&
                             (  ((u_num < 0 && u_den < 0 )  && (-u_num < -u_den))||
                                ((u_num > 0 && u_den > 0 )  && (u_num < u_den))))){
                            if(firstdraw){
                                printf("%d\t%d\t%d\n",index++,sp[s].x + pp[p].x, sp[s].y + pp[p].y);
                                _twin_path_smove(path, sp[s].x + pp[p].x, sp[s].y + pp[p].y);
                            }
                            else{
                                printf("%d\t%d\t%d\n",index++,sp[s].x + pp[p].x, sp[s].y + pp[p].y);
                                _twin_path_sdraw(path, sp[s].x + pp[p].x, sp[s].y + pp[p].y);
                            }
+                           firstdraw = false;
                        }
                    }
                    prev_s = s;
                    prev_p = p;  
-                   firstdraw = false;
                }
            }

        case TwinCapButt:
-           p = ptarget - 1;
+           p = ptarget;
            /* fall through … */

GPU-friendly Stroke Expansion

GPU-friendly Stroke Expansion demo

The way of Render Strokes

1. local techniques that break down the stroke into individual closed primitives.
2. distance field techniques (or similar techniques based on point inclusion).

This paper focus on stroke expansion, the generation of an outline that, when filled, represents the rendered stroke.
Implementing stroke rendering using stroke expansion enables a unified approach to rendering both filled and stroked paths downstream. Among other benefits, such an approach avoids a large number of draw calls when stroked and filled paths are finely interleaved, as is often the case.

strong correctness

the computation of the outline of a line swept along the segment, maintaining normal orientation, combined with stroke caps and joins.

Weak correctness

only requires the parallel curves (also known as “offset curves”) of path segments

Solution

Converting stroked primitives to filled primitives

https://github.com/diegonehab/stroke-to-fill/blob/e54cde26fea5f72645838a48ffcb9fff4c204596/rvg-cubic-bezier-evolute-approximator.h#L130

auto evolute = make_cubic_bezier_evolute_approximator<
RVG_CUBIC_BEZIER_APPROXIMATION_SAMPLES>(m_max_radius, s, ds, dds, m_sink);
auto ni = (ri/len(di))*perp(di);
auto qi = pi + ni;
auto nf = (rf/len(df))*perp(df);
auto qf = pf + nf;
m_sink.linear_segment(pi, qi);
evolute.approximate_partition(ti, qi, perp(di), tf, qf, perp(df),
m_splits, m_tol);
m_splits.clear();
m_sink.linear_segment(qf, pf);

stroke 
|- path
|- style
|- tolerance
|- strong correct (opts)

ctx
|- strong correct (opts)
|- join_thresh
|- tolerance
|- output
|- forward_path
|- backward_path
|- start_point
|- start_norm
|- start_tan
|- last_pt: Default::default(),
|- last_tan: Default::default(),

Euler Spirals

        for es in CubicToEulerIter::new(c, es_tol) {
            self.forward_path
                .do_euler_seg(&es, -0.5 * style.width, lower_tol, self.opts.strong);
            self.backward_path
                .do_euler_seg(&es, 0.5 * style.width, lower_tol, self.opts.strong);
        }

evolute

1. 在走訪 outlier 的時候，先判斷三個點
   $nod{e}_{i-1}$
   $nod{e}_i$
   $nod{e}_{i+1}$ 形成的為鈍角三角形，如果為鈍角三角形繼續走訪。
2. 如果是銳角三角形且
   $nod{e}_a$ 不存在，就先紀錄當前的點
   $nod{e}_i$ 為
   $nod{e}_a$。回到第一步。
3. 如果是銳角三角形且
   $nod{e}_a$ 存在，就先紀錄當前的點
   $nod{e}_i$ 為
   $nod{e}_b$。
4. 將
   $nod{e}_c$ =
   $nod{e}_{(a+b)/2}$
5. 開始進行 evolute。
6. 透過
   $nod{e}_a$ 與
   $nod{e}_c$ 算其中垂線。 再透過
   $nod{e}_b$ 與
   $nod{e}_c$ 算其中垂線，中垂線的交點為圓心再計算半徑 radius。
7. 圓心，
   $nod{e}_a$，
   $nod{e}_b$，
   $nod{e}_c$，旋轉角度 rotate_angle = angle -30
   

因為希望整個弧長座落在 240 度到 360 度，無論原先弧長的長度，因此要移動圓心，圓心不會總是 stroke 路徑上的點。

舉例：

8. 移動圓心的方式：
   

• 先計算移動距離：

  $nod{e}_a$，
  $nod{e}_b$ 中點
  $nod{e}_{mid}$直接移動
  $nod{e}_a$，
  $nod{e}_b$距離的
  $\frac{\sqrt{3}}{2}$

• 計算單位向量：（dx,dy)=（中點

  $nod{e}_{mid}$ 到 圓心的向量）/ （中點
  $nod{e}_{mid}$ 到 圓心的距離）
  

• 計算移動向量：單位向量*移動距離

• $nod{e}_{mid}$ + 移動向量 = 新的圓心。

9. deltoid_curve 從 start_angle = 0,end_angle = 240

10. $t\in$[start_angle,end_angle],
    
    $x=$ radius
    $\times \frac{2\times \cos (t)}{3}+$ radius
    $\times \frac{\cos (2\times t)}{3}$
    
    $y=$ radius
    $\times \frac{2\times \sin (t)}{3}-$ radius
    $\times \frac{\sin (2\times t)}{3}$

11. 旋轉
    

    $x=x\times \cos ($ rotate_angle
    $)-y\times \sin ($ rotate_angle
    $)$
    
    $y=y\times \sin ($ rotate_angle
    $)+y\times \cos ($ rotate_angle
    $)$

12. 位移
    

    $x=x+$ center_x
    
    $y=y+$ center_y

twin_spoint_t* find_center(twin_spoint_t *a, 
                twin_spoint_t *b, 
                twin_spoint_t *mid) {
    
    twin_dfixed_t slope1, intercept1, slope2, intercept2;
    twin_spoint_t center;
    line_parameters(a, b, &slope1, &intercept1);
    line_parameters(b, mid, &slope2, &intercept2);

    if (slope1 == NULL) {
        center.x = twin_dfixed_to_sfixed(intercept1);
        center.y = twin_dfixed_to_sfixed(
                        twin_dfixed_mul(slope2,twin_sfixed_to_dfixed(center.x)) 
                        + intercept2);
    } else if (slope2 == NULL) {
        center.x = twin_dfixed_to_sfixed(intercept2);
        center.y = twin_dfixed_to_sfixed(
                        twin_dfixed_mul(slope1,twin_sfixed_to_dfixed(center.y))
                        + intercept1);
    } else {        
        center.x = twin_dfixed_to_sfixed(
                        twin_dfixed_div(twin_sfixed_to_dfixed(intercept2 - intercept1),
                                        twin_sfixed_to_dfixed(slope1 - slope2)));
        center.y = twin_dfixed_to_sfixed(   
                        twin_dfixed_mul(slope1,twin_sfixed_to_dfixed(center.x))
                        + intercept1);
    }
    return &center;
}

---

TODO: Implement drop shadows using blurring techniques

X 中的混成器與 Composite 擴展

Add StackBlur for imgproc

Guassian Blur

premultiply_alpha

[ Pixmap: 10x10x1000000 ]
argb32_source_argb32                     0.103462 sec
argb32_over_argb32                       0.174524 sec
[ Pixmap: 100x100x20000 ]
argb32_source_argb32                     0.108871 sec
argb32_over_argb32                       0.177337 sec
[ Pixmap: 200x200x10000 ]
argb32_source_argb32                     0.141014 sec
argb32_over_argb32                       0.342229 sec
[ Pixmap: 1200x800x200 ]
argb32_source_argb32                     0.174455 sec
argb32_over_argb32                       0.991179 sec

first, instead of doing 2D blurring in one go, do it in two passes: one horizontal, one vertical. This works because the second time we blur, each pixel already has the average of all pixels in the other direction.

uint8_t kernel[5]= {0x01, 0x04, 0x06, 0x04, 0x01}; // 0x16
水平 kenerl[i] * pixel >> 4 = new pixel map
垂直 new pixel map * kernel >> 4 = target pixel map

uint8_t kernel[7] = {0x01, 0x06, 0x0f, 0x14, 0x0f, 0x06, 0x01}; // 0x64
如果是 0x64 每次都要 >> 6 ，原先每個為 0xff * kernel[i]
在 >> 8 以內，都行。
所以最多可以
uint8_t kernel[9] = {0x01, 0x08, 0x1c, 0x38, 0x46, 0x38, 0x1c, 0x08, 0x01}
kernel 9 sum is 256.
Therefore, the sum of kernel map is 256*256 = 2^16

[ Pixmap: 10x10x1000000 ]
argb32_source_argb32                     0.111397 sec
argb32_over_argb32                       0.191449 sec
[ Pixmap: 100x100x20000 ]
argb32_source_argb32                     0.152879 sec
argb32_over_argb32                       0.216831 sec
[ Pixmap: 200x200x10000 ]
argb32_source_argb32                     0.151917 sec
argb32_over_argb32                       0.358300 sec
[ Pixmap: 1200x800x200 ]
argb32_source_argb32                     0.384055 sec
argb32_over_argb32                       1.474838 sec

if (radius == 2){ // 2^8
        *r += ((((v & 0x00ff0000) >> 16) * (twin_argb32_t) wght) << 8) & 0x00ff0000;
        *g += (((v & 0x0000ff00) >> 8) * (twin_argb32_t) wght) & 0x0000ff00;
        *b += ((((v & 0x000000ff) >> 0) * (twin_argb32_t) wght) >> 8) & 0x000000ff;
    }
    else if(radius == 3){ // 2^12
        *r += ((((v & 0x00ff0000) >> 16) * (twin_argb32_t) wght) << 10) & 0x00ff0000;
        *g += (((v & 0x0000ff00) >> 8) * (twin_argb32_t) wght << 2) & 0x0000ff00;
        *b += ((((v & 0x000000ff) >> 0) * (twin_argb32_t) wght) >> 6) & 0x000000ff;
    }
    else if(radius == 4){ // 2^16
        *r += ((((v & 0x00ff0000) >> 16) * (twin_argb32_t) wght)) & 0x00ff0000;
        *g += (((v & 0x0000ff00) >> 8) * (twin_argb32_t) wght >> 8) & 0x0000ff00;
        *b += ((((v & 0x000000ff) >> 0) * (twin_argb32_t) wght) >>16) & 0x000000ff;
    }
    return;

Multiple_buffering
Compositing_window_manager Wiki
策略一 ： stacking

the window manager tells the back window to repaint itself

策略二 ： compositing

the window manager maintains an off-screen memory buffer containing the full appearance of each window, including the back window.

With a stacking manager, the repainting process can become corrupted when a program that is slow, unresponsive or buggy does not respond to messages in a timely manner. A malicious program can cause the system to appear unstable by simply neglecting to repaint its window.

物件 A
A = malloc
A -> twin_toplevel_t = twin_toplevel_create(twin_screen_t)
twin_toplevel_t
twin_window_t -> client_data : toplevel

With the window system supporting a single screen containing many windows, the toolkit extends this model by creating a single top-level widget.
This top-level widget contains a single box for layout purposes.
Each box can contain a number of widgets or other boxes.
Layout within each box is done either horizontally or vertically with an algorithm which comes from the Layout Widget[3] that the author developed for Xt[1] library.
Each widget has a natural size and stretch in both directions.
The natural size and stretch of a box is computed from the objects it contains.
This forms the sole geometry management mechanism within the toolkit and is reasonably competent at both constructing a usable initial layout and adapting to externally imposed size changes

single screen
|- many windows
|- single top-level widget
    |- single box
        |- a number of widgets 
        |- or other boxes.

Input would run in one thread, events were dispatched without queuing directly to the receiving
object.

Alpha and the History of Digital Compositing
Real-Time Rendering
Interpreting Alpha

邊界處理

void twin_screen_update(twin_screen_t *screen)

這段程式碼很厲害，它只要更新 damge 的部份而已，而不是整個 screen 去更新。但這樣就沒有一整塊最上層的 pixel map 讓我去讀，然後複製到陰影上再做模糊。

---

Setting

 target-y :=
 target.o-y :=
-TARGET_LIBS :=
+TARGET_LIBS := -lpthread

Blurring and Sharpening Images

四層pixelmap, 較遠的pixelmap移動的較少 較近的pucelmap移動的較多 ，這樣就有遠近的感覺
parallex 視差

距離遠的變化少
距離短的變化大
content level

1. SBS resolution 差
2. TB
3. Interleave
   frame 跟 frame 不會動的就不用重畫 ，要利用這樣的資訊去省掉pixel 腳踏車跟人也不用重畫只有移動。只有近的東西有動，但腳踏車還是一樣的不用重畫只是copy paste。
   紀錄當下的左眼差。
   只要研究delta = frame 跟frame的差值。跟twin screen update一樣。

3D projection 到的東西都要是2D矩形
跟世界地圖展開的方式一樣。保持畫面的連續性，這樣就可以用delta 去 render 這樣才會省硬體的資源。如果都隨便render到不同地方。畫一隻鳥飛過去 下一個frame就要重畫。
視野範圍是錐形,送了一堆pixel給user ，user只能看到一點點而已。field of view。

dynamic adaptive stream
content level
projections
transportation level
view points
看要給哪個視野版本給別人看
現實上要做到不一樣。

view point 改善計算成本。

twin.h

    /*
     * Repaint function
     */
    twin_put_begin_t put_begin;
    twin_put_span_t put_span;
    void *closure;

put_begin backend 會用到。

/*
 * twin_put_begin_t: called before data are drawn to the screen
 * twin_put_span_t: called for each scanline drawn
 */
typedef void (*twin_put_begin_t)(twin_coord_t left,
                                 twin_coord_t top,
                                 twin_coord_t right,
                                 twin_coord_t bottom,
                                 void *closure);
typedef void (*twin_put_span_t)(twin_coord_t left,
                                twin_coord_t top,
                                twin_coord_t right,
                                twin_argb32_t *pixels,
                                void *closure);

twin.h pixelmap 有時候可以代表 window 。

/*
 * A rectangular array of pixels
 */
typedef struct _twin_pixmap {
    /*
     * When representing a window, this point
     * refers to the window object
     */
    twin_window_t *window;
} twin_pixmap_t;

因此判斷該pixelmap 是不是 window 的最底層的方式可以直接檢查 pixmap->windown == NULL。所以在 twin_window_create() 建立完最底層的 pixmap 後會再指回去。
window->pixmap = twin_pixmap_create(format, width, height);
window->pixmap->window = window;

取子視窗的第幾的 row 。
twin_coord_t p_y = y % screen->background->height;

    for (p_left = left; p_left < right; p_left += p_this) {
        dst.argb32 = span + (p_left - left);
        m_left = p_left % p_width;
        p_this = p_width - m_left;
        if (p_left + p_this > right)
            p_this = right - p_left;
        src.p = twin_pixmap_pointer(screen->background, m_left, p_y);
                    bop32(dst, src, p_this);
    }

p_this 是範圍。

    dst.argb32 = span + 0;
    m_left = left % p_width;
    p_this = p_width - m_left;
    if (left + p_this > right)
        p_this = right - left;
    src.p =
        twin_pixmap_pointer(screen->background, m_left, p_y);
    bop32(dst, src, p_this);

如果把 for 迴圈拿掉結果是一樣。

原先 drop_shadow 的程式碼

會將背景所有的 pixel map 都 render 在最上層的 pixel map 的 shadow 範圍。

static void twin_drop_shadow(twin_screen_t *screen,
                             twin_pixmap_t *active_pix,
                             twin_pixmap_t *prev_active_pix)
{
    twin_pixmap_t *p = NULL;
    twin_pointer_t dst;
    twin_source_u src;
    twin_coord_t start, end, overlap, y, src_y, ori_wid, ori_hei, tgt_start,
        _tgt_start, tgt_end, src_start, src_end;

    twin_src_op pop32 = _twin_argb32_over_argb32,
                bop32 = _twin_argb32_source_argb32;

    /* Remove the drop shadow from the previously active pixel map. */
    if (prev_active_pix) {
        src.c = 0x00000000;
        for (y = 0; y < prev_active_pix->height; y++) {
            if (y < prev_active_pix->height - prev_active_pix->window->shadow_y)
                twin_cover(
                    prev_active_pix, src.c,
                    prev_active_pix->width - prev_active_pix->window->shadow_x,
                    y, prev_active_pix->window->shadow_x);
            else
                twin_cover(prev_active_pix, src.c, 0, y,
                           prev_active_pix->width);
        }
        prev_active_pix->shadow = false;
    }

    /* Mark the previously active pixel map as damaged to update its changes. */
    if (prev_active_pix && active_pix != prev_active_pix)
        twin_pixmap_damage(prev_active_pix, 0, 0, prev_active_pix->width,
                           prev_active_pix->height);

    /*
     * The shadow effect of the window only becomes visible when the window is
     * active.
     */
    if (active_pix) {
        active_pix->shadow = true;
        ori_wid = (*active_pix).width - (*active_pix).window->shadow_x;
        ori_hei = (*active_pix).height - (*active_pix).window->shadow_y;
        tgt_start = (*active_pix).x;
        tgt_end = (*active_pix).x + (*active_pix).width;
        for (y = 0; y < (*active_pix).height; y++) {
            /*
             * Take the screen's background as the bottom layer of the shadow.
             */
            if (screen->background) {
                p = screen->background;
                src_y = (*active_pix).y + y;
                if (src_y > 0 && src_y < (*p).height - 1) {
                    if (y < ori_hei)
                        _tgt_start = tgt_start + ori_wid;
                    else
                        _tgt_start = tgt_start;

                    overlap = tgt_end - _tgt_start;
                    src_start = _tgt_start;

                    dst = twin_pixmap_pointer(active_pix,
                                              _tgt_start - tgt_start, y);
                    src.p = twin_pixmap_pointer(p, src_start, src_y);
                    bop32(dst, src, overlap);
                }
            }
            /*
             * Render each layer of the pixmap under the shadow pixel map onto
             * the shadow.
             */
            for (p = screen->bottom; p; p = p->up) {
                /* Only render pixel maps beneath the visible drop shadow. */
                if (p->shadow)
                    break;

                /*
                 * Identify the areas where the currently pixel map overlaps
                 * with the drop shadow pixel map.
                 */
                src_y = (*active_pix).y + y - (*p).y;
                if (src_y < 0 || src_y >= (*p).height - 1)
                    continue;

                if (y < ori_hei)
                    _tgt_start = tgt_start + ori_wid;
                else
                    _tgt_start = tgt_start;

                src_start = (*p).x;
                src_end = (*p).x + (*p).width;

                if (src_start > tgt_end || _tgt_start > src_end)
                    continue;
                start = max(src_start, _tgt_start);
                end = min(src_end, tgt_end);
                overlap = end - start;
                dst = twin_pixmap_pointer(active_pix, start - tgt_start, y);
                src.p = twin_pixmap_pointer(p, start - src_start, src_y);
                pop32(dst, src, overlap);
            }
        }
        /*
         * Create a darker border of the active window that gives a more
         * dimensional appearance.
         */

        /* Create a temporary pixmap to handle the shadow function. */
        twin_pixmap_t *tmp_px = twin_pixmap_create(
            active_pix->format, active_pix->width, active_pix->height);
        tmp_px->window = active_pix->window;

        /* The shift offset and color of shadow can be selected by user. */
        twin_shadow_border(tmp_px, 0xff000000, 1, 1);

        /* Add a frosted glass effect to the shadow of the active window. */
        /* Right side of the active window */
        twin_stack_blur(tmp_px, 10, ori_wid,
                        ori_wid + active_pix->window->shadow_x, 0, ori_hei);
        /* Bottom side of the active window */
        twin_stack_blur(tmp_px, 10, 0, active_pix->width, ori_hei,
                        ori_hei + active_pix->window->shadow_y);

        for (twin_coord_t y = 0; y < active_pix->height; y++) {
            dst = twin_pixmap_pointer(active_pix, 0, y);
            src.p = twin_pixmap_pointer(tmp_px, 0, y);
            pop32(dst, src, active_pix->width);
        }
        twin_pixmap_destroy(tmp_px);
    }
}

---

twin_sdl_poll
裡面
twin_screen_dispatch

參考

glhull
handwriter.ttf
StrokeStyles: Stroke-Based Segmentation and Stylization of Fonts

Pixelmap

Pixel-art_scaling_algorithms
Depixelizing Pixel Art