# 2016q3 Homework1 (raytracing)
contributed by <`janetwei`>
reviewed by <`ierosodin`>
- 可以利用gnuplot畫出優化前優化後的差異
- 利用pthread實作raytracing的平行化
- SIMD來加速整個資料的處理
- 程式碼和結果不要用圖片表達
## 開發環境
先載入ㄧ些開發工具

$ sudo apt-get update
$ sudo apt-get install graphviz
$ sudo apt-get install imagemagick

## Raytracing
光線追蹤 (Ray tracing) 是三維電腦圖形學中的描繪演算法 (rendering algorithm)，跟蹤從眼睛發出的光線，而非是光源發出的光線，因此與物理現象有些差距


convert ppm to png
check圖片輸出是否符合預期

以 `make PROFILE=1` 重新編譯程式碼，並且學習 gprof

可以發現某些function被呼叫了千萬次
可以試著從這方面去優化程式
## graphviz
因為graphviz需要使用特定格式,所以需要安裝gprof2dot轉換格式
但因為需要用到pip,所以先安裝pyphon

$sudo apt install python python-pip
$sudo pip install --upgrade pip
$sudo pip install gprof2dot

將gprof所得出的結果用graphviz表示

這是用編譯器最佳化選項得出的結果
時間直接變成0.759s
因為將-O0- 變成 -Ofast-

但是我們的目標則是手動優化程式的結果比電腦更好（目標好遠大阿QQ）

## dot_product優化(loop unrolling)

由gprof分析 dot_product被呼叫了將近七千萬次，因此第一個想從這裡下手

double dot_product(const double *v1, const double *v2)
{
    double dp = 0.0;
    for (int i = 0; i < 3; i++)
        dp += v1[i] * v2[i];
    return dp;
}


將迴圈展開

double dot_product(const double *v1, const double *v2)
{
    double dp = 0.0;
    for (int i = 0; i < 3; i++)
        dp = v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2];
    return dp;
}

在這個做完之後，我才發現這個叫做loop unrolling

發現已經減少到6.62s了，所以將其他function也用loop unrolling

>> 避免用圖片來表達程式碼和輸出的文字訊息，日後才能方便搜尋和複製，敬請配合 [name=jserv]

結果已經減少到5.51s了

## static inline
static只對定義(有寫出程式碼)有用
inline 對編譯器有提示的作用,將function called轉成以下的程式碼,利用inline將程式碼展開,少了程式呼叫的成本

但是我們在程式一開始時設定為-O0,最佳化被關掉
所以inline沒有作用,因此我們要強制inline作用(force inline)
將`$ static inline`變成`$ static__attribute__((always_inline))`

已經將時間減少為2.638s

Each sample counts as 0.01 seconds.
  %   cumulative   self              self     total           
 time   seconds   seconds    calls   s/call   s/call  name    
 30.18      0.54     0.54 13861875     0.00     0.00  raySphereIntersection
 29.62      1.07     0.53 13861875     0.00     0.00  rayRectangularIntersection
 11.74      1.28     0.21  2110576     0.00     0.00  compute_specular_diffuse
  7.26      1.41     0.13  1048576     0.00     0.00  ray_color
  6.15      1.52     0.11  2110576     0.00     0.00  localColor
  4.47      1.60     0.08  1048576     0.00     0.00  rayConstruction
  4.19      1.68     0.08  4620625     0.00     0.00  ray_hit_object
  2.24      1.72     0.04  1241598     0.00     0.00  reflection
  1.68      1.75     0.03        1     0.03     1.79  raytracing
  1.12      1.77     0.02  1241598     0.00     0.00  protect_color_overflow
  0.56      1.78     0.01  1241598     0.00     0.00  refraction
  0.56      1.79     0.01  1204003     0.00     0.00  idx_stack_push
  0.28      1.79     0.01   113297     0.00     0.00  fresnel
  0.00      1.79     0.00  2558386     0.00     0.00  idx_stack_empty
  0.00      1.79     0.00  2520791     0.00     0.00  idx_stack_top
  0.00      1.79     0.00  1048576     0.00     0.00  idx_stack_init
  0.00      1.79     0.00    37595     0.00     0.00  idx_stack_pop
  0.00      1.79     0.00        3     0.00     0.00  append_rectangular
  0.00      1.79     0.00        3     0.00     0.00  append_sphere
  0.00      1.79     0.00        2     0.00     0.00  append_light
  0.00      1.79     0.00        1     0.00     0.00  calculateBasisVectors
  0.00      1.79     0.00        1     0.00     0.00  delete_light_list
  0.00      1.79     0.00        1     0.00     0.00  delete_rectangular_list

可以知道接下來從raySphereIntersection著手優化

## OpenMD
此網站上有著簡單介紹平行化處理的方法
https://kheresy.wordpress.com/2006/06/09/
對於一般單一執行緒（single thread）的程式，多核心的處理器並沒有辦法提升它的處理效能；不過對於多執行緒（multi thread）的程式，就可以透過不同的核心同時計算，來達到加速的目的了！

如果要在for迴圈平行化處理,在前面加上一行`$ #pragma omp parallel for`即可

如果想利用 OpenMP 把 mian() 裡面的迴圈平行化處理呢？只需要修改成下面的樣子：

#include <omp.h>
#include <stdio.h>
#include <stdlib.h>

void Test( int n )
{
 for( int i = 0; i < 10000; ++ i )
 {
  //do nothing, just waste time
 }
 printf( "%d, ", n );
}

int main(int argc, char* argv[])
{
 #pragma omp parallel for
 for( int i = 0; i < 10; ++ i )
  Test( i );
 system( "pause" );
}

只要兩個或以上個for迴圈裏面的變數不是互相獨立,而且不是在平行化處理裡面宣告的變數
要設定成private:在每個thread裡面都有個變數的副本

在function raytracing中在3個for迴圈前加上對3個變數private
以及發現下面也可以做Loop unrolling

void raytracing(uint8_t *pixels, color background_color,
                rectangular_node rectangulars, sphere_node spheres,
                light_node lights, const viewpoint *view,
                int width, int height)
{
    point3 u, v, w, d;
    color object_color = { 0.0, 0.0, 0.0 };

    /* calculate u, v, w */
    calculateBasisVectors(u, v, w, view);

    idx_stack stk;
    

    int factor = sqrt(SAMPLES);
    #pragma omp parallel for num_threads(512)\
    private(stk), private(d) \
    private(object_color)
    for (int j = 0; j < height; j++) {
        for (int i = 0; i < width; i++) {
            double r = 0, g = 0, b = 0;
            /* MSAA */
          //  for (int s = 0; s < SAMPLES; s++) {
                idx_stack_init(&stk);
                rayConstruction(d, u, v, w,
                                i * factor + 0 / factor,
                                j * factor + 0 % factor,view,
                                width * factor, height * factor);
                if (ray_color(view->vrp, 0.0, d, &stk, rectangulars, spheres,
                              lights, object_color,
                              MAX_REFLECTION_BOUNCES)) {
                    r += object_color[0];
                    g += object_color[1];
                    b += object_color[2];
                } else {
                    r += background_color[0];
                    g += background_color[1];
                    b += background_color[2];
                }
                pixels[((i + (j * width)) * 3) + 0] = r * 255 / SAMPLES;
                pixels[((i + (j * width)) * 3) + 1] = g * 255 / SAMPLES;
                pixels[((i + (j * width)) * 3) + 2] = b * 255 / SAMPLES;
                
                idx_stack_init(&stk);
                rayConstruction(d, u, v, w,
                                i * factor + 1 / factor,
                                j * factor + 1 % factor,
                                view,
                                width * factor, height * factor);
                if (ray_color(view->vrp, 0.0, d, &stk, rectangulars, spheres,
                              lights, object_color,
                              MAX_REFLECTION_BOUNCES)) {
                    r += object_color[0];
                    g += object_color[1];
                    b += object_color[2];
                } else {
                    r += background_color[0];
                    g += background_color[1];
                    b += background_color[2];
                }
                pixels[((i + (j * width)) * 3) + 0] = r * 255 / SAMPLES;
                pixels[((i + (j * width)) * 3) + 1] = g * 255 / SAMPLES;
                pixels[((i + (j * width)) * 3) + 2] = b * 255 / SAMPLES;

                idx_stack_init(&stk);
                rayConstruction(d, u, v, w,
                                i * factor + 2 / factor,
                                j * factor + 2 % factor,
                                view,
                                width * factor, height * factor);
                if (ray_color(view->vrp, 0.0, d, &stk, rectangulars, spheres,
                              lights, object_color,
                              MAX_REFLECTION_BOUNCES)) {
                    r += object_color[0];
                    g += object_color[1];
                    b += object_color[2];
                } else {
                    r += background_color[0];
                    g += background_color[1];
                    b += background_color[2];
                }
                pixels[((i + (j * width)) * 3) + 0] = r * 255 / SAMPLES;
                pixels[((i + (j * width)) * 3) + 1] = g * 255 / SAMPLES;
                pixels[((i + (j * width)) * 3) + 2] = b * 255 / SAMPLES;

                idx_stack_init(&stk);
                rayConstruction(d, u, v, w,
                                i * factor + 3 / factor,
                                j * factor + 3 % factor,
                 view,
                                width * factor, height * factor);
                if (ray_color(view->vrp, 0.0, d, &stk, rectangulars, spheres,
                              lights, object_color,
                              MAX_REFLECTION_BOUNCES)) {
                    r += object_color[0];
                    g += object_color[1];
                    b += object_color[2];
                } else {
                    r += background_color[0];
                    g += background_color[1];
                    b += background_color[2];
                }
                pixels[((i + (j * width)) * 3) + 0] = r * 255 / SAMPLES;
                pixels[((i + (j * width)) * 3) + 1] = g * 255 / SAMPLES;
                pixels[((i + (j * width)) * 3) + 2] = b * 255 / SAMPLES;
            //}
        }
    }
}

結果變成1.327s了

janet@janet-MacBookAir:~/raytracing$ ./raytracing 
# Rendering scene
Done!
Execution time of raytracing() : 1.327288 sec

因為執行`make PROFILE=1`會比直接`make`多花時間
因此就試試看直接`make`

janet@janet-MacBookAir:~/raytracing$ ./raytracing 
# Rendering scene
Done!
Execution time of raytracing() : 0.788209 sec

之後得到的時間為0.788s

## SIMD
可以同一時間點處理好幾筆資料