【python】手写光线追踪教程

最新推荐文章于 2023-01-11 11:13:51 发布
最新推荐文章于 2023-01-11 11:13:51 发布
阅读量3.4k 收藏 22
点赞数 10
文章标签： c# python 数学 图形学
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/qq_51214471/article/details/121360908
版权

请添加图片描述

本文未经允许禁止转载
B站：https://space.bilibili.com/455965619
作者：Heskey0

path tracer based on 《PBRT》

pbrt book

一.introduction to sampling theory

1. what is sampling?

impulse train：

冲激链

sampling process corresponds to multiplying the function by a “impulse train” function, an infinite sum of equally spaced delta functions.

采样

采样的图示

《PBRT》A digital image is represented as a set of pixel values, typically aligned on a rectangular grid. When a digital image is displayed on a physical device, these values are used to determine the spectral power emitted by pixels on the display.

《PBRT》the pixels that constitute an image are point samples of the image function at discrete points on the image plane.

there is no “area” associated with a pixel.

when sampling the film signal

pos = camera_pos
ray_dir = ti.Vector([
            (2 * fov * (u) / resolution[1] - fov * resolution[0] / resolution[1] - 1e-5),
            2 * fov * (v) / resolution[1] - fov - 1e-5, -1.0
        ]).normalized()

请添加图片描述

then we need anti-aliazing

pos = camera_pos
ray_dir = ti.Vector([
            (2 * fov * (u + ti.random()) / resolution[1] - fov * resolution[0] / resolution[1] - 1e-5),
            2 * fov * (v + ti.random()) / resolution[1] - fov - 1e-5, -1.0
        ]).normalized()

请添加图片描述

二.sampling

Preview (CDF sampling technique)

There are many techniques for generating random variates from a specified probability distribution such as the normal, exponential, or gamma distribution. However, one technique stands out because of its generality and simplicity: the inverse CDF sampling technique.

inverse CDF

1. Uniformly Sampling a Hemisphere (multidimensional sampling technique)

hemisphere sampling

a uniform distribution means that the density function is a constant, so we know that p(x) = c

半球均匀采样

so p(ω) = 1/2*pi

then p(θ, φ) = sinθ/2*pi

半球均匀采样-边缘密度

半球均匀采样-条件密度

Notice that the density function for φ itself is uniform

then use the 1D inversion technique to sample each of these PDFs in turn

半球均匀采样4

半球均匀采样5

半球均匀采样-Final

2. sample area light

sample light

def sample_area_light(hit_pos, pos_normal):
    # sampling inside the light area
    x = ti.random() * light_x_range + light_x_min_pos
    z = ti.random() * light_z_range + light_z_min_pos
    on_light_pos = ti.Vector([x, light_y_pos, z])
    return (on_light_pos - hit_pos).normalized()

请添加图片描述

3. introduction to importance sampling

why we need importance sampling?

the Monte Carlo estimator converges more quickly if the samples are taken from a distribution p(x) that is similar to the function f(x) in the integrand.

Monte Carlo estimator

《PBRT》：We will not provide a rigorous proof of this fact but will instead present an informal and intuitive argument.

then we try to analyze the importance sampling method

蒙特卡洛积分

we have three terms

  • BRDF
  • incident radiance ( infeasible )
  • cosine term
4. cosine-weighted sampling

1637050619568

1637050643929

1637050669913

Malley’s method

So, We could compute the marginal and conditional densities as before, but instead we can use a technique known as Malley’s method to generate these cosine-weighted points.

Malleys Method

  1. cosine term

  2. 2D Sampling with Multidimensional Transformations

    • (1) sampling a unit disk (Concentric Mapping)

    • (2) project up to the unit hemisphere (cosine-weighted hemisphere sampling)

(1) sampling a unit disk

1637051228888

1637051236921

(2) projection

To complete the (r,φ)=(sinθ,φ)⇒(θ,φ) transformation, we need the determinant of the Jacobian

1637051213404

Why1637051266186

1637051283108

5. multiple importance sampling

BDPT only：

请添加图片描述

BDPT + MIS：

1637051397587

Why we need MIS?

请添加图片描述

请添加图片描述

1637051701718

1637051712205

1637051724783

  • balance heuristic

1637051757112

  • power heuristic (Veach determined empirically that β=2 is a good value.)

1637051764336

//Compute heuristic
def mis_power_heuristic(pf, pg):
    # Assume 1 sample for each distribution
    f = pf ** 2
    g = pg ** 2
    return f / (f + g)
    # return 1

//combine
@ti.func
def sample_light_and_cosineWeighted(hit_pos, hit_normal):
    cosine_by_pdf = ti.Vector([0.0, 0.0, 0.0])

    light_pdf, cosineWeighted_pdf = 0.0, 0.0

    # sample area light => dir, light_pdf; then dir => lambertian_pdf; then mis
    light_dir = sample_area_light(hit_pos, hit_normal)
    if light_dir.dot(hit_normal) > 0:
        light_pdf = compute_area_light_pdf(hit_pos, light_dir)
        cosineWeighted_pdf = compute_cosineWeighted_pdf(hit_normal, light_dir)
        if light_pdf > 0 and cosineWeighted_pdf > 0:
            l_visible = visible_to_light(hit_pos, light_dir)
            if l_visible:
                heuristic = compute_heuristic(light_pdf, cosineWeighted_pdf)
                DoN = dot_or_zero(light_dir, hit_normal)
                cosine_by_pdf += heuristic * DoN / light_pdf

    # sample cosine weighted => dir, lambertian_pdf; then dir => light_pdf; then mis
    cosineWeighted_dir = cosine_weighted_sampling(hit_normal)
    cosineWeighted_pdf = compute_cosineWeighted_pdf(hit_normal, cosineWeighted_dir)
    light_pdf = compute_area_light_pdf(hit_pos, cosineWeighted_dir)
    if visible_to_light(hit_pos, cosineWeighted_dir):
        heuristic = compute_heuristic(cosineWeighted_pdf, light_pdf)
        DoN = dot_or_zero(cosineWeighted_dir, hit_normal)
        cosine_by_pdf += heuristic * DoN / cosineWeighted_pdf

    # direct_li = mis_weight * cosine / pdf
    return cosine_by_pdf

import taichi as ti
import numpy as np

ti.init(arch=ti.cuda)
resolution = (940, 940)

eps = 0.0001  # 浮点数精度
inf = 1e10

mat_none = 0
mat_lambertian = 1
mat_specular = 2    # 镜面
mat_glass = 3       # 玻璃
mat_light = 4
mat_microfacet = 5
mat_glossy = 6

# 光区域为一块板
light_y_pos = 2.0 - eps
light_x_min_pos = -0.7
light_x_range = 1.4
light_z_min_pos = 0.6
light_z_range = 0.4
light_area = light_x_range * light_z_range
light_min_pos = ti.Vector([
    light_x_min_pos,
    light_y_pos,
    light_z_min_pos])
light_max_pos = ti.Vector([
    light_x_min_pos + light_x_range,
    light_y_pos,
    light_z_min_pos + light_z_range
])
light_color = ti.Vector([1, 1, 1])
light_normal = ti.Vector([0.0, -1.0, 0.0])  # 光源方向向下


# 1.7700 : 红宝石的折射率
refract_index = 1.7700

# right sphere
sp1_center = ti.Vector([0.5, 1.18, 1.40])
sp1_radius = 0.18
# left sphere
sp2_center = ti.Vector([-0.35, 0.65, 1.70])
sp2_radius = 0.15
# middle sphere(microfacet)
sp3_center = ti.Vector([-0.10, 0.35, 1])
sp3_radius = 0.35
sp3_microfacet_roughness = 0.5
# sp3_idx = 1.55  # 石英晶体折射率
sp3_idx = 2.4  # 钻石折射率
# right front sphere(microfacet)
sp4_center = ti.Vector([-0.05, 1, 1])
sp4_radius = 0.3
sp4_microfacet_roughness = 1


# 构造变换矩阵，用于box
def make_box_transform_matrices(rotate_rad, translation):
    c, s = np.cos(rotate_rad), np.sin(rotate_rad)
    rot = np.array([[c, 0, s, 0],
                    [0, 1, 0, 0],
                    [-s, 0, c, 0],
                    [0, 0, 0, 1]])  # 绕y轴旋转67.5°
    # rot = np.array([[1, 0, 0, 0],
    #                 [0, c, s, 0],
    #                 [0,-s, c, 0],
    #                 [0, 0, 0, 1]])  # 绕y轴旋转67.5°
    translate = np.array([  # 平移 (0.5, 0, 1.4)
        [1, 0, 0, translation.x],
        [0, 1, 0, translation.y],
        [0, 0, 1, translation.z],
        [0, 0, 0, 1],
    ])
    m = translate @ rot  # 平移 + 旋转
    m_inv = np.linalg.inv(m)  # 逆矩阵
    m_inv_t = np.transpose(m_inv)  # 转置矩阵
    return ti.Matrix(m_inv), ti.Matrix(m_inv_t)  # 旋转-22.5° + 平移 (0.5, 0, 1)


# right box
box1_min = ti.Vector([0.0, 0.0, 0.0])
box1_max = ti.Vector([0.35, 1.0, 0.35])
box1_rotate_rad = np.pi / 16
box1_m_inv, box1_m_inv_t = make_box_transform_matrices(box1_rotate_rad, ti.Vector([0.30, 0, 1.20]))  # box的transform的 逆矩阵, 逆转置矩阵
# left box
box2_min = ti.Vector([0.0, 0.0, 0.0])
box2_max = ti.Vector([0.4, 0.5, 0.4])
box2_rotate_rad = np.pi / 4
box2_m_inv, box2_m_inv_t = make_box_transform_matrices(box2_rotate_rad, ti.Vector([-0.75, 0, 1.70]))  # box的transform的 逆矩阵, 逆转置矩阵


'''
lambertian brdf
'''
# No absorbtion     没有吸收光谱，Albedo为1，对单位半球积分
lambertian_brdf = 1.0 / np.pi  # f(lambert) = k*c / π       # k = 1,  c = hit_color*light_color


'''
microfacet brdf
'''
# compute reflectance
# 计算反射比
@ti.func
def schlick(cos, eta):  # 入射角cosine, 折射率refractive index
    r0 = (1.0 - eta) / (1.0 + eta)
    r0 = r0 * r0  # 反射比 reflectance
    return r0 + (1 - r0) * ((1.0 - cos) ** 5)


# normal distribution function
@ti.func
def ggx(alpha, i_dir, o_dir, n_dir):  # roughness, incident, exit, normal
    m_dir = (i_dir + o_dir).normalized()
    cos_theta_square = m_dir.dot(n_dir)
    tan_theta_square = (1-cos_theta_square) / cos_theta_square
    root = alpha / cos_theta_square * (alpha*alpha + tan_theta_square)
    return root*root / np.pi


@ti.func
def ggx2(alpha, i_dir, o_dir, n_dir):
    m_dir = (i_dir + o_dir).normalized()
    NoM = n_dir.dot(m_dir)
    d = NoM*NoM * (alpha*alpha-1) + 1
    return alpha*alpha / np.pi*d*d


@ti.func
def smithG1(alpha, v_dir, n_dir):
    out = 0.0
    # compute tan_theta(v / n)
    cos_theta_square = v_dir.dot(n_dir) ** 2
    tan_theta_square = (1-cos_theta_square) / cos_theta_square
    tan_theta = ti.sqrt(tan_theta_square)
    if tan_theta == 0:
        out = 1
    else:
        root = alpha * tan_theta
        out = 2 / (1 + ti.sqrt(1.0 + root * root))

    return out


@ti.func
# shadowing-masking
def smith(alpha, i_dir, o_dir, n_dir):  # roughness, incident, exit, normal
    # m_dir = (i_dir + o_dir).normalized()
    # shadowing * masking
    return smithG1(alpha, i_dir, n_dir) * smithG1(alpha, o_dir, n_dir)


@ti.func
def compute_microfacet_brdf(alpha, idx, i_dir, o_dir, n_dir):
    micro_cos = o_dir.dot((i_dir + o_dir).normalized())
    # numerator and denominator
    D = ggx2(alpha, i_dir, o_dir, n_dir)
    G = smith(alpha, i_dir, o_dir, n_dir)
    F = schlick(micro_cos, idx)
    # print(D, G, F)

    numerator = D * G * F
    denominator = 4 * o_dir.dot(n_dir) * i_dir.dot(n_dir)
    cook_torrance = numerator / ti.abs(denominator)
    return cook_torrance


'''
basic functions
'''
# 反射
@ti.func
def reflect(d, n):
    # d and n are both normalized
    ret = d - 2.0 * d.dot(n) * n  # d - 2*|d|*|n|*n*cos<d,n>(theta) = d - 2 |d|*cos(theta) * (n/|n|)
    return ret  # reflect vector


# 折射
@ti.func
def refract(d, n, ni_over_nt):
    dt = d.dot(n)  # cos    # sin**2 = 1 - cos**2
    discr = 1.0 - ni_over_nt * ni_over_nt * (1.0 - dt * dt)  # discr:折射角的cos
    rd = (ni_over_nt * (d - n * dt) - n * ti.sqrt(discr)).normalized()
    return rd  # 是否有反射光, 反射光方向


# 点由矩阵变换
@ti.func
def mat_mul_point(m, p):
    hp = ti.Vector([p[0], p[1], p[2], 1.0])
    hp = m @ hp
    hp /= hp[3]
    return ti.Vector([hp[0], hp[1], hp[2]])


# [3] => ti.Vector(4);   m@v  # [4, 4]@[4]
# 忽略矩阵的第4行第4列, 忽略矩阵的平移
@ti.func
def mat_mul_vec(m, v):
    hv = ti.Vector([v[0], v[1], v[2], 0.0])
    hv = m @ hv
    return ti.Vector([hv[0], hv[1], hv[2]])


# 判断射线与球是否相交
@ti.func
def intersect_sphere(pos, d, center, radius):  # pos:light_position, d:ray_dir
    # 构建余弦定理三角形:判断光与球是否相交
    T = pos - center
    A = 1.0
    B = 2.0 * T.dot(d)
    C = T.dot(T) - radius * radius
    delta = B * B - 4.0 * A * C
    dist = inf
    hit_pos = ti.Vector([0.0, 0.0, 0.0])

    if delta > 0:  # 有解
        delta = ti.max(delta, 0)
        sdelta = ti.sqrt(delta)
        ratio = 0.5 / A
        ret1 = ratio * (-B - sdelta)  # 方程的解, 即三角形的边长(离入射光近的点)
        dist = ret1
        hit_pos = pos + d * dist

    return dist, hit_pos  # 光源到命中点的距离, 命中点坐标


# plane
@ti.func
def intersect_plane(pos, d, pt_on_plane, norm):  # position, ray_dir, offset, normal
    dist = inf
    hit_pos = ti.Vector([0.0, 0.0, 0.0])
    denom = d.dot(norm)
    if abs(denom) > eps:  # 光与平面不平行
        dist = norm.dot(pt_on_plane - pos) / denom
        hit_pos = pos + d * dist
    return dist, hit_pos  # 光源到命中点的距离, 命中点坐标


# 参考清华大学图形学课程中的基于slab的求交算法：Liang_Barsky算法
# aabb包围体 call by intersect_box and intersect_light
@ti.func
def intersect_aabb(box_min, box_max, o, d):  # box_min, box_max, pos(box空间), ray_dir(box空间)
    intersect = 1  # 光与box是否相交

    near_t = -inf
    far_t = inf
    near_face = 0
    near_is_max = 0

    for i in ti.static(range(3)):  # ti.static(range()) can iterate matrix elements
        if d[i] == 0:  # 光平行于包围体的一个面
            if o[i] < box_min[i] or o[i] > box_max[i]:
                intersect = 0
        else:
            i1 = (box_min[i] - o[i]) / d[i]  # 除以d[i] : 判断光是否正对box
            i2 = (box_max[i] - o[i]) / d[i]

            new_far_t = max(i1, i2)     # 光朝着正半轴时，为i2
            new_near_t = min(i1, i2)    # 光朝着正半轴时，为i1
            new_near_is_max = i2 < i1   # 光朝着负半轴时(near_t取i2)，为true

            far_t = min(new_far_t, far_t)  # far_t 取最小
            if new_near_t > near_t:  # near_t 取最大
                near_t = new_near_t
                near_face = int(i)  # 记录最小的i所在的维
                near_is_max = new_near_is_max  # 在当前维中near_t, i2<i1 ?

    near_norm = ti.Vector([0.0, 0.0, 0.0])
    if near_t > far_t:
        intersect = 0
    if intersect:
        for i in ti.static(range(2)):
            if near_face == i:
                near_norm[i] = -1 + near_is_max * 2     # near_is_max => return 1; else => return -1
    return intersect, near_t, far_t, near_norm  # 是否相交, 首先相交的平面的距离, 远平面, 近平面法线


# params: min, max, position, ray_dir
# box
@ti.func
def intersect_aabb_transformed(box_m_inv, box_m_inv_t, box_min, box_max, o, d):
    # 射线转换到包围体的local position
    obj_o = mat_mul_point(box_m_inv, o)
    obj_d = mat_mul_vec(box_m_inv, d)

    intersect, near_t, _, near_norm = intersect_aabb(box_min, box_max, obj_o, obj_d)
    # print(near_norm)
    if intersect and 0 < near_t:
        near_norm = mat_mul_vec(box_m_inv_t, near_norm)
    else:
        intersect = 0
    # out params: hit?, cur_dist, pnorm
    return intersect, near_t, near_norm


# light
@ti.func
def intersect_light(pos, ray_dir, tmax):
    # t:near intersect distance
    hit, t, far_t, near_norm = intersect_aabb(light_min_pos, light_max_pos, pos, ray_dir)
    if hit and 0 < t < tmax:
        hit = 1
    else:
        hit = 0
        t = inf
    return hit, t


# 光线与场景相交
@ti.func
def intersect_scene(pos, ray_dir):
    # closest:深度缓冲区
    closest, normal = inf, ti.Vector.zero(ti.f32, 3)
    # color, material
    c, mat = ti.Vector.zero(ti.f32, 3), mat_none

    # right sphere
    cur_dist, hit_pos = intersect_sphere(pos, ray_dir, sp1_center, sp1_radius)
    if 0 < cur_dist < closest:  # 深度测试
        closest = cur_dist
        normal = (hit_pos - sp1_center).normalized()
        c, mat = ti.Vector([1.0, 1.0, 1.0]), mat_glass

    # middle Sphere
    cur_dist, hit_pos = intersect_sphere(pos, ray_dir, sp3_center, sp3_radius)
    if 0 < cur_dist < closest:  # 深度测试
        closest = cur_dist
        normal = (hit_pos - sp3_center).normalized()
        c, mat = ti.Vector([102.0/255.0, 153.0/255.0, 255.0/255.0]), mat_microfacet

    # left Sphere
    cur_dist, hit_pos = intersect_sphere(pos, ray_dir, sp2_center, sp2_radius)
    if 0 < cur_dist < closest:  # 深度测试
        closest = cur_dist
        normal = (hit_pos - sp2_center).normalized()
        c, mat = ti.Vector([1.0, 1.0, 1.0]), mat_specular

    # left box
    hit, cur_dist, pnorm = intersect_aabb_transformed(box2_m_inv, box2_m_inv_t, box2_min, box2_max, pos, ray_dir)
    if hit and 0 < cur_dist < closest:  # 深度测试
        closest = cur_dist
        normal = pnorm
        c, mat = ti.Vector([0.8, 1, 1]), mat_lambertian

    # right box
    hit, cur_dist, pnorm = intersect_aabb_transformed(box1_m_inv, box1_m_inv_t, box1_min, box1_max, pos, ray_dir)
    if hit and 0 < cur_dist < closest:  # 深度测试
        closest = cur_dist
        normal = pnorm
        c, mat = ti.Vector([0.8, 1, 1]), mat_lambertian

    # left plane
    pnorm = ti.Vector([1.0, 0.0, 0.0])
    cur_dist, _ = intersect_plane(pos, ray_dir, ti.Vector([-1.1, 0.0, 0.0]), pnorm)
    if 0 < cur_dist < closest:  # 深度测试
        closest = cur_dist
        normal = pnorm
        c, mat = ti.Vector([60.0 / 255.0, 200.0 / 255.0, 60 / 255.0]), mat_lambertian
    # right plane
    pnorm = ti.Vector([-1.0, 0.0, 0.0])
    cur_dist, _ = intersect_plane(pos, ray_dir, ti.Vector([1.1, 0.0, 0.0]), pnorm)
    if 0 < cur_dist < closest:  # 深度测试
        closest = cur_dist
        normal = pnorm
        c, mat = ti.Vector([200.0 / 255.0, 30.0 / 255.0, 30 / 255.0]), mat_lambertian
    # bottom plane
    gray = ti.Vector([0.93, 0.93, 0.93])
    pnorm = ti.Vector([0.0, 1.0, 0.0])
    cur_dist, _ = intersect_plane(pos, ray_dir, ti.Vector([0.0, 0.0, 0.0]), pnorm)
    if 0 < cur_dist < closest:  # 深度测试
        closest = cur_dist
        normal = pnorm
        c, mat = gray, mat_lambertian
    # top
    pnorm = ti.Vector([0.0, -1.0, 0.0])
    cur_dist, _ = intersect_plane(pos, ray_dir, ti.Vector([0.0, 2.0, 0.0]), pnorm)
    if 0 < cur_dist < closest:  # 深度测试
        closest = cur_dist
        normal = pnorm
        c, mat = gray, mat_lambertian
    # far
    pnorm = ti.Vector([0.0, 0.0, 1.0])
    cur_dist, _ = intersect_plane(pos, ray_dir, ti.Vector([0.0, 0.0, 0.0]), pnorm)
    if 0 < cur_dist < closest:  # 深度测试
        closest = cur_dist
        normal = pnorm
        c, mat = gray, mat_lambertian
    # close
    pnorm = ti.Vector([0.0, 0.0, -1.0])
    cur_dist, _ = intersect_plane(pos, ray_dir, ti.Vector([0.0, 0.0, 3]), pnorm)
    if 0 < cur_dist < closest:  # 深度测试
        closest = cur_dist
        normal = pnorm
        c, mat = ti.Vector([0, 0, 0]), mat_lambertian

    # light
    hit_l, cur_dist = intersect_light(pos, ray_dir, closest)
    if hit_l and 0 < cur_dist < closest:  # 深度测试
        # no need to check the second term
        closest = cur_dist
        normal = light_normal
        c, mat = light_color, mat_light

    return closest, normal, c, mat


# 判断ray_dir是否与光源相交
@ti.func
def visible_to_light(pos, ray_dir):
    # eps*ray_dir to prevent rounding error
    a, b, c, mat = intersect_scene(pos + eps * ray_dir, ray_dir)
    return mat == mat_light


@ti.func
def dot_or_zero(n, l):
    return max(0.0, n.dot(l))











# TODO:begin
# '''
# sampling functions

# multiple importance sampling
@ti.func
def compute_heuristic(pf, pg):
    # Assume 1 sample for each distribution
    f = pf ** 2
    g = pg ** 2
    return f / (f + g)


# 已知sample dir
# area light pdf
@ti.func
def compute_area_light_pdf(pos, ray_dir):
    hit_l, t = intersect_light(pos, ray_dir, inf)
    pdf = 0.0
    if hit_l:  # ray_dir命中了灯光
        l_cos = light_normal.dot(-ray_dir)  # 光源的方向 与 ray_dir 的夹角cosine
        if l_cos > eps:  # 光源 与 ray_dir 同向
            tmp = ray_dir * t
            dist_sqr = tmp.dot(tmp)
            pdf = dist_sqr / (light_area * l_cos)
    return pdf


# 已知sample dir
# cosine weighted sampling
@ti.func
def compute_cosineWeighted_pdf(normal, sample_dir):
    return dot_or_zero(normal, sample_dir) / np.pi  # p(theta, phi) = cos(theta) * sin(theta) / pi


# 未知sample dir
# sample light
@ti.func
def sample_area_light(hit_pos, pos_normal):
    # sampling inside the light area
    x = ti.random() * light_x_range + light_x_min_pos
    z = ti.random() * light_z_range + light_z_min_pos
    on_light_pos = ti.Vector([x, light_y_pos, z])
    return (on_light_pos - hit_pos).normalized()


# 未知sample dir
# Cosine-Weighted Sampling
@ti.func
def cosine_weighted_sampling(normal):
    r, phi = 0.0, 0.0  # 圆上的 (r, theta) 在半球里实际上是 (sin(theta), phi) ，将其变换到 (theta, phi)
    sx = ti.random() * 2.0 - 1.0  # -1 ~ 1 random
    sy = ti.random() * 2.0 - 1.0  # -1 ~ 1 random
    # 1.concentric sample
    # sample on a unit disk
    if sx != 0 or sy != 0:
        if abs(sx) > abs(sy):
            r = sx
            phi = np.pi / 4 * (sy / sx)
        else:
            r = sy
            phi = np.pi / 4 * (2 - sx / sy)

    # 2.apply Malley's method
    # project disk to hemisphere

    # 由normal为中心轴,u和v为水平轴建立笛卡尔坐标系
    # 不需要关心normal和vector.up的关系，vector.up的引入是为了辅助建立起坐标系(u,v,normal)
    u = ti.Vector([1.0, 0.0, 0.0])
    if abs(normal[1]) < 1 - eps:
        u = normal.cross(ti.Vector([0.0, 1.0, 0.0]))  # normal x vector.up = sin(eta)
    v = normal.cross(u)  # normal x u = |u| = sin(eta)

    # theta : vector.up 与 normal 的夹角
    # u,v垂直, 长度均为sin(phi), 均在微平面上

    xy = r * ti.cos(phi) * u + r * ti.sin(phi) * v  # 采样时的x,y,normal坐标系转换到u,v,normal坐标系(采样点随之旋转并变为sin(eta)倍)
    zlen = ti.sqrt(max(0.0, 1.0 - xy.dot(xy)))  # zlen:采样线沿normal的长度

    return xy + zlen * normal  # sample dir


# 两种pdf相乘, 结果为对光采样
# sample direct light
@ti.func
def sample_light_and_cosineWeighted(hit_pos, hit_normal):
    cosine_by_pdf = ti.Vector([0.0, 0.0, 0.0])

    light_pdf, cosineWeighted_pdf = 0.0, 0.0

    # sample area light => dir, light_pdf; then dir => lambertian_pdf; then mis
    light_dir = sample_area_light(hit_pos, hit_normal)
    if light_dir.dot(hit_normal) > 0:
        light_pdf = compute_area_light_pdf(hit_pos, light_dir)
        cosineWeighted_pdf = compute_cosineWeighted_pdf(hit_normal, light_dir)
        if light_pdf > 0 and cosineWeighted_pdf > 0:
            l_visible = visible_to_light(hit_pos, light_dir)
            if l_visible:
                heuristic = compute_heuristic(light_pdf, cosineWeighted_pdf)
                DoN = dot_or_zero(light_dir, hit_normal)
                cosine_by_pdf += heuristic * DoN / light_pdf

    # sample cosine weighted => dir, lambertian_pdf; then dir => light_pdf; then mis
    cosineWeighted_dir = cosine_weighted_sampling(hit_normal)
    cosineWeighted_pdf = compute_cosineWeighted_pdf(hit_normal, cosineWeighted_dir)
    light_pdf = compute_area_light_pdf(hit_pos, cosineWeighted_dir)
    if visible_to_light(hit_pos, cosineWeighted_dir):
        heuristic = compute_heuristic(cosineWeighted_pdf, light_pdf)
        DoN = dot_or_zero(cosineWeighted_dir, hit_normal)
        cosine_by_pdf += heuristic * DoN / cosineWeighted_pdf

    # direct_li = mis_weight * cosine / pdf
    return cosine_by_pdf



@ti.func
def sample_ray_dir(indir, normal, hit_pos, mat):
    u = ti.Vector([0.0, 0.0, 0.0])  # 用于下一次追踪的ray_dir
    pdf = 1.0
    if mat == mat_lambertian:
        u = cosine_weighted_sampling(normal)  # sample brdf : return ray_dir
        pdf = max(eps, compute_cosineWeighted_pdf(normal, u))  # 计算在该方向采样射线的pdf
    elif mat == mat_glossy:
        pass
    elif mat == mat_microfacet:
        # TODO:对cosine项采样
        u = cosine_weighted_sampling(normal)  # sample brdf : return ray_dir
        pdf = max(eps, compute_cosineWeighted_pdf(normal, u))  # 计算在该方向采样射线的pdf

    elif mat == mat_specular:  # 反射, pdf = 1
        u = reflect(indir, normal)
    elif mat == mat_glass:  # 折射, 反射, pdf = 1
        cos = indir.dot(normal)  # indir和normal的夹角 (indir和normal为单位向量)
        ni_over_nt = refract_index  # ni / nt = 折射率
        outn = normal
        if cos > 0.0:
            outn = -normal
            cos = refract_index * cos  # 出射角度
        else:
            ni_over_nt = 1.0 / refract_index
            cos = -cos  # indir转180°

        refl_prob = schlick(cos, refract_index)  # Fresnel reflectance
        if ti.random() < refl_prob:  # 反射的能量
            u = reflect(indir, normal)
        else:  # 折射的能量
            u = refract(indir, outn, ni_over_nt)
    return u.normalized(), pdf  # 用于下一次追踪的ray_dir, pdf


pixels = ti.Vector.field(3, dtype=ti.f32, shape=resolution)

camera_pos = ti.Vector([0.0, 0.6, 3.0])
fov = 0.8

max_bounce = 10


@ti.kernel
def render():
    for u, v in pixels:  # 遍历像素
        pos = camera_pos
        ray_dir = ti.Vector([
            (2 * fov * (u + ti.random()) / resolution[1] - fov * resolution[0] / resolution[1] - 1e-5),
            2 * fov * (v + ti.random()) / resolution[1] - fov - 1e-5, -1.0
        ]).normalized()

        final_throughput = ti.Vector([0.0, 0.0, 0.0])  # 累加到pixels
        throughput = ti.Vector([1.0, 1.0, 1.0])  # Lighting : (r, g, b)

        # 追踪开始
        bounce = 0
        while bounce < max_bounce:  # bounce的最大次数
            bounce += 1
            # closest:光源到物体的距离
            closest, hit_normal, hit_color, mat = intersect_scene(pos, ray_dir)  # 光发出后碰到场景

            # 0.命中灯光或无材质, 则中断追踪
            if mat == mat_none:
                final_throughput += throughput * 0
                break
            if mat == mat_light:
                final_throughput += throughput * light_color
                break

            hit_pos = pos + closest * ray_dir
            ray_dir_i = -ray_dir

            # 1.计算采样后的ray_dir, pdf

            # 2.lambertian : sample direct light [ mis(sample area light, sample brdf)=> Li ]
            if mat == mat_lambertian:  # lambertian模型
                final_throughput += light_color * throughput * lambertian_brdf * hit_color * sample_light_and_cosineWeighted(hit_pos, hit_normal)
                # Sample Direct Light Only
                # throughput *= sample_light_and_cosineWeighted(hit_pos, hit_normal, hit_color)

            # 2.lambertian : sample cosine-Weighted
            ray_dir, pdf = sample_ray_dir(ray_dir, hit_normal, hit_pos, mat)  # 由反射更新ray_dir
            pos = hit_pos + eps * ray_dir
            if mat == mat_lambertian:  # lambertian
                # f(lambert) * max(0.0, cos(n,l)) / pdf
                # throughput : Li or Lo
                throughput *= (lambertian_brdf * hit_color) * dot_or_zero(hit_normal, ray_dir) / pdf

            # 3.specular全反射
            if mat == mat_specular:
                throughput *= hit_color
            # 4.glass折射btdf
            if mat == mat_glass:
                throughput *= hit_color

            # 5.microfacet
            if mat == mat_microfacet:
                # compute_microfacet_brdf params:(alpha, idx, i_dir, o_dir, n_dir)
                cook_torrance_brdf = compute_microfacet_brdf(sp3_microfacet_roughness, sp3_idx, ray_dir_i, ray_dir, hit_normal)
                # print(lambertian_brdf, cook_torrance_brdf)

                microfacet_brdf = lambertian_brdf + cook_torrance_brdf  # TODO:BUG 黑屏

                throughput *= (microfacet_brdf * hit_color) * dot_or_zero(hit_normal, ray_dir) / pdf

            # 6.glossy
            if mat == mat_glossy:  #
                throughput *= (lambertian_brdf * hit_color) * dot_or_zero(hit_normal, ray_dir) / pdf




        pixels[u, v] += final_throughput


gui = ti.GUI('Path Tracing', resolution)
i = 0

while gui.running:
    # if gui.get_event(ti.GUI.PRESS):
    #     if gui.event.key == 'w':
    #         gui.clear()
    #         i = 0
    #         interval = 10
    #         # pixels = ti.Vector.field(3, dtype=ti.f32, shape=resolution)  # 屏幕像素缓冲 [800, 800] 元素为(r, g, b)
    #         count_var = ti.field(ti.i32, shape=(1,))
    #         box1_rotate_rad += np.pi/8

    if gui.get_event(ti.GUI.PRESS):
        if gui.event.key == 'w':
            img = pixels.to_numpy()
            img = np.sqrt(img / img.mean() * 0.24)
            fname = f'cornell_box.png'
            ti.imwrite(img, fname)
            print("图片已存储")

    render()
    interval = 10  # render()10次, 绘1次图
    if i % interval == 0 and i > 0:
        img = pixels.to_numpy()
        img = np.sqrt(img / img.mean() * 0.24)
        gui.set_image(img)

        gui.show()
    i += 1

本文未经允许禁止转载
B站：https://space.bilibili.com/455965619
作者：Heskey0

Heskey0
关注
  • 10
    点赞
  • 踩
  • 22
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
pyrate:基于Python的光学光线追踪
05-15
y酸盐 使用Python进行光学设计。 介绍 我们的目标是为各向同性，均质各向异性和非均质各向同性GRIN介质提供一种易于使用的光学光线示踪剂，该介质可以主要由python脚本控制，并提供足够的接口类和功能以集成到GUI中。 要求 您需要安装NumPy，SciPy，Yaml，Sympy和matplotlib的Python 2.7或3.x才能运行吡喃酸酯。 在Ubuntu，Mint和Debian中，您可以使用Python 2.7 $ sudo apt-get install python python-numpy python-scipy python-matplotlib python-yaml python-sympy 或者 $ sudo apt-get install python3 python3-numpy python3-scipy python3-matplotlib
探索光学新纪元：RayTracing——简单易用的Python光线追踪模块
最新发布
gitblog_00004的博客
06-07 334
探索光学新纪元：RayTracing——简单易用的Python光线追踪模块 项目地址:https://gitcode.com/DCC-Lab/RayTracing 在数字时代，理解并模拟光的行为是科学研究与工程设计中的重要环节。而今天，我们向您隆重推荐一个专为教学和研究定制的开源项目——RayTracing，这是一个基于Python的简单且实用的光线追踪工具。 项目介绍 RayTracing由DC...
用python模拟光学追迹，并给出示例
weixin_42576467的博客
01-11 848
在 Python 中可以使用各种库来模拟光学追迹。其中常用的有 RayTracing、OpticalRayTracer 和 PyTracer。 下面是一个使用 RayTracing 库进行简单光学追迹的示例代码： from raytracing import Scene, Sphere, PointLight # 创建场景scene = Scene() # 添加球体 sphere = Sphe...
pyoptics:用于教育和小型项目设计的python中的几何光学光线追踪
05-13
焦光学 用于教育和小型项目设计的python中的几何光学光线追踪。 这是一个简单的，具有alpha质量的项目，没有当前的未来开发计划。 另请参阅：GeOptics 特征 使用小型光学零件库（镜头，镜子，棱镜）进行光线追踪 Shott和其他用于色散建模的常见玻璃的目录 用于交互式放置光学元件和可视化光线追踪输出的GUI 例子 在这里学习内容的最快方法是阅读和/或运行示例： dispersion.py ：将白光分解为彩虹的棱镜 mirror.py ：反射平行光线的镜子的简单组合 2p_telescope.py ：在激光扫描显微镜设置中测试焦点的波长依赖性 2p_telescope2.py ：激光扫描显微镜设置的动画模拟 使用鼠标平移/缩放模型。
用Python和OpenGL实现体光线投射算法（Win11）
Alexabc3000的专栏
08-27 724
用Python和OpenGL实现体光线投射算法
学习光线追踪(12)---画一面墙
ARTELE的博客
05-17 415
0.简介 最近事情很多，这算是才有些空写写，先拿一个简单的缓一缓，不然之前写的都快忘了。上次画出了三角形，又添加了纹理，这次画一个墙面做背景。 1.墙面 一面墙就是一个矩形平面，两个三角形拼接起来就可以了，但是三角形一个一个的声明比较麻烦，所以就新建一个类，专门绘制这种多个三角形组合。 class Triangles : public Polygon { vector<Trian...
逆采样(Inverse Sampling)和拒绝采样(Reject Sampling)原理详解
anshuai_aw1的博客
12-05 1万+
在蒙特卡罗方法中，有一个关键的问题需要解决，即如何基于概率密度函数去采的nnn个xxx的样本集。 逆采样(Inverse Sampling)和拒绝采样(Reject Sampling)就是用于解决这个问题的。 其中，对于常见的分布，如均匀分布，高斯分布，指数分布，t分布，F分布，Beta分布，Gamma分布等，可以采用逆采样的方法进行采样；不过很多时候，我们的xxx的概率分布不是常见的分布，这些分...
lighttrace_光线追踪_python_
09-29
在这个“lighttrace_光线追踪_python_”项目中，我们将专注于使用Python编程语言实现一个简单的二维光线追踪程序。 光线追踪的核心原理是模拟光线从眼睛出发，通过虚拟相机系统进入场景，与场景中的各个物体相互...
基于python生成手写字体
12-11
使用hangright库生成手写字 也可用于艺术字生成 也可以仿照他人的手写笔记生成字体
光线追踪算法_100行python代码pdf.pdf
12-27
光线追踪算法_100行python代码pdf.pdf 光线追踪算法（Ray Tracing）是一种用于生成照片级真实感图像的技术。该技术通过模拟光线在三维场景中的行为，生成高质量的图像。以下是关于光线追踪算法的知识点： 1. 场景...
光线追踪算法实现C++实现
09-03
课程作业写的光线追踪算法 包含说明文档C++代码编写场景文件包含 球体 三角形 平面 光源 材料质 镜面 漫反射
Python-用PyTorch实现的光线追踪
08-11
用PyTorch实现的光线追踪
基于python的三维射线追踪库-ttcrpy详解（1）
Neverlevsun的博客
04-26 5114
三维射线追踪的python库——ttcrpy详解（1）
【Ray Tracing】光线追踪基本概念与代码实现
热门推荐
北境の守望者
08-09 2万+
Ref Ref 视觉技术的圣杯：光线追踪如何再现真实世界？ : 超级棒的综述报道，一文入坑 用JavaScript玩转计算机图形学(一)光线追踪入门:Milo Yip 大神的经典之作，带JS交互 一个系列，按 Milo Yip 的思路编排，使用 C++ 实现，带源码 光线追踪(RayTracing)算法理论与实践(一)入门 光线追踪(RayTracing)算法理论与实践(二...
【计算机图形学】从0实现光线追踪 · 后篇
丧心病狂Loli控的博客
12-14 533
回顾 在上一节中，我们已经完成了： 架了相机(Camera)，摆放好场景(Scene)； 射出射线、计算交点(IntersectResult)； 实现了渲染器(Render)，即将交点信息返回到画布(Canvas)上   接下来，我们将完成： 手写材质(Material)，其中，材质最重要的方法是取样(sample)，即获得颜色 手写光线追踪过程，这是一个递归 渲染出成品       颜色（Color） 计算出交点处的颜色(color)，返回并绘制到画布
问题十二：怎么用ray tracing画第一张图
图形跟班
01-13 2539
ray 类：（位于ray.h中。ray.cpp为空） #ifndef RAY_H #define RAY_H #include "vec3.h" class ray { public: ray() {} ray(const vec3& a, const vec3& b) { A = a; B = b; } vec3 orgin() con
python光线追踪
08-14
光线追踪（Ray Tracing）是一种用于生成逼真图像的渲染技术。在光线追踪中，通过跟踪光线在场景中的交互，来计算光线与物体之间的相互作用，从而生成最终的图像。 在 Python 中，有一些库可以用来实现光线追踪，例如 PyRay 和 PyTrace。这些库提供了一些基本的功能，如创建场景、定义光源和物体、计算光线与物体的交点等。 下面是一个简单的示例代码，使用 PyRay 库实现一个简单的光线追踪： ```python import pyray def main(): # 创建场景 scene = pyray.Scene() # 定义光源 light = pyray.Light(position=(1, 1, 1), intensity=1) # 定义物体 sphere = pyray.Sphere(center=(0, 0, 0), radius=1, color=(255, 0, 0)) # 将物体添加到场景中 scene.add_object(sphere) # 设置光源 scene.set_light(light) # 渲染图像 image = scene.render(width=800, height=600) # 保存图像 image.save("output.png") if __name__ == "__main__": main() ``` 这个示例中创建了一个简单的场景，含一个红色的球体和一个光源。通过调用 `scene.render()` 方法，可以渲染出一个指定大小的图像，并保存为文件。 需要注意的是，光线追踪是一种计算密集型的算法，对于复杂场景或高分辨率图像可能需要较长的渲染时间。此外，还可以通过调整光线追踪的参数和使用更高级的技术来提高渲染效果和速度。

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值