Raytracing - HackMD

# Raytracing contributed by <`zhanyangch`>, <`hunng`>, <`yanang`> ## 預期目標 * Task: 比照 [B02: raytracing](https://hackmd.io/s/HyuBWDwYl) 要求，以 SIMD 改善 math-toolkit.h 定義的若干數學操作函式，進而推廣到整個光線追蹤的演算法 (可在不影響正確性的前提下，將原本的實作替換)。 * Task: 參考 Joshua Barczak 的 [SIMD Raytracing 實作](http://www.joshbarczak.com/blog/?p=787)，改寫原有的程式碼，但是 model 應該和 [B02: raytracing](https://hackmd.io/s/HyuBWDwYl) 一致 ## math-toolkit.h with SSE Streaming SIMD Extensions：8個 128 位元暫存器，每個暫存器可以容納 2 個 64位元雙精度浮點數 Intrinsics：是一種外表類似一般的函數，但是實際上會被 compiler 直接編譯成組合語言的東西。 * 常用指令 ```clike __m128d _mm_loadu_pd (double const* mem_addr) __m128d _mm_load_sd (double const* mem_addr)//高位為0 __m128d _mm_add_pd (__m128d a, __m128d b) __m128d _mm_div_pd (__m128d a, __m128d b) __m128d _mm_mul_pd (__m128d a, __m128d b) void _mm_storeu_pd (double* mem_addr, __m128d a) void _mm_store_pd1 (double* mem_addr, __m128d a)//皆填入 a[63:0] __m128d _mm_sqrt_pd(__m128d a) _mm_shuffle_pd (__m128d a, __m128d b, int imm8) __m128d _mm_unpackhi_pd (__m128d a, __m128d b) double _mm_cvtsd_f64(__m128d a) ``` * math-toolkit-sse.h 實做 ```clike __attribute__((always_inline)) static inline void normalize(double *v) { __m128d I0 = _mm_loadu_pd(v);//v0 v1 __m128d I1 = _mm_loadu_pd(v+1);//v1 v2 __m128d I2 = _mm_mul_pd (I0,I0);//v0^2 v1^2 __m128d I3 = _mm_mul_pd (I1,I1);//v1^2 v2^2 __m128d I4 = _mm_add_pd (I2,I3);//v0^2+v1^2 v1^2+v2^2 __m128d T0 = _mm_unpackhi_pd(I3,I3); // v2^2 v2^2 __m128d I5 = _mm_add_pd(I4,T0);//v0^2+v1^2+v2^2 __m128d I6 = _mm_unpacklo_pd(I5,I5); __m128d I7 =_mm_sqrt_pd(I6); I0=_mm_div_pd(I0,I7); I1=_mm_div_pd(I1,I7); _mm_storeu_pd(v,I0); _mm_storeu_pd(v+1,I1); } __attribute__((always_inline)) static inline double length(const double *v) { __m128d I0 = _mm_loadu_pd(v);//v0 v1 __m128d I1 = _mm_loadu_pd(v+1);//v1 v2 __m128d I2 = _mm_mul_pd (I0,I0);//v0^2 v1^2 __m128d I3 = _mm_mul_pd (I1,I1);//v1^2 v2^2 __m128d I4 = _mm_add_pd (I2,I3);//v0^2+v1^2 v1^2+v2^2 __m128d T0 = _mm_unpackhi_pd(I3,I3); // v2^2 v2^2 __m128d I5 = _mm_add_pd(I4,T0);//v0^2+v1^2+v2^2 __m128d I6 = _mm_sqrt_pd(I5); return _mm_cvtsd_f64(I6); } __attribute__((always_inline)) static inline void add_vector(const double *a, const double *b, double *out) { __m128d I0 = _mm_loadu_pd(a);//a0 a1 __m128d I1 = _mm_loadu_pd(a+1);//a1 a2 __m128d I2 = _mm_loadu_pd(b);//b0 b1 __m128d I3 = _mm_loadu_pd(b+1);//b1 b2 __m128d T0 = _mm_add_pd(I0,I2); __m128d T1 = _mm_add_pd(I1,I3); _mm_storeu_pd(out,T0); _mm_storeu_pd(out+1,T1); } __attribute__((always_inline)) static inline void subtract_vector(const double *a, const double *b, double *out) { __m128d I0 = _mm_loadu_pd(a);//a0 a1 __m128d I1 = _mm_loadu_pd(a+1);//a1 a2 __m128d I2 = _mm_loadu_pd(b);//b0 b1 __m128d I3 = _mm_loadu_pd(b+1);//b1 b2 __m128d T0 = _mm_sub_pd(I0,I2); __m128d T1 = _mm_sub_pd(I1,I3); _mm_storeu_pd(out,T0); _mm_storeu_pd(out+1,T1); } __attribute__((always_inline)) static inline void multiply_vectors(const double *a, const double *b, double *out) { __m128d I0 = _mm_loadu_pd(a);//a0 a1 __m128d I1 = _mm_loadu_pd(a+1);//a1 a2 __m128d I2 = _mm_loadu_pd(b);//b0 b1 __m128d I3 = _mm_loadu_pd(b+1);//b1 b2 __m128d T0 = _mm_mul_pd(I0,I2); __m128d T1 = _mm_mul_pd(I1,I3); _mm_storeu_pd(out,T0); _mm_storeu_pd(out+1,T1); } static inline void multiply_vector(const double *a, double b, double *out) { __m128d I0 = _mm_loadu_pd(a);//a0 a1 __m128d I1 = _mm_loadu_pd(a+1);//a1 a2 __m128d I2 = _mm_load_pd1(&b);//b0 b1 __m128d T0 = _mm_mul_pd(I0,I2); __m128d T1 = _mm_mul_pd(I1,I2); _mm_storeu_pd(out,T0); _mm_storeu_pd(out+1,T1); } __attribute__((always_inline)) static inline void cross_product(const double *v1, const double *v2, double *out) { __m128d I0 = _mm_loadu_pd(v1);//v1[0] v1[1] __m128d I1 = _mm_loadu_pd(v1+1);//v1[1] v1[2] __m128d I2 = _mm_loadu_pd(v2);//v2[0] v2[1] __m128d I3 = _mm_loadu_pd(v2+1);//v2[1] v2[2] __m128d I4 = _mm_shuffle_pd(I1,I1,1);//v1[2] __m128d I5 = _mm_shuffle_pd(I3,I3,1);//v2[2] __m128d T1 = _mm_mul_pd(I0,I3); __m128d T2 = _mm_mul_pd(I1,I2); __m128d T3 = _mm_sub_pd(T1,T2);// out[2] out[0] __m128d T4 = _mm_mul_pd(I4,I2); __m128d T5 = _mm_mul_pd(I5,I0); __m128d T6 = _mm_sub_pd(T4,T5);// out[1] __m128d T7 =_mm_shuffle_pd(T3,T6,1);//out[0] out[1] _mm_storeu_pd(out,T7); *(out+2)=_mm_cvtsd_f64(T3); } __attribute__((always_inline)) static inline double dot_product(const double *v1, const double *v2) { __m128d I0 = _mm_loadu_pd(v1);//v1[0] v1[1] __m128d I1 = _mm_loadu_pd(v1+1);//v1[1] v1[2] __m128d I2 = _mm_loadu_pd(v2);//v2[0] v2[1] __m128d I3 = _mm_loadu_pd(v2+1);//v2[1] v2[2] __m128d T1 = _mm_mul_pd(I0,I2); __m128d T2 = _mm_mul_pd(I1,I3); __m128d T3 = _mm_add_pd(T1,T2); __m128d T4 = _mm_unpackhi_pd(T2,T2); __m128d T5 = _mm_add_pd(T3,T4); return _mm_cvtsd_f64(T5); } __attribute__((always_inline)) static inline void scalar_triple_product(const double *u, const double *v, const double *w, double *out) { cross_product(v, w, out); multiply_vectors(u, out, out); } __attribute__((always_inline)) static inline double scalar_triple(const double *u, const double *v, const double *w) { double tmp[3]; cross_product(w, u, tmp); return dot_product(v, tmp); } ``` ## 執行測試 * 正確性 ```shell $ ./raytracing-sse # Rendering scene Done! Execution time of raytracing() : 2.419169 sec $ diff baseline.ppm out.ppm ``` * 時間比較在此的 orig 為 loop unrolling 的版本 1. 不使用 openmp ``` orig:2.247861 sec sse:5.866375 sec ``` 2. 使用 openmp num_threads (2) ``` orig:1.289200 sec sse:2.836198 sec ``` * 分析使用 SSE 的版本反而比原版慢，推測可能的原因有 1. SSE 的之資料長度為 128bit 運算 3個 double 需要切個成兩部份計算，造成指令數目增加 2. 載入跟儲存的部份佔全部指令的比例太高，使用　_mm_load _mm_store ，取代 _mm_loadu _mm_storeu 應可以加快載入跟儲存的速度，但記憶體需對齊 16byte ，目前試著加入　__attribute__((aligned(16))但未能成功 * 修改 raytracing.c 為了解決載入跟儲存的部份佔全部指令的比例太高的問題，我們決定採用幾個合併多個功能的函數 ``` clike void multiply_vector_add(const double *a, double b,const double *c,double *out) //b*a+c void subtract_vector_normalize(const double *a, const double *b, double *out) void multiply_vector_normalize(const double *a, double b, double *out) void add_vectors4(const double *a, const double *b, const double *c,const double *d,double *out)//a+b+c+d void cross_product_normalize(const double *v1, const double *v2, double *out) ``` 1. 不使用 openmp ``` sse　with additional functions: 4.545110 sec ``` 2. 使用 openmp num_threads (2) ``` sse with additional functions: 2.658897 sec ``` 不使用 openmp 的相對改善幅度較大，但皆無法超過原始版本。 ## math-toolkit.h with AVX **針對 math-toolkit.h 進行 avx 的改寫** 起先因為是計算的向量皆為三項的依照叫淺顯的寫法就是直接使用 `_mm256_set_pd` 以下方的 subtrack_vector 為例 ```c __attribute__((always_inline)) static inline void subtract_vector(const double *a, const double *b, double *out) { register __m256d ymm0, ymm1; double temp[4] = {0, 0, 0, 0}; ymm0 = _mm256_set_pd (*temp, *(a+2), *(a+1), *a); // ymm0 = {a0, a1, a2, 0} ymm1 = _mm256_set_pd (*temp, *(b+2), *(b+1), *b); // ymm1 = {b0, b1, b2, 0} ymm0 = _mm256_sub_pd (ymm0, ymm1); // ymm0 = {a0-b0, a1-b1, a2-b2, 0} _mm256_storeu_pd (temp, ymm0); *out = *temp; *(out+1) = *(temp+1); *(out+2) = *(temp+2); } ``` 但因為 `_mm256_set_pd` 的 latency 較高所以我們修改 `primitives.h` 中 point3, point4 的存儲數 3 -> 4 ```c typedef double point3[4]; typedef double point4[4]; ``` 改用 `_mm256_loadu_pd` 作為載入資料手段 ```c __attribute__((always_inline)) static inline void subtract_vector_avx(const double *a, const double *b, double *out) { register __m256d ymm0, ymm1; ymm0 = _mm256_loadu_pd (a); // ymm0 = {a0, a1, a2, X} ymm1 = _mm256_loadu_pd (b); // ymm1 = {b0, b1, b2, X} ymm0 = _mm256_sub_pd (ymm0, ymm1); // ymm0 = {a0+b0, a1+b1, a2+b2, X} _mm256_storeu_pd (out, ymm0); // temp = {a0+b0, a1+b1, a2+b2, X} } ``` 則可以讓整體程式碼較為簡潔但缺點就是需要注意其最後一項的值是否會影響其結果 (像是後面用到的 `hadd`) 為了解決這個問題，我們改用了 `_mm256_maskload_pd` ```c __attribute__((always_inline)) static inline void subtract_vector(const double *a, const double *b, double *out) { register __m256d ymm0, ymm1; ymm0 = _mm256_maskload_pd (a, mask); // ymm0 = {a0, a1, a2, 0} ymm1 = _mm256_maskload_pd (b, mask); // ymm0 = {b0, b1, b2, 0} ymm0 = _mm256_sub_pd (ymm0, ymm1); // ymm0 = {a0-b0, a1-b1, a2-b2, 0} _mm256_maskstore_pd (out, mask, ymm0); } ``` 若要使用 mask 則需要指定一個全域變數 `mask` 並在使用raytracing時叫用 `init_mask()` ```c void init_mask() { long int value = 0xAAAAAAAAAAAAAAAA; mask = _mm256_set_epi64x (0, value, value, value); } ``` 但這又衍生出了 latency 的問題若每個 load, store 都使用 maskload, maskstore 會降低整體的速度要衡量每個 load, store 是否需要使用 mask 而在 dot_product 中會使用到 `_mm256_hadd_pd` ```c __attribute__((always_inline)) static inline double dot_product(const double *a, const double *b) { double temp[4] = {0, 0, 0, 0}; register __m256d ymm0, ymm1, ymm2; ymm0 = _mm256_maskload_pd (a, mask); // ymm0 = {a0, a1, a2, 0} ymm1 = _mm256_maskload_pd (b, mask); // ymm0 = {b0, b1, b2, 0} ymm2 = _mm256_setzero_pd(); // ymm2 = {0, 0, 0, 0} ymm0 = _mm256_mul_pd (ymm0, ymm1); // ymm0 = {a0*b0, a1*b1, a2*b2, 0} ymm0 = _mm256_hadd_pd (ymm0, ymm2); // ymm0 = {a0*b0+a1*b1, 0, a2*b2, 0} ymm1 = _mm256_permute2f128_pd (ymm0, ymm0, 131); // ymm0 = { a2*b2, 0, 0, 0} ymm0 = _mm256_add_pd (ymm0, ymm1); // ymm0 = {a0*b0+a1*b1+a2*b2, 0, 0, 0} _mm256_storeu_pd (temp, ymm0); return *temp; } ``` 原本是使用以下寫法，利用一個 128 bit 的暫存，然後在衍生至 256 bit ```c const __m128d valupper = _mm256_extractf128_pd(ymm0, 1); // valupper = {a2*b2, 0} ymm1 = _mm256_castpd128_pd256 (valupper); // ymm2 = {a2*b2, 0, X, X} X means undefined ``` 但考量到整體可讀性及正確性，改用了 `_mm256_permute2f128_pd` ```c ymm2 = _mm256_permute2f128_pd (ymm3, ymm3, 131); // ymm0 = { v2^2, 0, 0, 0} ``` 而 cross_product 因為使用了叫複雜的內容所以所有值都是利用設定的並為了讓整體運算結果會剛好符合順序載入的順序就沒有那麼直觀 ```c __attribute__((always_inline)) static inline void cross_product(const double *a, const double *b, double *out) { register __m256d ymm0, ymm1, ymm2, ymm3; double temp = 0; ymm0 = _mm256_set_pd (*(a+1), *a, *(a+2), *(a+1)); // ymm0 = {a1, a2, a2, a0} ymm1 = _mm256_set_pd (*b, *(b+1), *(b+1), *(b+2)); // ymm1 = {b2, b1, b0, b2} ymm2 = _mm256_set_pd (temp, temp, *a, *(a+2)); // ymm2 = {a0, a1, 0, 0} ymm3 = _mm256_set_pd (temp, temp, *(b+2), *b); // ymm3 = {b1, b0, 0, 0} ymm0 = _mm256_mul_pd (ymm0, ymm1); // ymm0 = {a1*b2, a2*b1, a2*b0, a0*b2} ymm2 = _mm256_mul_pd (ymm2, ymm3); // ymm2 = {a0*b1, a1*b0, 0, 0} ymm3 = _mm256_hsub_pd (ymm0, ymm2); _mm256_maskstore_pd (out, mask, ymm3); } ``` 在整體完成後因為在實際在 raytracing 的程式中看到了速度大幅的降低所以我們對舊版與 avx 版本的 math-toolkit.h 進行測試 ![](https://i.imgur.com/bzfzvPp.png) 發現只有 normalize 與 length 有較佳的表現 ## Raytracing.c * all function - [x] raySphereIntersection - [x] rayRectangularIntersection - [ ] localColor - [x] compute_specular_diffuse - [x] reflection - [x] refraction - [x] fresnel - [ ] ray_hit_object - [x] rayConstruction - [x] calculateBasisVectors - [ ] protect_color_overflow - [ ] ray_color ## 修改內容 [Joshua Barczak的文章筆記](https://hackmd.io/s/rJhfg1IlZ) 為了使一次能運算 4 個 double 我們首先修改了 point3 的使用改為使用 point3v 其內容為 ```c typedef struct { __m256d x; __m256d y; __m256d z; }point3v; ``` 一個 point3v 將會對應儲存四個需要運算的點 x 第 0 個值, y 第 0 個值, z 第 0 個值為第 0 個點 x 第 1 個值, y 第 1 個值, z 第 1 個值為第 1 個點 x 第 2 個值, y 第 2 個值, z 第 2 個值為第 2 個點 x 第 3 個值, y 第 3 個值, z 第 3 個值為第 3 個點接下來修改了 math-toolkit.h 使其的運算參數從 point3 -> point3v 一次處理四個點的運算舉 cross_product 為例 ```c __attribute__((always_inline)) static inline void cross_product(const point3v *a, const point3v *b, point3v *out) { __m256d x1 = __mm256_mul_pd(a->y, b->z); __m256d x2 = __mm256_mul_pd(a->z, b->y); __m256d y1 = __mm256_mul_pd(a->z, b->x); __m256d y2 = __mm256_mul_pd(a->x, b->z); __m256d z1 = __mm256_mul_pd(a->x, b->y); __m256d z2 = __mm256_mul_pd(a->y, b->x); out->x = __mm256_sub_pd(x1, x2); out->y = __mm256_sub_pd(y1, y2); out->z = __mm256_sub_pd(z1, z2); } ``` 再來就是 Raytracing.c 的改寫 Raytracing.c 內許多函式的寫法為利用 if 判斷後 return 但因為現在為計算四個點，所以當判斷後並不能直接 return 需要考量其他三個點的情況而我們目前的作法是先全部運算後在利用 mask 去還原舊的值（該 return 的點）以 refraction 為例，以下為原檔案 ```c= static void refraction(point3 t, const point3 I, const point3 N, double n1, double n2) { double eta = n1 / n2; double dot_NI = dot_product(N,I); double k = 1.0 - eta * eta * (1.0 - dot_NI * dot_NI); if (k < 0.0 || n2 <= 0.0) t[0] = t[1] = t[2] = 0.0; else { point3 tmp; multiply_vector(I, eta, t); multiply_vector(N, eta * dot_NI + sqrt(k), tmp); subtract_vector(t, tmp, t); } } ``` 修改過後 ```c= static void refraction (point3v *t, const point3v *I, const point3v *N, double n1, double n2) { __m256d n2v = _mm256_set1_pd(n2); __m256d eta = _mm256_set1_pd(n1/n2); __m256d dot_NI = dot_product(N, I); __m256d k = _mm256_set_pd(1); __m256d eta2 = _mm256_mul_pd(eta, eta); // eta2 = eta * eta __m256d dot_NI2 = _mm256_mul_pd(dot_NI, dot_NI); // dot_NI2 = dot_NI * dot_NI dot_NI2 = _mm256_sub_pd(k, dot_NI2); // dot_NI2 = 1 - dot_NI * dot_NI eta2 = _mm256_mul_pd(eta2, dot_NI2); // eta2 = eta * eta * ( 1 - dot_NI * dot_NI) k = _mm256_sub_pd(k, eta2); // k = 1 - eta * eta * ( 1 - dot_NI * dot_NI) #define _CMP_LT_OS 0x01 /* Less-than (ordered, signaling) */ #define _CMP_LE_OS 0x02 /* Less-than-or-equal (ordered, signaling) */ __m256d mzero = _mm256_setzero_pd(); __m256i iall1 = _mm256_set1_epi64x(-1); // iall1 = 0xFF...F __m256d dall1 = _mm256_castsi256_pd(iall1); __m256d if1 = _mm256_cmp_pd(k, mzero, _CMP_LT_OS); // k < 0.0 __m256d if2 = _mm256_cmp_pd(n2v, mzero, _CMP_LE_OS); // n2 <= 0.0 __m256d ifa = _mm256_or_pd(k, n2v); // if (k < 0.0 || n2 <= 0.0) __m256d notifa = _mm256_xor_pd(ifa, dall1); point3v tmp; multiply_vector(I, eta, t); __m256d midpar = _mm256_mul_pd(eta, dot_NI); k = _mm256_sqrt_pd(k); // k = sqrt(k) midpar = _mm256_add_pd(midpar, k); multiply_vector(N, midpar, tmp); subtract_vector(t, tmp, t); t->x = _mm256_and_pd(t->x, notifa); t->y = _mm256_and_pd(t->y, notifa); t->z = _mm256_and_pd(t->z, notifa); } ``` * 處理 if 的方法是利用 `_mm256_cmp_pd` 並適時利用每個 bit 皆為 1 的 `__m256d` 和 `_mm256_xor_pd` 將 0->1 1->0 而每個 bit 皆為 1 的 `__m256d` 則透過整數 `__m256i` 設定為 -1 後轉換 (16, 17 行) * 19~22 為處理 `if (k < 0.0 || n2 <= 0.0)` (原程式 7 行) 而因為這裡為 if 成立就設定為0，所以利用 xor 轉換 bit 最後將求出來的值 and 這個產生出來的 mask (notifa) `t->x = _mm256_and_pd(t->x, notifa);` (32 行) * 但這個方法因為每次使用這個函式就必定會運算假如四個點皆不需運算，則會造成嚴重的資源浪費這也是 [Joshua Barczak](http://www.joshbarczak.com/blog/?p=787) 後面提到他使用到 Ray Reordering 的原因 * 而目前遇到的問題是 * 如何處理遞迴呼叫的問題 e.g., ray_color 會呼叫自己 * 且因為整體 struct 尚未設計完全，寫出的函式沒辦法測試，無法確認其正確性 ## 其他 raytrace 實做 - [ ] 探討 [jaybird19/simd-raytracer](https://github.com/jaybird19/simd-raytracer) simd 切割及multi-threaded方式 ### LinearAllocator ```c++ allocator.init(main_arena_size + ALIGNMENT_PADDING, true); ``` 先一次宣告大量的記憶體 ```c++ pixmap = init_itr<real>((real*)allocator.allocate(pixels_heap_size, SSE_ALIGNMENT), pixels_heap_size); ``` 為每一個資料分配位置 ```c++ //Align address if necessary if (alignment != 0 && (reinterpret_cast<size_t>(aligned_address) & (alignment - 1))) { aligned_address = reinterpret_cast<void*>((reinterpret_cast<size_t>(aligned_address) + (alignment - 1)) & ~(alignment - 1)); alignment_padding = static_cast<uchar>((size_t)aligned_address - (size_t)m_top); } m_allocated_memory += size + alignment_padding; ``` * 先看 if 在哪些情況下成立,alignment 為要對齊的大小 ,通常不為 0,(aligned_address) & (alignment - 1)可以視為取 alignment 的餘數, 注意 alignment 為2^n,才能這樣使用,我們帶數字來看,假設 alignment 為　16 ```clike alignment=0x100;//alignment-1=0xFF aligned_address=0xA00; (aligned_address) & (alignment - 1)//=0 不作 align aligned_address=0xA0A; (aligned_address) & (alignment - 1)//=0x0A 需要 align ``` * 如果需要 align 則找到下一個對齊的位置，並將設置 padding 使下一次 m_allocated_memory 數值正確 ### SIMD trace_section 為 single thread時主要的函式,其針對每4個 pixel(x,y) 產生對應的 rdir 結構，並且呼叫trace_ray 進行計算。 ```clike static void VECTOR_CC trace_section(sse_ray rorigin, uint32 width, uint32 height, int ystart) { sse_ray rdir; sse_pixel pixel; uint32 disp = width*ystart*3; for (uint32 y = ystart; y < height; ++y) { for (uint32 x = 0; x < width - 3; x += 4) { rdir.x = _mm_load_ps(&rtm.cam_ray.x[SsePos(y*width + x)]); rdir.y = _mm_load_ps(&rtm.cam_ray.y[SsePos(y*width + x)]); rdir.z = _mm_load_ps(&rtm.cam_ray.z[SsePos(y*width + x)]); pixel = trace_ray(rorigin, rdir, 0); _mm_store_ps(&rtm.pixmap[SsePos(disp)], pixel.r); disp += 4; _mm_store_ps(&rtm.pixmap[SsePos(disp)], pixel.g); disp += 4; _mm_store_ps(&rtm.pixmap[SsePos(disp)], pixel.b); disp += 4; } } } ``` * trace_ray 中利用 intersect_ray_spheres(rorigin, rdir); 先確定有哪些 spheres 需要被計算(此函式可得到4個點分別需要被計算的 spheres 之 index)在進行計算 ```C++ static ray_sphere_interesection VECTOR_CC intersect_ray_spheres(sse_ray rorigin, sse_ray rdir) { ray_sphere_interesection rs; rs.tpoint = _mm_set_ps1(FLT_MAX); rs.index = _mm_set_ps1(-1.0f); sse_ray rlength; __m128 no_intersection; for (uint32 i = 0; i < rtm.num_of_spheres; ++i) { sse_ray sphere_pos; sphere_pos.x = _mm_set_ps1(rtm.spheres.x[i]); sphere_pos.y = _mm_set_ps1(rtm.spheres.y[i]); sphere_pos.z = _mm_set_ps1(rtm.spheres.z[i]); __m128 radius2 = _mm_set_ps1(rtm.spheres.radius[i]); radius2 = _mm_mul_ps(radius2, radius2); rlength.x = _mm_sub_ps(sphere_pos.x, rorigin.x); rlength.y = _mm_sub_ps(sphere_pos.y, rorigin.y); rlength.z = _mm_sub_ps(sphere_pos.z, rorigin.z); __m128 lendotdir = dot(rlength, rdir); no_intersection = _mm_cmplt_ps(lendotdir, _mm_set_ps1(0.0f)); __m128 d2 = _mm_sub_ps(dot(rlength, rlength), _mm_mul_ps(lendotdir, lendotdir)); no_intersection = _mm_or_ps(no_intersection, _mm_cmpgt_ps(d2, radius2)); if (_mm_test_all_ones(_mm_castps_si128(no_intersection))) continue; __m128 dprojrdir = _mm_sub_ps(radius2, d2); dprojrdir = _mm_sqrt_ps(dprojrdir); __m128 tnear = _mm_sub_ps(lendotdir, dprojrdir); __m128 tfar = _mm_add_ps(lendotdir, dprojrdir); __m128 tpoint = _mm_cmpgt_ps(tnear, _mm_set_ps1(0.0f)); tnear = _mm_blendv_ps(tfar, tnear, tpoint); tpoint = _mm_cmplt_ps(tnear, rs.tpoint); __m128 add_intersection = _mm_andnot_ps(no_intersection, tpoint); rs.tpoint = _mm_blendv_ps(rs.tpoint, tnear, add_intersection); rs.index = _mm_blendv_ps(rs.index, _mm_set_ps1(static_cast<float>(i)), add_intersection); } return rs; } ``` * fresnel effect 計算出新的原點及方向後遞迴呼叫trace_ray ```clike if (depth < TRACER_DEPTH) { __m128 facingratio = dot(rdir, nhit); facingratio = _mm_mul_ps(facingratio, _mm_set_ps1(-1.0f)); __m128 fresneleffect = _mm_sub_ps(_mm_set_ps1(1.0f), facingratio); fresneleffect = _mm_mul_ps(fresneleffect, fresneleffect); fresneleffect = _mm_mul_ps(fresneleffect, fresneleffect); const real mixing_value = 0.2f; fresneleffect = _mm_mul_ps(fresneleffect, _mm_set_ps1(1.0f - 0.2f)); fresneleffect = _mm_add_ps(fresneleffect, _mm_set_ps1(mixing_value)); sse_ray refldir; __m128 rdirdotnhit = dot(rdir, nhit); refldir.x = _mm_sub_ps(rdir.x, _mm_mul_ps(_mm_mul_ps(nhit.x, rdirdotnhit), _mm_set_ps1(2.0f))); refldir.y = _mm_sub_ps(rdir.y, _mm_mul_ps(_mm_mul_ps(nhit.y, rdirdotnhit), _mm_set_ps1(2.0f))); refldir.z = _mm_sub_ps(rdir.z, _mm_mul_ps(_mm_mul_ps(nhit.z, rdirdotnhit), _mm_set_ps1(2.0f))); normalize(refldir); sse_pixel reflpixel = trace_ray(phit, refldir, depth + 1); __m128 fresnel_r = _mm_mul_ps(reflpixel.r, fresneleffect); __m128 fresnel_g = _mm_mul_ps(reflpixel.g, fresneleffect); __m128 fresnel_b = _mm_mul_ps(reflpixel.b, fresneleffect); fresnel_r = _mm_mul_ps(sphere_color.r, _mm_andnot_ps(_mm_castsi128_ps(test_intersection), fresnel_r)); fresnel_g = _mm_mul_ps(sphere_color.g, _mm_andnot_ps(_mm_castsi128_ps(test_intersection), fresnel_g)); fresnel_b = _mm_mul_ps(sphere_color.b, _mm_andnot_ps(_mm_castsi128_ps(test_intersection), fresnel_b)); pixel.r = _mm_add_ps(pixel.r, fresnel_r); pixel.g = _mm_add_ps(pixel.g, fresnel_g); pixel.b = _mm_add_ps(pixel.b, fresnel_b); } } ``` ### multi-threaded * 利用 threadpool 將任務以 height_section = height / num_of_tasks 切割，測試在不同 thread 及 tasks 數量下的差異 ``` RayTracing Started Width: 1920 rays Height: 1020 rays Tracing 9 spheres Ray Tracer(0)[ Generated Camera Rays]: 0.0249563 s Ray Tracer(1)[ Raytracing Complete]: 0.190369 s Ray Tracer(2)[ BMP and PPM Images Saved to Disk]: 0.353309 s Have a nice day! ``` * 使用內建的模型來測試時間，有平行化的部份為Ray Tracer(1)[ Raytracing Complete] 我們探討這部份的時間 * 在 Task 為 32 下改變 thread 的數目，可以發現 2 thread 以後時間差異不大，這是受限於硬體的核心。![](https://i.imgur.com/cWbH7Y6.png) * 在 thread 為 2 下改變 task 的數目，在 32 個 task 下時間最小，但在 8 個task 以上的差異不大。 ![](https://i.imgur.com/2JOyWtL.png)