HackMD
Raytracing
contributed by <zhanyangch>, <hunng>, <yanang>
預期目標
- Task: 比照 B02: raytracing 要求，以 SIMD 改善 math-toolkit.h 定義的若干數學操作函式，進而推廣到整個光線追蹤的演算法 (可在不影響正確性的前提下，將原本的實作替換)。
- Task: 參考 Joshua Barczak 的 SIMD Raytracing 實作，改寫原有的程式碼，但是 model 應該和 B02: raytracing 一致
math-toolkit.h with SSE
Streaming SIMD Extensions：8個 128 位元暫存器，每個暫存器可以容納 2 個 64位元雙精度浮點數
Intrinsics：是一種外表類似一般的函數，但是實際上會被 compiler 直接編譯成組合語言的東西。
- 常用指令
__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 實做
__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);
}
執行測試
- 正確性
$ ./raytracing-sse
# Rendering scene
Done!
Execution time of raytracing() : 2.419169 sec
$ diff baseline.ppm out.ppm
- 時間比較
在此的 orig 為 loop unrolling 的版本
- 不使用 openmp
orig:2.247861 sec
sse:5.866375 sec
- 使用 openmp num_threads (2)
orig:1.289200 sec
sse:2.836198 sec
- 分析
使用 SSE 的版本反而比原版慢，推測可能的原因有
- SSE 的之資料長度為 128bit 運算 3個 double 需要切個成兩部份計算，造成指令數目增加
- 載入跟儲存的部份佔全部指令的比例太高，使用 _mm_load _mm_store ，取代 _mm_loadu _mm_storeu 應可以加快載入跟儲存的速度，但記憶體需對齊 16byte ，目前試著加入 attribute((aligned(16))但未能成功
- 修改 raytracing.c
為了解決載入跟儲存的部份佔全部指令的比例太高的問題，我們決定採用幾個合併多個功能的函數
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)
- 不使用 openmp
sse with additional functions: 4.545110 sec
- 使用 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 為例
__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
typedef double point3[4];
typedef double point4[4];
改用 _mm256_loadu_pd 作為載入資料手段
__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
__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()
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
__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
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
ymm2 = _mm256_permute2f128_pd (ymm3, ymm3, 131); // ymm0 = { v2^2, 0, 0, 0}
而 cross_product 因為使用了叫複雜的內容
所以所有值都是利用設定的
並為了讓整體運算結果會剛好符合順序
載入的順序就沒有那麼直觀
__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
進行測試
發現只有 normalize 與 length 有較佳的表現
Raytracing.c
- all function
- raySphereIntersection
- rayRectangularIntersection
- localColor
- compute_specular_diffuse
- reflection
- refraction
- fresnel
- ray_hit_object
- rayConstruction
- calculateBasisVectors
- protect_color_overflow
- ray_color
修改內容
為了使一次能運算 4 個 double 我們首先修改了 point3 的使用
改為使用 point3v
其內容為
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 為例
__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 為例，以下為原檔案
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);
}
}
修改過後
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 後面提到他使用到 Ray Reordering 的原因 -
而目前遇到的問題是
- 如何處理遞迴呼叫的問題
e.g., ray_color 會呼叫自己 - 且因為整體 struct 尚未設計完全，寫出的函式沒辦法測試，無法確認其正確性
- 如何處理遞迴呼叫的問題
其他 raytrace 實做
- 探討 jaybird19/simd-raytracer simd 切割及multi-threaded方式
LinearAllocator
allocator.init(main_arena_size + ALIGNMENT_PADDING, true);
先一次宣告大量的記憶體
pixmap = init_itr<real>((real*)allocator.allocate(pixels_heap_size, SSE_ALIGNMENT), pixels_heap_size);
為每一個資料分配位置
//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
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 進行計算。
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)在進行計算
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
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 以後時間差異不大，這是受限於硬體的核心。
-
在 thread 為 2 下改變 task 的數目，在 32 個 task 下時間最小，但在 8 個task 以上的差異不大。
