2020q3 Homework10 (quiz 10)

contributed by < JimmyLiu0530 >

tags: `進階電腦系統理論與實作`

測驗1: RGB to Gray

給定每個像素 (pixel) 為 32-bit 的 RGBA 的位元圖 (bitmap)，RGBA 代表紅綠藍三原色的首字母，Alpha 值則表示顏色的透明度/不透明度，其轉換為灰階影像 (gray scale) 的函式為：


















void rgba_to_bw(uint32_t *bitmap, int width, int height, long stride) {
    int row, col;
    uint32_t pixel, r, g, b, a, bw;
    
    for (row = 0; row < height; row++) {
        for (col = 0; col < width; col++) {
            pixel = bitmap[col + row * stride / 4];
            
            a = (pixel >> 24) & 0xff;
            r = (pixel >> 16) & 0xff;
            g = (pixel >> 8) & 0xff;
            b = pixel & 0xff;
            bw = (uint32_t) (r * 0.299 + g * 0.587 + b * 0.114);
            
            bitmap[col + row * stride / 4] = (a << 24) + (bw << 16) + (bw << 8) + (bw);
        }
    }
}

人類視網膜包含三類錐形對光敏感的細胞，分別對不同波長範圍的光有反應:

Image Not Showing Possible Reasons

The image file may be corrupted
The server hosting the image is unavailable
The image path is incorrect
The image format is not supported

Learn More →

上圖中的縱軸代表人眼對某對應波長的光的敏感度。其中 1 Å (埃) =
$10^{- 10}$ 米 =
$0.1 n m$ [出處]

人眼吸收綠色比其他顏色敏感，也可說人眼最容易捕捉到綠色，所以當影像變成灰階時，僅僅將紅色、綠色、藍色加總取平均，不足以反映出人眼所見。常見的方法是將

R e d \times 77, G r e e n \times 151, B l u e \times 28

，這三個除數的總和為 256 (即

2^{8}

)，可使除法變簡單 (等同於右移 8 位元)。

如果我們將所有的運算都建立表格，於是

bwPixel = table[rgbPixel & 0x00ffffff] + rgbPixel & 0xff000000;

\to

16 MB; 表格太大

如果先計算針對「乘上 0.299」一類的運算，計算後再建立表格呢？

bw = (uint32_t) mul_299[r] + (uint32_t) mul_587[g] + (uint32_t) mul_144[b];
bitmap[col + row * strike / 4] = (a << 24) + (bw << 16) + (bw << 8) + bw;

\to

降到 32 KB 以內; cache friendly

接著嘗試使用 NEON 指令集加速。將 RGB 三色的權重存入 r3, r4, r5 等暫存器。

vdup.8 (Vector Duplicate)，分別複製到大小為 8 bit 的 NEON register d0 - d2

    mov     r3, #77
    mov     r4, #151
    mov     r5, #28
    vdup.8  d0, r3
    vdup.8  d1, r4
    vdup.8  d2, r5

vld4.8 (Vector Load)，載入 pixel 的資料到 4 個 8-bit 的 NEON register d4-d7，其中那個 4 為 interleave，因為我們有 ARGB，所以 gap = 4。

再來就是計算 weighted average。可用 Vector Multiply 和 Vector Multiply Accumulate 指令:

   @ (alpha,R,G,B) = (d7,d6,d5,d4)
   vld4.8      {d4-d7}, [r0]!
   vmull.u8    q10, d6, d0
   vmlal.u8    q10, d5, d1
   vmlal.u8    q10, d4, d2

將值除以 256 就是我們要的灰階值。

vrshrn (Vector Shift Right by immediate value)

vrshrn.u16 d4, q10, #8

最後儲存結果。

vst (Vector Store)

vst4.8 {d4-d7}, [r3]!

效能評比如下: (對照 v0 和 v4)

Image Not Showing Possible Reasons

The image file may be corrupted
The server hosting the image is unavailable
The image path is incorrect
The image format is not supported

Learn More →

從這張表即可理解，使用 NEON 指令加速只用原本 4.6％的時間即可完成 RGBA 轉灰階的處理。

Neon Intrinsics 允許我們避免直接撰寫組合語言，而是有如 C 函式呼叫，例如 vmull.u8 和 vmlal.u8 分別對應以下 intrinsics: (提示: 在 Neon Intrinsics 網頁搜尋 NEON intrinsics 可獲得更詳盡的描述)

考慮以下使用 float 來保存像素的程式碼:
























#include <stdint.h>
#define MIN(a, b) (((a) < (b)) ? (a) : (b))
#define MAX(a, b) (((a) > (b)) ? (a) : (b))

uint8_t *to_rgb(const float *m, int w, int h, uint8_t *rgb)
{
    int size = w * h;
    const float *ptr0 = m + (0 * size);
    const float *ptr1 = m + (1 * size);
    const float *ptr2 = m + (2 * size);
    uint8_t *ret = rgb;

#define SATURATE_CAST_UCHAR(X) (uint8_t) MIN(MAX((int) (X), 0), 255);
    for (int remain = size; remain > 0; remain--) {
        rgb[0] = SATURATE_CAST_UCHAR(*ptr0);
        rgb[1] = SATURATE_CAST_UCHAR(*ptr1);
        rgb[2] = SATURATE_CAST_UCHAR(*ptr2);
        rgb += 3;
        ptr0++, ptr1++, ptr2++;
    }
#undef SATURATE_CAST_UCHAR

    return ret;
}

函式 to_rgb 可將原本透過連續的 float 空間保存的資料轉為 RGB 像素，每個像素佔用 3 個位元組，即 24 位元。針對 Aarch64 (Arm 64-bit) 架構 (假設為 little-endian)，下方的程式碼則可將 RGB 轉換為灰階:







































#include <arm_neon.h>
float *rgb_to_gray(const uint8_t *rgb, int w, int h, float *m)
{
    const uint8_t Y_shift = 8;
    const uint8_t R2Y = 77, G2Y = 151, B2Y = 28;
    float *ptr = m;
    int size = w * h;

    int nn = size >> 3;
    int remain = size - (nn << 3);

    uint8x8_t _R2Y = vdup_n_u8(R2Y);
    uint8x8_t _G2Y = vdup_n_u8(G2Y);
    uint8x8_t _B2Y = vdup_n_u8(B2Y);

    for (; nn > 0; nn--) {
        uint8x8x3_t _rgb = vld3_u8(rgb);

        uint16x8_t y16 = vmull_u8(_rgb.val[0], _R2Y);
        y16 = vmlal_u8(y16, _rgb.val[1], _G2Y);
        y16 = vmlal_u8(y16, _rgb.val[2], _B2Y);
        y16 = vshrq_n_u16(y16, Y_shift);

        float32x4_t _ylow = vcvtq_f32_u32(vmovl_u16(vget_low_u16(y16)));
        float32x4_t _yhigh = vcvtq_f32_u32(vmovl_u16(vget_high_u16(y16)));

        vst1q_f32(ptr + OFFSET1, _ylow);
        vst1q_f32(ptr + OFFSET2, _yhigh);

        rgb += 3 * 8, ptr += 8;
    }

    for (; remain > 0; remain--) {
        *ptr = (rgb[0] * R2Y + rgb[1] * G2Y + rgb[2] * B2Y) >> Y_shift;
        rgb += 3, ptr++;
    }

    return m;
}

請補完程式碼，使上述程式碼得以符合預期地運作。

程式說明

line 4-7 初始化等等會用到的變數，像是 RGB 的權重、位移量…等。

因為本題是在 64-bit 的環境下執行，因此暫存器也是 64-bit，代表一輪可以處理 8 個像素。因此， nn 紀錄需要做的輪數，剩餘不足一輪的會在 line 33-35 單獨處理。

line 12-14 將 RGB 各顏色的權重分別載入到 64-bit 暫存器，其中載入的方法為複製 8 份 1-byte 權重(e.g. R2Y、G2Y、B2Y)到暫存器中。

每輪執行以下步驟:

先將 8 個像素一次載入到三個暫存器。

    uint8x8x3_t _rgb = vld3_u8(rgb);

載入的方式是採 interleaving，也就是說從 rgb 取到的第一個 byte 會載入到第一個暫存器，取到的第二個 byte 會載入到第二個暫存器，取到的第三個 byte 會載入到第三個暫存器，取到的第四個 byte 會載入到第一個暫存器。

計算出灰階值

    uint16x8_t y16 = vmull_u8(_rgb.val[0], _R2Y);
    y16 = vmlal_u8(y16, _rgb.val[1], _G2Y);
    y16 = vmlal_u8(y16, _rgb.val[2], _B2Y);
    y16 = vshrq_n_u16(y16, Y_shift);

vmull_u8 會將第一個 argument 與第二個 argument 的各個 byte 相乘，因此乘出來的結果就會變成 16 bits。全部的結果就會變成 16 x 8 = 128 bits。vmulalu8 就只是 vmull_u8 乘完再與另一個暫存器相加而已。
最後，將計算出來的加權值除以 256，也就是向右位移 8 後，即為灰階值。

存回結果

    float32x4_t _ylow = vcvtq_f32_u32(vmovl_u16(vget_low_u16(y16)));
    float32x4_t _yhigh = vcvtq_f32_u32(vmovl_u16(vget_high_u16(y16)));

    vst1q_f32(ptr + OFFSET1, _ylow);
    vst1q_f32(ptr + OFFSET2, _yhigh);

    rgb += 3 * 8, ptr += 8;

現在要將每個像素的 RGB 從原本 16-bit unsigned 轉成 32-bit float，因為我們原本就是使用 float 來保存資料。
轉換的方法如下:

將原本存於 y16 中 8 個像素的資料分成兩群各自去做轉換，因此 _ylow、_yhigh 裏頭各自有 4 個像素的資料。
最後將 _ylow、_yhigh 存回 ptr 中。顯然，_ylow 要存回的初始位置就是 ptr，所以OFFSET1 = 0。因為 ylow 佔 32*4 bits，也就是 4 個 float，所以 _yhigh 要從 ptr + 4 開始存，因此 OFFSET1 = 4。

延伸問題

1. 整合 libattopng，輸出 RGBA 和灰階的 PNG 圖片

可能會用到 ImageMagick 事先處理輸入檔案

參照 Bomb Lab 實作紀錄得知如何安裝 GNU Toolchain 及 QEMU

測驗2: LeetCode 1617. Count Subtrees With Max Distance Between Cities

給定

n

個城市，編號從

1

到

n

，再給予一個大小為

n - 1

的陣列 edges，其中

e d g e s [i] = [u_{i}, v_{i}]

表示城市

u_{i}

和

v_{i}

之間存在雙向邊。題目保證任兩個城市間有唯一的路徑，也就是說所有城市構成一棵樹 (tree)。

上圖的聯通關係為:

1, 4

1, 3

3, 2

1, 3, 2

4, 1, 3

4, 1, 3, 2

一棵子樹 (subtree) 是上述所有城市的子集。在子集中，任意城市間可透過子集的其他城市和邊到達。兩個子樹被認為不一樣的條件是，至少有個城市在其中一棵子樹中存在，卻又不存在於另一棵子樹中。

給定

d

範圍介於

1

到

n - 1

，請找出城市間最大距離恰好為

d

的所有子樹數目。請你回傳一個大小為

n - 1

的陣列，其中第

d

個元素（下標從 1 開始）是城市間最大距離恰好等於

d

子樹數目。

兩個城市間距離定義是，它們之間需要經過的邊的數目。

第一組輸入:

$n = 4$
$e d g e s = [[1, 2], [2, 3], [2, 4]]$

預期輸出：

$[3, 4, 0]$

既然
$n = 4$ ，於是輸出的陣列會有 3 個元素

子樹 {1,2}, {2,3} 和 {2,4} 最大距離都為 1

子樹 {1,2,3}, {1,2,4}, {2,3,4} 和 {1,2,3,4} 最大距離都為 2

不存在城市間最大距離為 3 的子樹。

第二組輸入:

$n = 2$
$e d g e s = [[1, 2]]$

預期輸出：

$[1]$

既然
$n = 2$ ，於是輸出的陣列只會有 1 個元素，且題目已保證任意城市之間只有唯一的路徑

第三組輸入:

$n = 3$
$e d g e s = [[1, 2], [2, 3]$

預期輸出：

$[2, 1]$

既然
$n = 3$ ，於是輸出的陣列會有 2 個元素

子樹 {1,2} 和 {2,3} 最大距離都為 1

子樹 {1,2,3} 最大距離為 2

限制:

$2 \leq n \leq 15$
$e d g e s . l e n g t h == n - 1$
$e d g e s [i] . l e n g t h == 2$
$1 \leq u_{i}, v_{i} \leq n$
題目保證
$(u_{i}, v_{i})$ 所表示的邊互不相同

思路:

因為
$n \leq 15$ ，所以就往狀態壓縮 + dynamic programming 的技巧去思考，但本題狀態間的轉移無法進行，只能窮舉所有情況。
題目的樹可看成是無向無環連通圖（二者等價）

在圖論中，環是一條只有第一個和最後一個頂點重複的非空路徑

不具備環的圖被稱作無環圖

不具備有向環的有向圖被稱做有向無環圖 (DAG,directed acyclic graph)

無環的連通圖被稱作樹

窮舉狀態，首先要判斷是否無效
1. 該狀態僅有
  $1$ 個節點，因題目要求的
  $d$ 最小是
  $1$ ，即距離為
  $1$ ，所以可縮減
2. 狀態下不是一個子樹，即非一個連通圖，那怎麼判斷是否是連通圖呢？可利用 BFS 或 DFS 判斷所能到達的節點的個數，是否和集合內元素個數 (二進位
  $1$ 的個數)相同即可
判斷該狀態是有效後，如何求最大距離？所謂的最大距離是無向無環連通圖的最大直徑

以下程式碼是可能的解法:






























































#include <stdint.h>
#include <stdlib.h>
#include <string.h>

#define MIN(a, b) (((a) < (b)) ? (a) : (b))

#define ITERATE(i, begin, end) \
    for (int i = (begin), i##_end_ = (end); i < i##_end_; i++)

int *countSubgraphsForEachDiameter(int n,
                                   int **edges,
                                   int edgesSz,
                                   int *edges0sz,
                                   int *rsz)
{
    uint8_t d[n][n]; /* distance */
    memset(d, 0xFF, sizeof(d));
    uint16_t conn[n];
    ITERATE (i, 0, n)
        conn[i] = 1 << i;

    for (int z = 0; z < edgesSz; z++) {
        int *e = edges[z];
        int u = e[0] - 1, v = e[1] - 1;
        d[u][v] = d[v][u] = 1;
        conn[u] |= 1 << v, conn[v] |= 1 << u;
    }

    ITERATE (k, 0, n)
        ITERATE (i, 0, n)
            ITERATE (j, 0, n)
                d[i][j] = MIN(d[i][j], d[i][k] + d[k][j]);
    int *rv = calloc(n - 1, sizeof(int));

    ITERATE (S, 1, (1 << n)) {
        int ctmp = 1 << __builtin_ctz(S);
        while (1) {
            int conn_nxt = ctmp;
            ITERATE (d, 0, n)
                if ((ctmp >> d) & 1)
                    conn_nxt |= conn[d];
            conn_nxt &= S;
            if (COND1)
                break;
            ctmp = conn_nxt;
        }
        if (COND2)
            continue;
        uint8_t ret = 0;
        ITERATE (i, 0, n)
            if ((S >> i) & 1)
                ITERATE (j, 0, i)
                    if (((S >> j) & 1) && (ret < d[i][j]))
                        ret = d[i][j];
        if (!ret)
            continue;
        rv[ITER]++;
    }

    *rsz = n - 1;
    return rv;
}

參考資訊:

程式說明

首先，介紹一些變數的意義:

d[i][j]:節點 i 與節點 j 之間的最短距離
conn[i]: 與節點 i 相鄰的點，用二進位表示

根據 input，將 conn、d 做初始化，然後計算兩點間的最短距離。

上面提到，我們採用窮舉所有子樹來求解，不過我們要先去檢查該子樹是否為連通? 作法說明如下:

S 代表一顆子樹，從這子樹中 index 最小的點開始，不斷新增與此點相鄰的點進 conn_nxt，因為我們只在乎這顆子樹的點而已，所以其他不在這子樹的點我們都忽略，對應到程式碼

    conn_nxt &= S;

透過不斷的從現有的點中往外擴展，直到收斂，因此，COND1 為 conn_nxt == ctmp。
此時，若該子樹為連通，則 ctmp 一定包含該子樹所有的點。反之，如果
ctmp 沒有包含該子樹所有的點，也就是說 ~ctmp & S != 0，則該子樹不連通，因此 COND2 == ~ctmp & S。

接著確認是連通圖後，針對該子樹去找出這棵樹的 diameter，即

m a x_{u, v \in V} d [u] [v]

，其中

V

為該子樹的點集合。找法說明如下:

看過所有該子圖中的點與子圖中其他點的距離，就可以找出最大的距離，即 diameter。
最後，根據該子樹的 diameter 去遞增相對應 rv[] 的位置，因此，ITER == --ret。

指出改進方式
更有效率找 all pair shortest path 的方法-BFS:
$O (n^{3})$ to
$O (n^{2})$

2020q3 Homework10 (quiz 10)

tags: 進階電腦系統理論與實作

測驗1: RGB to Gray

程式說明

延伸問題

1. 整合 libattopng，輸出 RGBA 和灰階的 PNG 圖片

測驗2: LeetCode 1617. Count Subtrees With Max Distance Between Cities

程式說明

延伸問題

1. 用 C99/C11 (可搭配 GNU extension) 重新實作符合 LeetCode 要求但效率優於上述表現的實作，一併分析時間複雜度。

Read more

台灣萬萬稅

CS:APP 第一章重點整理 - Introduction

GDB介紹

Ubuntu 20.04 wifi 選項消失之解決方法

tags: `進階電腦系統理論與實作`