Halide Tutorial 非官方中文翻譯 - Part 4

Halide Tutorial 非官方中譯 - Part 1
Halide Tutorial 非官方中譯 - Part 2
Halide Tutorial 非官方中譯 - Part 3

==

Lesson 9 - Multi-pass Funcs, update definitions, and reductions

// Halide tutorial lesson 9:


// 宣告後續會使用的變數
Var x("x"), y("y");

// 讀取灰階影像做為輸入影像
Buffer input = load_image("images/gray.png");

// 能夠定義 multiple  pass Func. 首先定義玩具範本
{

// 第一的定義如同先前所看過的
// - 一個自 Var 到 Expr 的 mapping
Func f;
f(x, y) = x + y;
// 我們稱這個定義為"純"(pure)定義

// 但之後的定義包含在等號兩邊需要計算的表示式. 最簡單的範例是只修改一個點
f(3, 7) = 42;

// 我們稱這些額外的定義為"更新"(update)定義, 或是"歸納"(reduction)定義.
// 一個歸納定義是一個遞迴地參考自身同位置目前數值的更新定義
f(x, y) = f(x, y) + 17;

// 若限縮更新範圍在單一橫列, 就能遞迴參考在同一直行的數值
f(x, 3) = f(x, 0) * f(x, 10);
// 同樣地, 若限縮更新範圍在單一直行, 就能遞迴參考在同一橫列的數值
f(0, y) = f(0, y) / f(3, y);

// 通則是: 每個用在一個更新定義中的變數, 必須在等號兩邊出現在純定義中對應不變的位置// 因此下列皆為合法的定義
f(x, 17) = x + 8;
f(0, y) = y * 8;
f(x, x + 1) = x + 8;
f(y/2, y) = f(0, y) * 17;

// 但下列都會產生錯誤

// f(x, 0) = f(x + 1, 0);
// 右邊的 f 的第一個參數應為 x 而非 x+1

// f(y, y + 1) = y + 8;
// 左邊 f 的第二個參數應為 y 而非 y+1

// f(y, x) = y - x;
// 左邊的 f 的參數位置錯誤

// f(3, 4) = x + y;
// 右邊出現了變數, 然而左邊並有

// 實現這函數僅是確認它編譯過,
// 從第二個之後的定義將強迫在長大於寬的區域上實現
f.realize(100, 101);

// 對每個 f 的實現而言, 每一步都在下一個開始之前完成.
// 接著透過追蹤簡單例子的讀取與寫入
Func g("g");
g(x, y) = x + y;   // Pure definition
g(2, 1) = 42;      // First update definition
g(x, 0) = g(x, 1); // Second update definition

g.trace_loads();
g.trace_stores();

g.realize(4, 4);
// Click to show output ...

// 下圖為視覺化結果


// 讀取 log 後可以觀察到每個 pass 是依序填入
// 等義 C code 如下:
int result[4][4];
// Pure definition
for (int y = 0; y < 4; y++) {
    for (int x = 0; x < 4; x++) {
        result[y][x] = x + y;
    }
}
// 第一個更新定義
result[1][2] = 42;
// 第二個更新定義
for (int x = 0; x < 4; x++) {
    result[0][x] = result[1][x];
}

}

// 將更新定義放入 loops
{

// 首先用此純定義
Func f;
f(x, y) = (x + y)/100.0f;

// 接著想要能更新前 50 列
// 當然能夠直接加 50 個更新定義

// f(x, 0) = f(x, 0) * f(x, 0);
// f(x, 1) = f(x, 1) * f(x, 1);
// f(x, 2) = f(x, 2) * f(x, 2);
// ...
// f(x, 49) = f(x, 49) * f(x, 49);

// 或等義於在 C++ 程式碼加入在編譯時期的 loop
// for (int i = 0; i < 50; i++) {
//   f(x, i) = f(x, i) * f(x, i);
// }

// 但更具可管理性與更加彈性的方式是將這 loop 放入編譯所產生的程式碼中
// 這方式是藉由定義一個 "reduction domain" 並且在一個更新定義中使用所達到
RDom r(0, 50);
f(x, r) = f(x, r) * f(x, r);
Buffer<float> halide_result = f.realize(100, 100);

// 下圖為視覺化結果


 Your browser does not support the video tag :(

// 等義的 C code:
float c_result[100][100];
for (int y = 0; y < 100; y++) {
    for (int x = 0; x < 100; x++) {
        c_result[y][x] = (x + y)/100.0f;
    }
}
for (int x = 0; x < 100; x++) {
    for (int r = 0; r < 50; r++) {
        // reduction domain 的 loop 發生在任何純定義中作為更新步驟變數的 loop之中
        c_result[r][x] = c_result[r][x] * c_result[r][x];
    }
}

// 檢查結果是否一致
for (int y = 0; y < 100; y++) {
    for (int x = 0; x < 100; x++) {
        if (fabs(halide_result(x, y) - c_result[y][x]) > 0.01f) {
            printf("halide_result(%d, %d) = %f instead of %f\n",
                   x, y, halide_result(x, y), c_result[y][x]);
            return -1;
        }
    }
}

}

// 接著以實際用途來應用更新定義: 計算 histogram
{

// 一些影像上的操作無法單以自輸出座標到所寫入數值的純函數所表示
// 經典的例子是 histogram 的計算.
// 最直接的方式是不斷地在輸入影像上, 更新 histogram bucket
// 在 Halide 中如此實作:
Func histogram("histogram");

// 所以的 Histogram buckets 預設為 0
histogram(x) = 0;

// 在輸入影像上定義一個多維的 reduction domain
RDom r(0, input.width(), 0, input.height());

// 對於 reduction domain 的每個點, 將對應輸入影像點的數值的 histogram bucket 增加
histogram(input(r.x, r.y)) += 1;

Buffer halide_result = histogram.realize(256);

// 等義的 C code:
int c_result[256];
for (int x = 0; x < 256; x++) {
    c_result[x] = 0;
}
for (int r_y = 0; r_y < input.height(); r_y++) {
    for (int r_x = 0; r_x < input.width(); r_x++) {
        c_result[input(r_x, r_y)] += 1;
    }
}

// 檢查結果
for (int x = 0; x < 256; x++) {
    if (c_result[x] != halide_result(x)) {
        printf("halide_result(%d) = %d instead of %d\n",
               x, halide_result(x), c_result[x]);
        return -1;
    }
}

}

// 更新步驟的排程
{

// 如同以往, 更新步驟中的純變數依然可以 parallelized, vectorized 與 split

// 對 reduction domain 中的變數作 vectorizing, splitting 與 parallelize  較為棘手
// 後續的章節會說明

// 考量以下定義:
Func f;
f(x, y) = x*y;
// 將第0列的數值設為第8列的數值
f(x, 0) = f(x, 8);
// 將第1行的數值設為第8行的數值加上2
f(0, y) = f(8, y) + 2;

// 在每個 stage 中的純變數能獨立地被排程.
// 為了控制純定義, 能如同先前般地排程.
// 下列的程式碼對純定義做 vectorization 與 parallelization
f.vectorize(x, 4).parallel(y);

// 為了排程目的這裡使用 Func::update(int)來取得更新步驟所用的 handle
// 接著沿著 x 作 vectorization.
// 由於更新定義沒有使用到 y 變數, 這裡無法對 y 作操作
f.update(0).vectorize(x, 4);

// 對於第二更新步驟作以 4 做分割大小的平行化
Var yo, yi;
f.update(1).split(y, yo, yi, 4).parallel(yo);

Buffer halide_result = f.realize(16, 16);

// 下圖為視覺化結果

 Your browser does not support the video tag :(

// 以下為等義的 C code:
int c_result[16][16];

// 純定義的步驟, x  方向做 vectorization, 於 y 方向做 parallelization
for (int y = 0; y < 16; y++) { // Should be a parallel for loop
    for (int x_vec = 0; x_vec < 4; x_vec++) {
        int x[] = {x_vec*4, x_vec*4+1, x_vec*4+2, x_vec*4+3};
        c_result[y][x[0]] = x[0] * y;
        c_result[y][x[1]] = x[1] * y;
        c_result[y][x[2]] = x[2] * y;
        c_result[y][x[3]] = x[3] * y;
    }
}

// 第一次更新, 沿著 x 方向做 vectorization
for (int x_vec = 0; x_vec < 4; x_vec++) {
    int x[] = {x_vec*4, x_vec*4+1, x_vec*4+2, x_vec*4+3};
    c_result[0][x[0]] = c_result[8][x[0]];
    c_result[0][x[1]] = c_result[8][x[1]];
    c_result[0][x[2]] = c_result[8][x[2]];
    c_result[0][x[3]] = c_result[8][x[3]];
}

// 第二次更新, 以分割大小為 4 的方式, 沿著 y 做 parallelization
for (int yo = 0; yo < 4; yo++) { // Should be a parallel for loop
    for (int yi = 0; yi < 4; yi++) {
        int y = yo*4 + yi;
        c_result[y][0] = c_result[y][8] + 2;
    }
}

// 確認 C 與 Halide 結果一致:
for (int y = 0; y < 16; y++) {
    for (int x = 0; x < 16; x++) {
        if (halide_result(x, y) != c_result[y][x]) {
            printf("halide_result(%d, %d) = %d instead of %d\n",
                   x, y, halide_result(x, y), c_result[y][x]);
            return -1;
        }
    }
}

}

// 這些僅涵蓋如何在單個 Func 的更新步驟中的變數做排程
// 但像是 producer-consumer 牽涉到 compute_at 與 store_at 的關係又如何?
// 接著透過 producer-consumer 成對的方式, 來確認以 producer 作歸納的方式
{

// 由於一個更新是在儲存陣列上做的 multiple pass
// 因此將其作為 inline 沒有意義. 因此預設排程作可能上最為接近的動作
// 是在 consumer 最內層中計算這些數值
// 考慮此無關緊要的例子:
Func producer, consumer;
producer(x) = x*2;
producer(x) += 10;
consumer(x) = 2 * producer(x);
Buffer<int> halide_result = consumer.realize(10);

// 以下為視覺化結果




// 等義的 C code 如下:
int c_result[10];
for (int x = 0; x < 10; x++)  {
    int producer_storage[1];
    // 純定義計算 producer 的步驟
    producer_storage[0] = x * 2;
    // producer 的更新步驟
    producer_storage[0] = producer_storage[0] + 10;
    // 純定義計算 consumer 的步驟
    c_result[x] = 2 * producer_storage[0];
}

// 確認結果一致
for (int x = 0; x < 10; x++) {
    if (halide_result(x) != c_result[x]) {
        printf("halide_result(%d) = %d instead of %d\n",
               x, halide_result(x), c_result[x]);
        return -1;
    }
}

// 對於其他 compute_at/store_at 選項,
// 歸納會置於 consumer 巢狀 loop 中某處預期的位置

}

// 接著考量在 producer-consumer 組合中, 以 consumer 做歸納
// 這會有些複雜
{

{
    // Case 1: The consumer 在純定義步驟中參考 producer
    Func producer, consumer;
    // 這裡 producer 也是純定義.
    producer(x) = x*17;
    consumer(x) = 2 * producer(x);
    consumer(x) += 50;

    // 在這 case 中有效的 producer 排程包含預設排程 - inline 以及:    //
    // 1) producer.compute_at(x), 將 producer 計算置於 consumer 純定義中 x loop 內
    //    // 2) producer.compute_root(), 事先計算所有的 producer
    //
    // 3) producer.store_root().compute_at(x), 在 consumer 做外層的 loop 配置空間
    // 但在所需要數值的 loop 中填入
    //
    // 這裡使用選項 1

    producer.compute_at(consumer, x);

    Buffer<int> halide_result = consumer.realize(10);
    // 下圖為視覺化結果

    

    // 等義 C code:
    int c_result[10];
    // 計算 consumer 的純定義步驟
    for (int x = 0; x < 10; x++)  {
        // producer 的純定義步驟
        int producer_storage[1];
        producer_storage[0] = x * 17;
        c_result[x] = 2 * producer_storage[0];
    }
    // consumer 的更新步驟
    for (int x = 0; x < 10; x++) {
        c_result[x] += 50;
    }

    // 所有的純定義步驟都在更新步驟前計算
    // 所以會有著兩個在 x 方向上的 loop

    // 確認結果一致
    for (int x = 0; x < 10; x++) {
        if (halide_result(x) != c_result[x]) {
            printf("halide_result(%d) = %d instead of %d\n",
                   x, halide_result(x), c_result[x]);
            return -1;
        }
    }
}

{
    // Case 2: consumer 僅在更新步驟中參考 producer
    Func producer, consumer;
    producer(x) = x * 17;
    consumer(x) = 100 - x * 10;
    consumer(x) += producer(x);    // 相同地, 在 consumer 每個 x 座標中計算 producer
    // producer 程式碼會置於consumer 更新步驟中
    // 因為那是唯一使用到 producer 的步驟
    producer.compute_at(consumer, x);    // 注意, 這並不表示:
    //
    // producer.compute_at(consumer.update(0), x).
    //
    // 排程是藉由 Func 中的 Vars 所達成
    // 而 Func 內的 Vars 會在純定義與更新步驟中所共用

    Buffer<int> halide_result = consumer.realize(10);

    // 下圖為視覺化結果


    // 等義的 C code:
    int c_result[10];
    // consumer 的純定義步驟
    for (int x = 0; x < 10; x++)  {
        c_result[x] = 100 - x * 10;
    }
    // consumer 的更新步驟
    for (int x = 0; x < 10; x++) {
        // producer 的純定義步驟
        int producer_storage[1];
        producer_storage[0] = x * 17;
        c_result[x] += producer_storage[0];
    }


    // 確認結果一致
    for (int x = 0; x < 10; x++) {
        if (halide_result(x) != c_result[x]) {
            printf("halide_result(%d) = %d instead of %d\n",
                   x, halide_result(x), c_result[x]);
            return -1;
        }
    }
}

{
    // Case 3: consumer 在使用共變數的多個步驟中參考使用 producer
    Func producer, consumer;
    producer(x) = x * 17;
    consumer(x) = 170 - producer(x);
    consumer(x) += producer(x)/2;
    // 同樣地, 在 consumer 中每個 x 座標計算 producer
    // 這會把 producer 程式碼置於 consumer 的純定義與更新步驟
    // 因此最後的結果是兩個個別實現, 並有著重複的冗工
    producer.compute_at(consumer, x);

    Buffer<int> halide_result = consumer.realize(10);

    // 下圖為視覺化結果


    // 等義的 C code:
    int c_result[10];
    // consumer 的純定義步驟
    for (int x = 0; x < 10; x++)  {
        // producer 的純定義步驟
        int producer_storage[1];
        producer_storage[0] = x * 17;
        c_result[x] = 170 - producer_storage[0];
    }
    // consumer 的更新步驟
    for (int x = 0; x < 10; x++) {
        // 複製 producer 純定義步驟
        int producer_storage[1];
        producer_storage[0] = x * 17;
        c_result[x] += producer_storage[0]/2;
    }

    // 確認結果一致
    for (int x = 0; x < 10; x++) {
        if (halide_result(x) != c_result[x]) {
            printf("halide_result(%d) = %d instead of %d\n",
                   x, halide_result(x), c_result[x]);
            return -1;
        }
    }
}

{
    // Case 4: consumer 在不使用共用變數的多個步驟參考使用 producer
    Func producer, consumer;
    producer(x, y) = (x * y) / 10 + 8;
    consumer(x, y) = x + y;
    consumer(x, 0) = producer(x, x);
    consumer(0, y) = producer(y, 9-y);    // 在這個 case 中因為會造成其中一個 producer 錯誤,
    // 因此 producer.compute_at(consumer, x) 與 producer.compute_at(consumer, y) 都無法使用
    // 所以必須使用 producer.compute_root() 來 inline    // 這裡想在 conumer 兩個更新步驟的內層 loop 中作 producer 計算
    // Halide 不允許在單一 Func 上套用不同的排程steps.
    // 但能夠藉由 producer 產生兩個 wrapper 並且在 wrapper 上作排程

    // Attempt 2:
    Func producer_1, producer_2, consumer_2;
    producer_1(x, y) = producer(x, y);
    producer_2(x, y) = producer(x, y);

    consumer_2(x, y) = x + y;
    consumer_2(x, 0) += producer_1(x, x);
    consumer_2(0, y) += producer_2(y, 9-y);

    // wrapper 函數提供兩個個別的 producer handle
    // 因此能夠以不同方式來排程
    producer_1.compute_at(consumer_2, x);
    producer_2.compute_at(consumer_2, y);

    Buffer<int> halide_result = consumer_2.realize(10, 10);
    // 下圖為視覺化結果

    // 等義的 C code:
    int c_result[10][10];
    // consumer 純定義的步驟
    for (int y = 0; y < 10; y++) {
        for (int x = 0; x < 10; x++) {
            c_result[y][x] = x + y;
        }
    }
    // 第一個 consumer 更新步驟
    for (int x = 0; x < 10; x++) {
        int producer_1_storage[1];
        producer_1_storage[0] = (x * x) / 10 + 8;
        c_result[0][x] += producer_1_storage[0];
    }
    // 第二個 consumer 更新步驟
    for (int y = 0; y < 10; y++) {
        int producer_2_storage[1];
        producer_2_storage[0] = (y * (9-y)) / 10 + 8;
        c_result[y][0] += producer_2_storage[0];
    }

    // 確認結果一致
    for (int y = 0; y < 10; y++) {
        for (int x = 0; x < 10; x++) {
            if (halide_result(x, y) != c_result[y][x]) {
                printf("halide_result(%d, %d) = %d instead of %d\n",
                       x, y, halide_result(x, y), c_result[y][x]);
                return -1;
            }
        }
    }
}

{
    // Case 5: 在一個 consumer 特定變數的歸納範圍排程下排程 producer
    // 這裡不只是在 consumer 純定義變數上限制 producer 的排程W
    // 若 producer 只用在特定歸納區域(RDom), 依然可以對其做排程
    Func producer, consumer;

    RDom r(0, 5);
    producer(x) = x % 8;
    consumer(x) = x + 10;
    consumer(x) += r + producer(x + r);

    producer.compute_at(consumer, r);

    Buffer<int> halide_result = consumer.realize(10);
    // 下圖為視覺化結果

    

    // 等義的 C CODE:
    int c_result[10];
    // consumer 純定義步驟
    for (int x = 0; x < 10; x++)  {
        c_result[x] = x + 10;
    }
    // consumer 的更新步驟
    for (int x = 0; x < 10; x++) {
        // 歸納區域的 loop 總是內層 loop.
        for (int r = 0; r < 5; r++) {
            // 已經將 producer 的儲存空間與計算做排程於此
            // 這裡只需要單一數值
            int producer_storage[1];
            // producer 純定義步驟
            producer_storage[0] = (x + r) % 8;

            // 在 consumer 更新步驟中使用 producer
            c_result[x] += r + producer_storage[0];
        }
    }

    // 確認結果一致
    for (int x = 0; x < 10; x++) {
        if (halide_result(x) != c_result[x]) {
            printf("halide_result(%d) = %d instead of %d\n",
                   x, halide_result(x), c_result[x]);
            return -1;
        }
    }
}

}

// 在 producer-consumer 中使用歸納的實際的例子
{

// 歸納的預設排程對於類似 convolution 的運算是不錯的.
// 像是下列以 clamp-to-edge 邊界條件於測試影像上計算 5x5 box-blur

// 首先加入邊界條件
Func clamped = BoundaryConditions::repeat_edge(input);

// 定義一個自 (-2, -2) 起始的 5x5 box
RDom r(-2, 5, -2, 5);

// 計算所有 5x5 pixel 的總和
Func local_sum;
local_sum(x, y) = 0; // 以 32-bit 整數計算總和
local_sum(x, y) += clamped(x + r.x, y + r.y);

// 將總和除以 25 來取得平均數
Func blurry;
blurry(x, y) = cast(local_sum(x, y) / 25);

Buffer<uint8_t> halide_result = blurry.realize(input.width(), input.height());

// 預設排程會將 'clamped' 以 inline 置入更新步驟 'local_sump' 中
// 因為 clamped 只有純定義, 因此預設排程為 fully-inlined.
// 接著會計算每個 x 座標的模糊所需的 local_sum,
// 因為歸納的預設排程是最內層計算.
// 以下為等義的 C code:

Buffer<uint8_t> c_result(input.width(), input.height());
for (int y = 0; y < input.height(); y++) {
    for (int x = 0; x < input.width(); x++) {
        int local_sum[1];
        // local_sum 的純定義步驟計算
        local_sum[0] = 0;
        // local_sum 的更新步驟
        for (int r_y = -2; r_y <= 2; r_y++) {
            for (int r_x = -2; r_x <= 2; r_x++) {
                // 邊界座標的 inlin 置於更新步驟中.
                int clamped_x = std::min(std::max(x + r_x, 0), input.width()-1);
                int clamped_y = std::min(std::max(y + r_y, 0), input.height()-1);
                local_sum[0] += input(clamped_x, clamped_y);
            }
        }
        // 模糊的純定義
        c_result(x, y) = (uint8_t)(local_sum[0] / 25);
    }
}

// 確認結果一致
for (int y = 0; y < input.height(); y++) {
    for (int x = 0; x < input.width(); x++) {
        if (halide_result(x, y) != c_result(x, y)) {
            printf("halide_result(%d, %d) = %d instead of %d\n",
                   x, y, halide_result(x, y), c_result(x, y));
            return -1;
        }
    }
}

}

// 歸納輔助 (reduction helpers)
{

// 在 Halide.h 中提供了計算小型歸納與最內層排程的數個歸納輔助.
// 最有用的一個是 "sum"
Func f1;
RDom r(0, 100);
f1(x) = sum(r + x) * 7;

// Sum 建立一個小型匿名 Func 來作歸納. 其等義於:
Func f2;
Func anon;
anon(x) = 0;
anon(x) += r + x;
f2(x) = anon(x) * 7;

// 即便透過 f1 參考一個歸納範圍, 它是個純函數.
// 歸納範圍已經被涵蓋來定義的內部匿名歸納

Buffer<int> halide_result_1 = f1.realize(10);
Buffer<int> halide_result_2 = f2.realize(10);

// 等義的 C code:
int c_result[10];
for (int x = 0; x < 10; x++) {
    int anon[1];
    anon[0] = 0;
    for (int r = 0; r < 100; r++) {
        anon[0] += r + x;
    }
    c_result[x] = anon[0] * 7;
}

// 確認結果是否一致
for (int x = 0; x < 10; x++) {
    if (halide_result_1(x) != c_result[x]) {
        printf("halide_result_1(%d) = %d instead of %d\n",
               x, halide_result_1(x), c_result[x]);
        return -1;
    }
    if (halide_result_2(x) != c_result[x]) {
        printf("halide_result_2(%d) = %d instead of %d\n",
               x, halide_result_2(x), c_result[x]);
        return -1;
    }
}

}

// 使用歸納輔助的複雜範例
{

// 其他歸納輔助包含了 "product", "minimum", "maximum", "argmin" 與 "argmax"
// 使用 argmin 與 argmax 需要了解下一節所介紹的 tuples
// 這裡使用 minimum 與 maximum 來計算灰階影像的局部擴散

// 首先加入輸入的邊界條件
Func clamped;
Expr x_clamped = clamp(x, 0, input.width()-1);
Expr y_clamped = clamp(y, 0, input.height()-1);
clamped(x, y) = input(x_clamped, y_clamped);

RDom box(-2, 5, -2, 5);
// 計算局部最大值減去局部最小值
Func spread;
spread(x, y) = (maximum(clamped(x + box.x, y + box.y)) -
                minimum(clamped(x + box.x, y + box.y)));

// 以 32 scanline 做分割來計算
Var yo, yi;
spread.split(y, yo, yi, 32).parallel(yo);

// 在每個 32 scanline 分割中沿著 x 方向作 vectorization.
// 相關不明確的 vectorization 是在 spread 中沿著 x 方向作計算 
// 其中包含了亦已 vectorize 的 minimum 與 maximum 輔助函數.
spread.vectorize(x, 16);

// 當需要在 circular buffer 的數值時, 藉由 padding 來套用邊界條件
clamped.store_at(spread, yo).compute_at(spread, yi);

Buffer<uint8_t> halide_result = spread.realize(input.width(), input.height());

// 等義的 C code 幾乎太過可怕而難以思考 (花了作者很長時間除錯)
// 這次為了同時測量 Halide 與 C 版本,
// 將使用 SSE intrinsics 來做 vectorization
// 以及 openmp 來平行化 loop (需要使用 -fopenmp 來編譯, 以取得正確的時間)
#ifdef __SSE2__

// 不用將配置輸出 buffer 的時間計入
Buffer<uint8_t> c_result(input.width(), input.height());

#ifdef _OPENMP
double t1 = current_time();
#endif

// 執行 100 次來計算平均時間
for (int iters = 0; iters < 100; iters++) {

    #pragma omp parallel for
    for (int yo = 0; yo < (input.height() + 31)/32; yo++) {
        int y_base = std::min(yo * 32, input.height() - 32);

        // 計算被限制在大小為 8 的 circular buffer
        // (大於 5 的最小 2 的冪次方). Each thread
        // 每個執行緒需要自己的配置, 因此在此處理

        int clamped_width = input.width() + 4;
        uint8_t *clamped_storage = (uint8_t *)malloc(clamped_width * 8);

        for (int yi = 0; yi < 32; yi++) {
            int y = y_base + yi;

            uint8_t *output_row = &c_result(0, y);

            // 計算此 scanline 的邊界限制, 並且跳過此 slice 已計算的列
            int min_y_clamped = (yi == 0) ? (y - 2) : (y + 2);
            int max_y_clamped = (y + 2);
            for (int cy = min_y_clamped; cy <= max_y_clamped; cy++) {
                // 使用 bitmasking 方式計算要填入 circular buffer 的哪一列:
                uint8_t *clamped_row =
                    clamped_storage + (cy & 7) * clamped_width;

                // 藉由對 y 座標 clamp 來計算要讀取輸入的哪一列:
                int clamped_y = std::min(std::max(cy, 0), input.height()-1);
                uint8_t *input_row = &input(0, clamped_y);

                // 搭配 padding 填入
                for (int x = -2; x < input.width() + 2; x++) {
                    int clamped_x = std::min(std::max(x, 0), input.width()-1);
                    *clamped_row++ = input_row[clamped_x];
                }
            }

            // 對輸出的純定義步驟沿著 x 方向的 vector 做計算
            for (int x_vec = 0; x_vec < (input.width() + 15)/16; x_vec++) {
                int x_base = std::min(x_vec * 16, input.width() - 16);

                // 配置 minimum 與 maximum 輔助函數所需的儲存空間
                // 單一 vector 大小已足夠.
                __m128i minimum_storage, maximum_storage;

                // The pure step for the maximum is a vector of zeros
                // maximum 的純定義步驟是數值為 0 的 vector
                maximum_storage = _mm_setzero_si128();

                // maximum 的更新步驟
                for (int max_y = y - 2; max_y <= y + 2; max_y++) {
                    uint8_t *clamped_row =
                        clamped_storage + (max_y & 7) * clamped_width;
                    for (int max_x = x_base - 2; max_x <= x_base + 2; max_x++) {
                        __m128i v = _mm_loadu_si128(
                            (__m128i const *)(clamped_row + max_x + 2));
                        maximum_storage = _mm_max_epu8(maximum_storage, v);
                    }
                }

                // minimum 的純定義步驟是數值為 1 的 vector .
                // 藉由與自身相較來建立
                minimum_storage = _mm_cmpeq_epi32(_mm_setzero_si128(),
                                                  _mm_setzero_si128());

                // minimum 的更新步驟
                for (int min_y = y - 2; min_y <= y + 2; min_y++) {
                    uint8_t *clamped_row =
                        clamped_storage + (min_y & 7) * clamped_width;
                    for (int min_x = x_base - 2; min_x <= x_base + 2; min_x++) {
                        __m128i v = _mm_loadu_si128(
                            (__m128i const *)(clamped_row + min_x + 2));
                        minimum_storage = _mm_min_epu8(minimum_storage, v);
                    }
                }

                // 計算擴散
                __m128i spread = _mm_sub_epi8(maximum_storage, minimum_storage);

                // 寫入
                _mm_storeu_si128((__m128i *)(output_row + x_base), spread);

            }
        }

        free(clamped_storage);
    }
}

// Skip the timing comparison if we don't have openmp
// enabled. Otherwise it's unfair to C.
// 若無使用 OpenMP 則忽略比較時間, 否則對 C 不公平
#ifdef _OPENMP
double t2 = current_time();

// 執行 Halide 版本並避免 JIT 編譯的 overhead.
// 此外也執行 100 次
for (int iters = 0; iters < 100; iters++) {
    spread.realize(halide_result);
}

double t3 = current_time();

// 印出時間的結果, 在作者的機器兩者對於 4-megapixel 的輸入都需要大約 3ms
// 這是合理的, 因為使用相同的 vectorization 與 parallelization 策略
// 然而 Halide 叫容易閱讀, 撰寫, 除錯, 修改與移植
printf("Halide spread took %f ms. C equivalent took %f ms\n",
       (t3 - t2)/100, (t2 - t1)/100);

#endif // _OPENMP

// 確認結果一致:
for (int y = 0; y < input.height(); y++) {
    for (int x = 0; x < input.width(); x++) {
        if (halide_result(x, y) != c_result(x, y)) {
            printf("halide_result(%d, %d) = %d instead of %d\n",
                   x, y, halide_result(x, y), c_result(x, y));
            return -1;
        }
    }
}

#endif // __SSE2__

}